Effect of reservoir size on water quality in coastal reservoirs during desalinization period using numerical model

After the industrial revolution and expansion of urbanization, the growing demand for water in the world was met by constructing large-scale surface and groundwater reservoirs. However, looking at the current state of the water crisis, it could be concluded that with construction of these structures, there is still a water shortage crisis. Research results show that by 2050, the number of countries with water shortages will increase to 54, affecting 4 billion of the estimated global population (Yang & Lin 2011). Therefore, it can be concluded that water scarcity is a major economic and social problem today and will soon become one of the main constraints for future economic development. Alternatively, looking at the process of dam construction in the world, it is clear that almost all large surface water storage reservoirs have already been built in suitable sites. In addition, few large dams have been built in the world since 2010, because of the negative impacts on the ecosystem (Yang 2018). Conversely, according to Herrera-León et al. (2018), in the process of economic and demographic growth, modern metropolises are gradually appearing in coastal areas. Rapid economic development in these areas has intensified the problem of water shortage. This situation is becoming more complicated by the irregular spatial and temporal distribution of freshwater resources in these areas.

By passing over the traditional building methods of dams and constructing surface reservoirs, it is necessary to go to new and more reliable methods of water supply in order to overcome the problems of water shortage in the present and near future. One of the most crucial methods is construction of coastal reservoirs to collect depleted floodwaters from rivers to the sea. A coastal reservoir consists of an impermeable sea wall (earthen or concrete) to control fresh water and prevent it from mixing with seawater (Yang & Ferguson 2010). These reservoirs are built at the river estuary (first generation) or on the adjacent shoreline (second generation) (Yang & Liu 2010; Kolathayar et al. 2019). Sivakumar et al. (2020) stated that this technology can convert floodwater into valuable water resources closer to the coastal population centres where more than 50% of global population live. In the world, the runoff lost to the sea is approximated as 45,500 km³/year, whereas the 7 billion people's drinking water is only 380 km³/year, or 0.8% of the runoff.

Many coastal reservoirs have been built in China, South Korea, Hong Kong and Singapore (Yuan et al. 2007; Liu et al. 2015). It seems that the worst enemy for a coastal reservoir is the water quality of the source water. With poor maintenance management, these reservoirs may become wastewater reservoirs, such as the Sihwa Lake in South Korea (Bae et al. 2010).

In these reservoirs, inlet and outlet water flow are adjusted using different gates. In this way, only an acceptable quality of water enters the reservoir for use during water shortage periods. Figure 1 shows a view of Plover Cove reservoir in Hong Kong.

Generally, limited studies have been carried out on the performance of coastal reservoirs. Yang & Lin (2011) introduced a new technology of coastal reservoir construction using flexible dams. They mentioned that the main problem of conventional coastal reservoirs is the reduction of water quality during the operation period due to infiltration of saline currents, which can be solved with this technology.

Yang et al. (2013) pointed out that different coastal reservoirs can be constructed in series and therefore, based on the quality of water output from each reservoir, different needs such as domestic water, agriculture and industry can be met. Sitharam (2017) conducted an economic analysis and showed that the construction costs and water price of a coastal reservoir are much lower than for other methods. Yang (2018) estimated the lifespan of the dam to be 150 years (considering the sedimentation of the reservoir and the useful lifespan of the construction materials) and concluded that by 2150 there will be no reservoir dam in the world for drinking water supply and introduced coastal reservoirs as a new solution to overcome this problem. Kolathayar et al. (2019) pointed out the pros and cons of coastal reservoirs compared to other water supply solutions. They showed that coastal reservoirs are more economical than other methods of water supply. They also pointed out that the most important challenges for operation of these reservoirs is to control salinization/seawater intrusion, which can be done with proper management of the reservoir. Rakhmawati & Armono (2021) analyzed the socio-economic of creating coastal reservoir in Welang Coastal Area, Indonesia and compared it with on-land reservoirs. Their analysis showed that the Benefit Cost Ratio (BCR) value for on-land reservoirs is 1 and for coastal reservoirs is 1.54. In addition, they claimed that with the use of the coastal reservoirs, the river runoff, water supply and the water quality in Welang river basins can be improved. Ji et al. (2022) discussed how to apply a coastal reservoir at the Shoalhaven River mouth to secure an additional water supply for ever-growing Greater Sydney. They concluded that the proposed reservoir could provide sufficient water to meet the increasing water demand caused by population growth with a reliability of 90% per year.

Limited studies have been carried out on the use of numerical methods to model the flow and quality of water in coastal reservoirs. Mao et al. (2005) simulated the desalination process of water in the Polder reservoir using pumping operations near the bed of the reservoir. They used the Delft3D numerical model in their study. The results of this study showed that this method of reducing salinity has good efficiency so that after three years of simulation in the hydrological conditions of the region, the salinity of the reservoir water was reduced by 60%. Chen (2014) used a 3D numerical model to analyze the residence time of saltwater intrusion in the Qingcaosha reservoir in the Yangtze River estuary. Results showed that the flow in the reservoir was mainly controlled by the wind-driven current and, to a lesser extent by the inlet/outlet current. Liu (2017) investigated hydrodynamic and numerical salinity simulations in the Lower Lakes through setting up 1D and 2D models by using MIKE software. A 1D model is applied for five barrage structures while a 2D model was used to reproduce the hydrodynamic processes and salinity changes in the Lower Lakes. He proposed two coastal reservoir designs in the numerical model. He concluded that the coastal reservoir could provide 150 GL/year water under three typical low to moderate flow conditions and the salinity in the reservoir was around 0.3 PSU. A report from Falconer et al. (2020) describes a general 3D numerical modelling for predicting hydrodynamic, water quality, sediment transport and morphological processes. This is followed by two coastal reservoir examples, including: (1) the Cardiff Bay project where a tidally flushed bay was closed off from the outside estuary to create a freshwater reservoir for urban regeneration and (2) the West Somerset Lagoon for tidal range power generation. Karimaei & Fouladfar (2021) used MIKE3 numerical model to simulate the flow and salinity of an understudied coastal reservoir along the Caspian Sea near Tajan River estuary in Iran. For the desalination process of the reservoir, the average monthly discharge of the river for a period of one year (first model) and a flood hydrograph with a return period of about 1,000 years (second model) were used. The results showed that, despite the first model, intense density stratification occurs when flood water hydrograph enters the reservoir. In addition, at the end of the simulation period, the maximum salinity in the whole reservoir decreased from the initial value of 12 g/l (which is related to the reference salinity of the Caspian Sea) to less than 2 g/l in the first model and 1 g/l in the second model.

Unfortunately, there are many uncertainties about the performance of coastal reservoirs and research in this regard still needs to go far way. The present study tried to investigate the effect of coastal reservoir size on reservoir water salinity, temperature and density during desalinization period.

With 890 km coastlines in the north (Caspian Sea) and 4,910 km in the south (Oman Sea and Persian Gulf) of Iran, this country has a distinguished geographical position in Asia. The existence of these long coastlines can provide suitable economic opportunities for seven coastal provinces and the whole parts of country.

However, there is a backlog mainly due to lack of infrastructure such as access to drinking water and sanitation in these areas, especially the southern coast of the country (Naseri & Najafi 2019). Conversely, a large part of the runoff in these areas enters the sea without any use. Studies have shown that in the Persian Gulf and Oman Sea watershed with a drainage area of 424,515 km² and an average annual rainfall of 134 mm in Chabahar to 1,457 mm in Koohrang, about 4.5 − 6 × 10⁶ m³ of runoff is discharged annually to the sea. Similarly, in the Caspian Sea watershed with a drainage area of 58,253 km² and a range of rainfall changes from a minimum of 260.5 mm in Bojnourd to a maximum of 1,863 mm in Anzali port, about 6.8 × 10⁶ m³ of runoff is discharged into the sea (Karimaei & Fouladfar 2022). In this way, due to lack of effective control of runoff, significant volumes of fresh water with good quality flow into the sea each year. However, almost all sites with ideal hydro-social and geotechnical combinations have already been dammed.

According to above issues, the construction of coastal reservoirs in the coastlines of the country can be reasonable and therefore the problem of water shortage in water crisis areas of the country can be solved to some extent. It is also possible to reduce treatment costs of drinking water by treating the water in these reservoirs. Alternatively, it seems that in some southern parts of the country, use of these reservoirs is the only solution to overcome the water crisis. In some areas where the problems of water shortage are so serious, efforts are being made to access deep water in the ground with a huge cost and low efficiency.

Tajan River is one of the most important rivers in the Caspian Sea basin, it rises from a height of nearly 3,250 m on the northern slopes of the Alborz Mountain range in the south part of Sari city, the capital of Mazandaran province, in Iran. Figure 2 shows a satellite image of this river estuary. Flow in this river is influenced by some upstream hydraulic structures, such as the Shahid Rajaee dam. In present study, the coastal reservoir bounds are located near the outlet of Tajan River (Figure 2(c)). The dimensions of the coastal reservoir in this study are different. Its width is equal to 1 and 2 km and its length (along the coastline) is equal to 4.5 and 9 km. Water in the reservoir is supplied from the river through an access flood channel with 20 m wide.

In May 2019, due to heavy rains in the upstream basins, a relatively large flood occurred in this river. Measurements showed that the flood peak discharge was equal to 180 m³/s and the volume of flood due to the shape of the hydrograph was more than 40 MCM at Kordkhil station. Figure 3 shows the recorded hydrograph flow at Kordkhiel station. Studies have shown that the return period of this flood is about 5 years (Regional Water Company of Mazandaran 2021).

Studies have shown that the maximum amount of total dissolved solids at Kordkhil station (near river estuary) is less than 500 mg/l (Karimaei & Fouladfar 2021). As an assumption, the salinity refers to all the dissolved solids amount in the water and this measurement was considered to be 500 mg/l or 0.5 PSU¹ in present study.

In this study, the MIKE software package developed by Danish Hydraulic Institute (DHI) was employed to conduct a series of numerical simulations. This suite of software applications is a reliable and common simulation in the field of sea and coast hydrodynamics. In the present study, MIKE3 which is the 3D module of this package was used. The reason for choosing this model is to simulate the water salinity, temperature and density gradient within the reservoir more accurately. MIKE3 flow software is based on a flexible mesh approach. The system is based on the numerical solution of the 3D incompressible Reynolds averaged Navier–Stokes equations invoking the assumptions of the Boussinesq approximation and hydrostatic pressure.

In present study, the developed model consisted of continuity, momentum, temperature, salinity and density equations and it is closed by a turbulent closure scheme. The spatial discretization of these equations was performed using unstructured triangular meshes in the horizontal plane and a cell-centred finite volume method. The spatial domain was discretized by subdivision of the continuum into non-overlapping element/cells.

Dimensions of coastal reservoirs depends on various factors such as water resources planning, water quality control, amount of inflow and etc. In this work, in order to study the effect of reservoir size on hydraulic and water salinity parameters during desalinization period, four types of reservoirs (models) with different dimensions were considered. The first model with surface area dimensions of 1 km × 4.5 km (R1 × 4.5), second model with dimensions of 2 km × 4.5 km (R2 × 4.5), third model with dimensions of 1 km × 9 km (R1 × 9) and fourth model with dimensions of 2 km × 9 km (R2 × 9) were evaluated. Boundaries and initial conditions in these models were similar, and the model outputs were flow (depth and velocity), and water quality parameters (water salinity, temperature, and density).

Figure 4 shows the general layout of the reservoir and associated various components. In this model, a desilting basin with 350 m width was considered. The length of the desilting basin was such that its connection with the main reservoir was established in a length of 200 m in the deep part of the main reservoir. In addition, in the eastern part of the reservoir in all models, an emergency spillway (ogee crest) with a length of 90 m (1 and 2% of the reservoir wall in different models) was introduced and the difference between the level of its crest and sea level (initial reservoir level) was considered as 0.5 m. To build the numerical models, first, hydrographical maps of the region were used to determine the bathymetry (depth changes) of the reservoir. Second, a numerical mesh was generated on the study area. Figure 5 shows the depth variation in the reservoir and performed triangular meshes. In this figure, the horizontal and vertical axes represent latitude and longitude coordinates, respectively. In order to achieve higher accuracy, meshes with different sizes at the inlet, the middle parts of the reservoir and near the spillway were determined.

At the inlet of the reservoir and near spillway, due to presence of more intense currents, smaller mesh size was considered rather than other regions. To analyze the sensitivity of the computational mesh size, the salinity of the reservoir at a midpoint of the reservoir surface was examined. As a sample, in the second model, the maximum mesh area was considered to be 5,000, 10,000 and 15,000 m² and salinity of the point in each model was calculated and is shown in Figure 6. As shown in this figure, the selected mesh sizes have little effects on salinity changes. As a result, the largest mesh area in the final models was limited to 10,000 m² at the middle part of the reservoir and 2,000 m² near the spillway and access channel. In addition, the results for different size of reservoirs showed that only the number of computational meshes should be increased for larger reservoirs. Furthermore, a quadrilateral type of mesh was used for mesh generation of the access channel. Finally, to investigate the change of water salinity in depth of the reservoir, five computational layers were determined. Further studies showed that as the reservoirs were shallow (with maximum depth of 12 m), increasing the number of layers had little effect on the results.

After mesh generation, three different boundary conditions were introduced to the models. These boundaries included: (1) the river boundary at the junction of the access waterway from the Tajan River to the reservoir (Figure 3). In addition, the salinity of river water was considered as 0.5 PSU. The river water temperature was also determined according to the European Meteorological Organization data. (2) The sea boundary is around the terminal structure of the spillway. As the Caspian Sea is not affected by tidal phenomenon, the sea level was considered as zero. (3) The land boundary around the reservoir, which was introduced to the models as a zero-vertical velocity.

Data from the European Meteorological Organization were used to meet precipitation and evaporation data in the region. These parameters were considered constant spatially but variable over time, which is a good approximation due to the small sizes of the reservoirs. For heat exchange of the reservoir, there was a need for climatic information of the region including air and river temperature, clearance coefficient and relative humidity. A clearance coefficient of 100% specifies a clear sky and 0% specifies cloudy weather. These data were extracted from satellite and European Centre for Medium-Range Weather Forecasts (ECMWF) data sets. The Mike3 software uses this information to calculate the exchange of reservoir temperature with the air and ultimately simulate water temperature in the reservoir. The initial water surface level and its salinity in the reservoir was similar to that of the Caspian Sea. The latter was equal to 12 PSU.

The simulation period was considered as 8 days from January 1, 2013, according to the duration time of the inflow hydrograph (Figure 3). Furthermore, with time step of 30 seconds, the total number of simulation steps was 23,040. The maximum running time with a system of CORE ™ i7–7,500 U and 8GB RAM Intel was about 4 hours. In this section, the analysis and discussions of the model results, including flow (water surface level and velocity), salinity, temperature and density parameters is presented separately.

Study on the water level showed that the value of this parameter was a function of reservoir size and simulation time. Figure 7 shows the temporal variations of water level in different coastal reservoirs at a midpoint of the reservoir. The position of this point is shown in the figure. From this figure, the maximum increase of water level in the smallest reservoir (R1 × 4.5) is about 2 m and in the largest reservoir (R2 × 9) is about 1.4 m. This can affect the upstream flooding region of the beach. In addition, the difference in water level in two reservoirs of R1 × 9 and R2 × 4.5 with similar surface areas is insignificant, although the increase in water level of R1 × 9 is slightly more (about 5%). The increase in water level was gradually discharged through the spillway. Results also showed that flow velocity gradually decreased from the sediment trap reservoir (desilting basin) to the eastern wall near the spillway. The maximum velocity was calculated as about 0.1 m/s around the middle part of the reservoir R1 × 4.5.

Salinity is the most important parameter for water quality in coastal reservoirs during the desalinization period. The lower the salinity of the water, the lower the cost of subsequent treatment. Gradually, with entry of fresh water into the reservoir, the salinity of the water decreased obviously the larger the reservoir, the longer the time for decreasing the salinity of water in the reservoir.

Figure 8 shows the salinity changes in the surface water after the passage of hydrograph peak discharge for the two reservoirs of R1 × 9 and R2 × 4.5 with the same surface area. In contrast with reservoir R2 × 4.5, the salinity of water in the desilting basin of reservoir R1 × 9 was similar to the salinity of the inflow (0.5 PSU). In addition, Figure 9 shows the salinity changes in the two reservoirs at the end of the simulation. This figure shows the better performance of the shallower reservoir R1 × 9 in reducing the salinity of the reservoir after the desalination process. The salinity of almost all parts of the reservoir R1 × 9 was less than 4 PSU, while the salinity of different points in reservoir R2 × 4.5 was about 6 PSU. Therefore, it seems that the construction of shallow coastal reservoirs can be more efficient during the desalination period, however, more research is needed. It should be noted that the evaporation phenomenon is a serious restriction for finalizing the reservoir depth.

To investigate this issue better, the temporal variation of salinity at two points on the surfaces of the reservoirs are shown in Figure 10. The coordinates of these points are also shown in the figure. This figure shows that the size of the reservoir is a crucial parameter in the process of reducing the salinity of the reservoir. Furthermore, no meaningful difference was observed between salinity of two points in a specific reservoir at the end of the simulation. According to this figure, the highest salinity at the end of the simulation is related to the reservoir R2 × 9 with the value of about 8 PSU and the smallest value is related to the reservoir R1 × 4 with about 2 PSU.

Next, Figure 11 shows the salinity profile along a transverse line in the middle of the reservoir after passing the hydrograph peak discharge. Generally, salinity stratification is obvious in all modelled reservoirs. However, due to faster layer mixing in shallow reservoirs, the stratification is more intense. Therefore, the salinity of a large part of the cross-sectional area in reservoir R2 × 4.5 (Figure 11(c)) and reservoir R2 × 9 (Figure 11(d)) remained unchanged and equal to 12 PSU.

Finally, Figure 12 shows the salinity variations in different layers of the reservoir. The position of the point is shown in Figure 10(b). This figure shows that the salinity difference in the upper and lower layers is mainly a function of the inflow discharge, so that the largest difference is during the passage of the hydrograph peak discharge. In this case, the maximum value of salinity difference in the surface and the bottom layers is calculated as 6 PSU for the reservoir R2 × 4.5. Finally, as the inflow discharge decreases, more intense layer mixing occurs, resulting in a constant value at the column of water.

Temperature affects the density variation of water in the reservoir during desalinization. Figure 13 shows the temperature of the surface water in different coastal reservoirs after passing the hydrograph peak discharge. In the worst case, the maximum temperature difference was less than 1 degree Celsius. However, results show that a shallower reservoir has a more uniform temperature gradient.

Figure 14 shows the temporal changes of temperature at a point in all models. According to this figure, daily changes in water temperature occur. In addition, the highest rate of temperature change in all reservoirs is related to the passage of the hydrograph peak discharge. Finally, results show that the values of temporal temperature changes in depth of the reservoirs are not significant and in the worst case are less than 0.5 degree Celsius.

After determination of water salinity and temperature, one can calculate water density. Figure 15 shows the temporal variations of water density at reservoir depth. From this figure, the water density has a decreasing trend with time mainly due to the reduction in water salinity in the reservoir. Hence, the value of water density decreases from 1,009.35 kg/m³ at the beginning of modelling to its lowest value of 1,001.23 kg/m³ for reservoir R1 × 4.5 and 1,006.52 kg/m³ for reservoir R2 × 9. However, water density increases slightly due to increasing water salinity from reservoir surface to bed. This increase is about 2 kg/m³ in the worst case, which is related to the deeper reservoirs R2 × 4.5 and R2 × 9.

In the present work, the effect of coastal reservoir size on flow and water quality parameters such as water level, flow velocity, salinity, temperature and water density were studied using the 3D numerical software MIKE3. The case study included four types of coastal reservoirs located on the shores of the Caspian Sea and near the mouth of the Tajan River. A flood with 5-years return period in Tajan River was considered as the inflow boundary condition for desalinization of the reservoirs. In each of these models, changes in different parameters of flow, salinity, temperature and water density were investigated. In the following, a summary of the results obtained for each parameter is presented.

• (1)

  The maximum water surface level and flow velocity were about 2 m and 0.1 m/s for the smallest reservoir R1 × 4.5.

• (2)

  Water salinity decreased gradually with progression of fresh water into the reservoir. In addition, the shallower the reservoir the higher the rate of salinity decrease. The highest salinity (about 8 PSU) at the end of the simulation was related to the reservoir R2 × 4.5 and the lowest value (about 2 PSU) was related to the smallest reservoir R1 × 4.5. The salinity stratification of shallow reservoirs was much more intense than deep reservoirs and mainly is a function of the flow discharge.

• (3)

  The distribution of water temperature in the shallow reservoir was more uniform than in the deep ones. However, the temperature difference in different parts of the reservoirs was not so large and in the worst case was less than 1 degree Celsius. Therefore, this parameter could not have a meaningful effect on water density.

• (4)

  The value of water density decreased from 1,009.35 kg/m³ at the beginning of simulation to its lowest value 1,001.23 kg/m³ in the smallest reservoir R1 × 4.5. Due to more salinity, the water density slightly increased in deeper layers.

In the present study, an attempt was made to answer some of the ambiguities about the desalination process of coastal reservoirs. Various issues remain about the performance and operation of these structures, such as seawater intrusion, river pollution depletion management, sedimentation and hydraulic conditions of the river, and access channel for further studies.

This work was supported by Shahid Rajaee Teacher Training University under contract number 13239.

All relevant data are included in the paper or its Supplementary Information.