Impacts of brine disposal from water desalination plants on the physical environment in the Persian/Arabian Gulf

The numerical simulations were run with a high-resolution (Δx = Δy = 1/36-degree or ≈2.8 km) implementation of the Hybrid Coordinate Ocean Model (HYCOM) (Bleck 2002, Halliwell 2014). The experiments considered the geographic domain delimited by latitudes 9.3°S–30.6°N and longitudes 33.0°E–76.9°E (figure 1). The vertical structure was discretized in 22 hybrid layers, all of them allowed to transform from isopycnic to terrain-following or to z-coordinates. The bathymetry was based on the NOAA-NGDC's ETOPO1 dataset (Amane and Eakins 2009). The shallowest depths were kept in 2 meters (regions with depths less than 2 meters were considered land). A one-way nesting scheme was used to provide the initial and boundary conditions. This approach, similar to the one used by (Campos et al 2020), is summarized as follows. First, a global 1/12-deg, 32-layers experiment with HYCOM was initialized with results downloaded from the HYCOM Consortium and forced for 21 years with products from the Comprehensive Ocean and Atmosphere Data Set (COADS) (Woodruff et al 1987). For lack of data storage space, a decision was made to save only the averages every six days of the global run products. Following, the available six-days averages of all relevant fields from the global run, starting on day 355 of the climatological year 21, were remapped onto the 22 isopycnic surfaces and interpolated to the 1/36-deg horizontal grid used in the present work. Then, the experiments with the nested model were executed using the interpolated global products as initial and boundary conditions, also forced with COADS climatological products. The COADS forcing fields used consisted of surface air temperature, net downward radiation flux, net downward shortwave flux, precipitation, vapor-mixing ratio, surface wind speed and wind stress. Similarly as in (Campos et al 2020), at each time step, the products of the model's integration were relaxed to the boundary conditions, in a 36 grid-points 'buffer' zone with an e-folding time of 53 days.

Zoom In Zoom Out Reset image size

A pair of twin-runs were performed for the entire oceanic region shown in the bottom panel (a) of figure 1, including the Arabian Sea, the Gulf and the Red Sea. In the first run, hereinafter referred as NODESAL, no anthropogenic salt input was considered. In the second, hereinafter referred as DESAL, selected desalination plants were included around the Gulf in the form of 'negative rivers' or 'rivers of salt', mainly along the Arabian Peninsula coastline. Both experiments included the most significant rivers in the modeled region. For more details on the rivers and desalination plants, see table 2. The climatological values for the riverine inflows (here, 'inflow/outflow' indicates water into/from the Gulf), in cubic meters per second (m³ s⁻¹), were obtained from a global database of monthly mean river discharge constructed at the U.S. Naval Research Laboratory (Barron and Smedstad 2002). In the entire area considered in the model, nine rivers with mean flow larger 128m³s⁻¹ were considered. From these, only three discharge their waters directly inside the Gulf: The Shatt-al-Arab, the Karun and Karkkeh. The combined mean inflow of these three rivers is 2,900m³ s⁻¹ on average, according to table 2. At each of the desalination plants, due to the lack of more accurate information in the literature, the freshwater outflow was prescribed in a somewhat arbitrary way, in general, reflecting the order of magnitude of the values listed in table 1. The total mean value of the ten plants listed in table 2 is of 120 m³ s⁻¹, which is higher than the total value of approximately 80 m³ s⁻¹ from the Gulf's twelve plants listed in table 1. It is likely that the actual values are either higher than those disclosed by the plants or will be in a short while.

In both experiments, tidal forcing was not considered in the simulations. Findings of previous studies show that, despite their relevance on the shorter term, tides do not generate significant residual currents in the longer time scales in the Gulf [e.g.:] (Pous et al 2012, 2015). However, this does not mean that residual currents resulting from non-linear interactions of the tidal flow with topography should be completely ignored. There are regions in the Gulf where the residual currents can exceed 0.15 m/s, or about 10% of the local barotropic velocity (Poul et al 2016). The effects of tides will be addressed in future experiments.

The products of a 6-years numerical integration, from Jan-01-0022 to Jan-01-0027, were saved daily, for both the internal (baroclinic) and external (barotropic) modes at each point of the entire grid.

One of the objectives of the present work is to assess the significance of any salinity increase associated with the extraction of freshwater at locations along the Gulf's shores. For doing that, a methodology to calculate the salt budget and any resulting trend was devised considering a semi-enclosed sea as represented in top right panel (c) of figure 1. The equations for conservation of volume (V) and total salt content (S) within the region can be written as follows.

$$\begin{eqnarray}&&{V}_{t}=\iint {v}_{H}{dxdz}+{EMP}+D+{Riv}+{{Rest}}_{V}\end{eqnarray} \tag{ 1 }$$

and

$$\begin{eqnarray}&&\iiint {S}_{t}{dV}+\overline{S}{V}_{t}=\iint {v}_{H}{S}_{H}{dxdz}+{SEMP}+{SD}+{{Rest}}_{S},\end{eqnarray} \tag{ 2 }$$

In these equations, V_t  = ∂V/∂t is the volume trend; S_t  = ∂S/∂t is the total salinity trend; $\overline{S}=1/V\iiint {SdV}$ is the volume-averaged salinity; EMP is the evaporation minus precipitation over the Gulf and D is the total volume of water extracted by the desalination plants. SEMP represents the net flux of salt associated with evaporation minus precipitation and SD is the total salt influx resulting from the desalination process. Rest_V and Rest_S are residual terms due to small scale mixing and truncation errors, v is the meridional component of the velocity and the subscript H indicates the variable at a vertical section across the Strait of Hormuz (see figure 1). Note that volume trend term (V_t ) is the link between the two equations.

As it is more intuitive to deal with volumes of water being extracted from or discharged into the sea, the salt budget equation (equation (2)) is transformed into an equivalent-freshwater (hereinafter simply freshwater) budget equation, dividing it by a reference salinity S₀:

$$\begin{eqnarray}&&{M}_{t}+\displaystyle \frac{1}{{S}_{0}}\overline{S}{V}_{t}=\displaystyle \frac{1}{{S}_{0}}\left(\iint {v}_{H}{S}_{H}{dxdz}+{SEMP}+{SD}+{{Rest}}_{S}\right),\end{eqnarray} \tag{ 3 }$$

where

$$\begin{eqnarray*}&&{M}_{t}=\displaystyle \frac{1}{{S}_{0}}\iiint {S}_{t}{dV}.\end{eqnarray*}$$

With the insertion of V_t from equation (1) into equation (3), the total freshwater trend, M_t , can be estimated by the equation:

$$\begin{eqnarray}&&{M}_{t}=\displaystyle \frac{1}{{S}_{0}}\left(\iint ({S}_{H}-\overline{S}){v}_{H}{dxdz}+({SEMP}-\overline{S}{EMP})-(\overline{S}D-{SD})\right)+\displaystyle \frac{1}{{S}_{0}}\left(\overline{S}{Riv}+({{Rest}}_{S}-\overline{S}{{Rest}}_{V})\right)\end{eqnarray} \tag{ 4 }$$

Here,

$$\begin{eqnarray*}&&{EMP}=\iint (E-P){dxdy},\end{eqnarray*}$$

where E—P (Evaporation minus Precipitation) is an output of the numerical model;

$$\begin{eqnarray*}&&D=\displaystyle \sum _{i=1}^{N}{D}_{i},\end{eqnarray*}$$

N is the number of desalination plants considered and D_i the volume processed by each plant;

$$\begin{eqnarray*}&&{SEMP}=\iint ({SS}(E-P)){dxdy},\end{eqnarray*}$$

SS is the surface salinity and

$$\begin{eqnarray*}&&{SD}=\displaystyle \sum _{i=1}^{N}{S}_{i}{D}_{i}.\end{eqnarray*}$$

Note that SRiv = 0 because the salinity S at the river mouth is equal to zero.

It is important to note that the terms multiplied by $\overline{S}$ in the above equations represent the contributions to the trend from the changes in volume. Although small, V_t is not necessarily equal to zero since there are changes in SSH. In addition, it is assumed that the smaller scale effects and truncation errors in the twin runs are equivalent (perhaps this is not entirely true but the differences between these terms are expected to by small enough to be neglected). Thus, the difference between the results of equation (4) for the two experiments (DESAL—NODESAL) eliminates the residual terms and yields an equation for computing ${\rm{\Delta }}{M}_{t}^{* }$, the freshwater loss due to the desalination process within the modeled region:

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\Delta }}{M}_{t}^{* } & = & \displaystyle \frac{1}{{S}_{0}}\left(\iint ({\rm{\Delta }}{S}_{H}{v}_{H}-{\rm{\Delta }}\overline{S}{v}_{H}){dxdz}\right)\\ & & +\displaystyle \frac{1}{{S}_{0}}\left({\rm{\Delta }}{SEMP}-\overline{S}{\rm{\Delta }}{EMP}+{\rm{\Delta }}\overline{S}{Riv}-({\rm{\Delta }}{SD}-{\rm{\Delta }}\overline{S}D)\right)\end{array}\end{eqnarray} \tag{ 5 }$$

Each experiment was started on the 1st of January of climatological year 22 and let to run for five years, to 1st of January of year 27. The time-history of each realization is represented by the time series of the basin-averaged sea surface height (SSH) and kinetic energy (KE), plotted in figure 2. The plots show that it takes about one year for the model to reach a state of reasonable equilibrium in both experiments. After a pronounced spike in KE in the first year of the run, due to the initialization process, the curves show a more uniform behavior, with strong seasonality, certainly associated with the forcing seasonal variability. They also suggest some energy in higher and lower frequencies. The intraseasonal variability could be associated with either the forcing or nonlinear mesoscale features. However, since the model is forced with climatological products, with no year-to-year changes, the interannual variability displayed in the KE and SSH curves are certainly associated only with the model's internal variability. It is also relevant to point out that, despite being small, there are some noticeable differences between NODESAL and DESAL, especially in KE.

Zoom In Zoom Out Reset image size

Because the main objective of this work is to investigate the differences between the two experiments, only a brief description of the model's general results is given here. As suggested by the KE and SSH curves in figure 2, in both experiments, with and without desalination, the circulation shows a marked seasonal variability. The spatial structure has a dominant two-layer structure, in good agreement with the general pattern described in the literature [e.g.:] (Reynolds 1993, Johns et al 2003, Campos et al 2020). Similarly to a previous work based on results of experiments with HYCOM (Campos et al 2020), in the upper-layer, defined by the model's top eight isopycnic surfaces, the horizontal circulation is dominated by a year-round cyclonic gyre, with less saline waters from the Gulf of Oman entering the Gulf through the Strait of Hormuz forming two branches. The northern branch flows along the Iranian coast, reaching the northwestern regions of the Gulf. The southern branch flows towards the coasts of the UAE and is the most affected by the seasonality, being stronger in the summer. In the winter, the top layer circulation has the two branches well defined but weaker mesoscale activity. Starting in April, strong meso-scale eddies start to form and remain until December, with a sequence of quasi-stationary cyclonic and anticyclonic vortices along the Gulf's deeper channel. In the bottom layer, the overall salinity is higher and the predominant direction of the flow is towards the Strait of Hormuz, in agreement with the inverse-estuary circulation patter described in previous works (Reynolds 1993, Johns et al 2003, Yao and Johns 2010, Campos et al 2020). This two-layer structure is well represented by the vertical structure of temperature, salinity and velocity on a section along 26°N, the 'Hormuz section', indicated in figure 1(b). As seen in figure 3, which represents the time averaged distributions, the velocity is positive (towards north) in the upper layer and negative (southward) in the bottom layer. The mean temperature decreases with depth while the salinity is lower in the region above the interface defined by the isopycnic layer k = 8 (σ_t  = 24.70 kg m⁻³) and has higher values below. It should be pointed out that, despite the fact that the salinity distribution shows a region of lower salinity near the easter side of the vertical section, the density (panel c) increases monotonically with depth, indicating stability of the water column.

Zoom In Zoom Out Reset image size

To assess the impact of the brine injection on the Gulf's residual circulation, differences between the two experiments (DESAL—NODESAL) for the seasonal means (Summer and Winter) and long term mean were taken. Figure 4 displays these differences for salinity (left panels) and velocity and SHH (right panels). In that figure, JAS and JFM seasons are, respectively, the averages for the months Jul-Aug-Sep and Jan-Feb-Mar of the last year of the run (Year 26). The 'Mean' is the time average for the period Jan-1-0023 to Jan-1-0027. In spite of the fact that the only difference between the two experiments was the extraction of freshwater in a few points, after five years there are relatively large differences in salinity, SSH and velocity, especially along the deeper channel of the Gulf, near the Iranian coast. These differences are prominent during JAS and appear closely related with the mesoscale activity that are more intense during the summer season.

Zoom In Zoom Out Reset image size

To some extent, the large differences shown in figure 4 represent an intriguing result. To better understand these impacts of the brine injection (or freshwater extraction), an animated plot (figure S1 (available online at stacks.iop.org/ERC/2/125003/mmedia), Supplementary Material) was produced with a high-frequency sequence of snapshots of the difference in SSH (DESAL—NODESAL) during the first couple of days of run. The 'movie' shows that the perturbation of the sea level caused by the injection of salt at some points along the Gulf shores generates anomalies that propagate as surface gravity waves, reaching regions as far as the northern end of the Red Sea in less than 24 h. These SSH anomalies then generate disturbances in other properties, as in salinity, for instance, according to the physical process described as follows.

Consider the equation for conservation of salinity:

$$\begin{eqnarray}&&\displaystyle \frac{\partial S}{\partial t}=-\left(u\displaystyle \frac{\partial S}{\partial x}+v\displaystyle \frac{\partial S}{\partial y}\right)+{K}_{H}\left(\displaystyle \frac{{\partial }^{2}S}{\partial {x}^{2}}+\displaystyle \frac{{\partial }^{2}S}{\partial {y}^{2}}\right)+\displaystyle \frac{\partial }{\partial z}\left({K}_{z}\displaystyle \frac{\partial S}{\partial z}\right),\end{eqnarray} \tag{ 6 }$$

where K_H and K_z are the horizontal and vertical components the turbulent diffusivity coefficient.

The integration of equation (6) from the bottom of an isopycnic layer (−h) to the sea surface (η) yields:

$$\begin{eqnarray}&&{\int }_{-h}^{\eta }\displaystyle \frac{\partial S}{\partial t}{dz}=-{\int }_{-h}^{\eta }\left(u\displaystyle \frac{\partial S}{\partial x}+v\displaystyle \frac{\partial S}{\partial y}\right){dz}+{\int }_{-h}^{\eta }{K}_{H}\left(\displaystyle \frac{{\partial }^{2}S}{\partial {x}^{2}}+\displaystyle \frac{{\partial }^{2}S}{\partial {y}^{2}}\right){dz}+{F}_{\eta }-{F}_{-h},\end{eqnarray} \tag{ 7 }$$

where F_η and F_−h are, respectively, the salt fluxes across the free surface and the bottom of the layer.

The zonal and meridional components of the geostrophic velocity can be written as:

$$\begin{eqnarray}&&{u}_{g}=-\displaystyle \frac{g}{f}\displaystyle \frac{\partial \eta }{\partial y}-\displaystyle \frac{1}{\rho f}\displaystyle \frac{\partial }{\partial y}{\int }_{-h}^{\eta }\rho (z){dz}\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&{v}_{g}=\displaystyle \frac{g}{f}\displaystyle \frac{\partial \eta }{\partial s}-\displaystyle \frac{1}{\rho f}\displaystyle \frac{\partial }{\partial x}{\int }_{-h}^{\eta }\rho (z){dz}\end{eqnarray} \tag{ 9 }$$

These two equations, equations (8) and (9), show that the geotrophic flow, the predominant component of the velocity in the meso-to-large scale dynamics, depends on the sea surface elevation η. On the other hand, the release of brine at a point changes the surface elevation, generating surface gravity waves that propagate with speed $c=\sqrt{{gH}}$. For instance, for H = 10 meters, c ≈ 10 m s⁻¹ or 864 km/day. This means that SSH anomalies originating in the Gulf would propagate with speeds of about a thousand kilometers per day, or faster. Consequently, according to equation (7), other thermodynamic properties would be affected as well. Contrary to what one would expect, instead attenuating with time, the perturbations excite normal modes of variability, through some sort of resonance, or trigger mesoscale activities by means hydrodynamic instability. Although no further investigation was done to test these hypotheses, the fact is that, instead of dissipating, these small initial SSH anomalies ended up in the sustained and significant perturbations in the Gulf's salinity, SSH and velocity distributions shown in figure 4. The differences are much well pronounced during the summer months, July trough September (JAS), with relatively strong anomalies associated with the centers of cyclonic and anticyclonic circulation along the axis of the Gulf's deeper channel. During JAS, the salinity anomalies vary from near −0.60 to more than 0.60 psu. The range in SSH differences is from −10.0 to 10.0 cm while velocity can be as higher as 0.20 m s⁻¹ over the deep channel. These differences are much weaker during the winter (JFM), especially in SSH and velocity. Although small, when compared with the values in JAS, the differences in the long-term means are not negligible.

The bottom left panel of figure 4 show that the difference (DESAL—NODESAL) in the mean surface salinity has a noticeable spatial variability. In particular, the areas surrounding Bahrain, Qatar and the northern coast of the UAE show more intense red colors than in other regions. This means that those areas may be more affected by the salt buildup resulting from the brine injection. To further investigate this possibility, the area averaged differences in the surface temperature and salinity distributions from the two experiments (DESAL—NODESAL) were plotted for six areas around of the desalination plants, as indicated in figure 5. For each of the areas labeled as Kuwait, KSA1/Bahrain, Qatar, UAE3, UAE1 and Iran, the 5-years mean and standard deviation were calculated, as listed in table 3. The results show that, in general, after five years of discharging brine in those specific areas, both mean changes of salt and temperature are very small. The highest values for the mean salinity, in the order of 0.1 psu are found in the KSA1/Bahrain and Qatar areas. Both of these locations show relatively small standard deviation. On the other hand, in the UAE1 region, in spite of the smaller mean salinity, in the order of 0.04 psu, the standard deviation is the highest, with a value of the order of 0.2 psu. In UAE1, there is also a relatively larger variability in the mean temperature, of about .016 °C as compared with the other sites, except for Kuwait. At this point one could argue that, considering that there is a constant input of salt, the mean values calculated for the entire period are not representative of the actual values at the end of the simulation. However, as clearly seen in the plots of figure 5, except for UAE1 and Iran, there are no noticeable trends in the time series. In the UAE1 area, considering only the last 24 months of the run, the mean temperature and salinity values are of the order of 0.2 °C and 0.2 psu, respectively. However, the plots for the UAE1 area also show considerable interannual variability. This raises a question on the statistical significance of the trend and suggests that the mean for only the last two years could be biased.

Zoom In Zoom Out Reset image size

As mentioned in section 1, the Gulf is not an enclosed sea. By means of an inverse-estuary circulation, it exchanges water with the Indian Ocean and, to compensate the volume decrease due to the extraction of freshwater by the desalination process, water must be pumped in through the Strait of Hormuz. However, this is not a simple and straightforward dynamic process. The small changes due to brine injection into the Gulf lead to changes in the overall temperature and salinity, which lead to changes in EMP and in the residual circulation. The net result is that, in spite of maintaining the volume constant, the basin wide mean salinity would be different in unpredictable ways. To evaluate the changes resulting from the introduction of some desalination plants, as described in the previous sections, equation (5) was used to calculate the trend in the total volume of equivalent-freshwater in the Gulf. For that, the constant value of 39.64 psu was adopted for the reference salinity S₀. The results are summarized in figure 6.

Zoom In Zoom Out Reset image size

On the left panels of figure 6, the individual terms on the right-hand-side of equation (5) are plotted: (a) the residual freshwater flux across 26°N; (b): the difference in EMP from the two experiments; (c): the net effect of the rivers discharge; and (d): the effect of brine injection. All values are in million cubic meters per day. On the right, the top panel (e) shows M^* _t , the sum of all RHS terms in equation (5), and in the lower panel (f) is the time integral of ${M}_{t}^{* }$, which represents the total freshwater loss due to the desalination process during the five years of the experiment. These results show that, when the changes associated with the volume trend (V_t ) are considered, the residual impacts of river inflow, EMP and the desalination plants outflow are relatively small, as compared with the lateral exchange at the Strait of Hormuz. The combined effect is a negative trend of 20.2 million cubic meters per day, resulting in a total loss of 36.4 km³ in five years. Considering the Gulf's mean volume of 8118 km³ (area of 2.39E+5 km² and mean depth of 33.96m), this freshwater loss corresponds to a salinity increase of 0.18 psu in five years. This value is similar to the result obtained by (Lee and Kaihatu 2018), but only half of the salinity increase reported by (Ibrahim and Eltahir 2019). However, one must consider that the present study is based on only five years, as compared with the 12 years in the cited works (Lee and Kaihatu 2018, Ibrahim and Eltahir 2019).