Modeling soil salinization at the downstream of a lowland reservoir

In arid inland regions, because the annual evapotranspiration potential is higher than the annual precipitation level, artificial irrigation is vital for agricultural production (Scanlon & Goldsmith 1997; Li et al. 2017; Liu et al. 2018). In such an environment, rainfall events are scarce, and available water is insufficient (Garg & Wani 2013; Yang et al. 2018; Yuan et al. 2018); therefore, the conservation and management of surface water resources are necessary. In China's Xinjiang province, lowland reservoirs (as opposed to mountain reservoirs) are built to meet the demands for irrigation, drinking, and industrial water. There are more than 578 reservoirs of various types in Xinjiang (Wang et al. 2014), most of which are lowland reservoirs, with a total capacity of approximately 1.35 × 10⁴ km³.

A recognized problem regarding lowland reservoirs is that they pose threats of increased soil salinization due to their high leakage rates which raise local groundwater (GW) levels (e.g., Cui & Shao 2005; Grunberger et al. 2008; Zarei et al. 2009; Nachshon 2016). The resulting salt crust on the surface can be more than 10 cm, up to tens of cm thick (Weisbrod & Dragila 2006; Hammecker et al. 2012). Thus, lowland reservoirs contributing to salinity is already a serious issue in arid and semi-arid regions (Jorenush & Sepaskhah 2003; Gassama et al. 2012; Hammecker et al. 2012). For example, in northwestern China, the saline area covers more than 2 million ha and accounts for approximately one-third of the total saline area in the country (Cui & Shao 2005; Yuan et al. 2018). Lowland reservoirs are advantageous because of low water costs and convenient transportation, so they cannot be replaced entirely by mountain reservoirs that do not cause soil salinization (Yan et al. 2015). Therefore, the impact of lowland reservoirs must be controlled by proper planning and management. This requires an understanding of the hydrogeological conditions near a lowland reservoir in both the GW and the vadose zones.

Modeling studies have been carried out to assess the relationship between soil salinization and the GW table (Post & Kooi 2003; Zhao et al. 2005; Goncalves et al. 2006; Magri et al. 2009; Hammecker et al. 2012; Mirlas 2012). Goncalves et al. (2006) used the HYDRUS-1D model to analyze water flow and solute transport and evaluate salinization and alkalization hazards. Mirlas (2012) combined MODFLOW with the geodatabase and provided reliable information on soil salinization factors in irrigated fields, in addition to predicting the timing, extent, and duration of shallow water tables. Hammecker et al. (2012) conducted numerical modeling of water flow, based on HYDRUS-3D, to show that managing the depth of flooding within a plot can significantly reduce the outbreak of saline plumes. HYDRUS-1D has been widely used to simulate all soil water hydrological processes (Xu et al. 2005; Zhang et al. 2011; Fan et al. 2012; Tafteh & Sepaskhah 2012), as well as multiple types of solute transport in the vadose zone (Simunek et al. 2003; Simunek & van Genuchten 2008; Simunek et al. 2016). HYDRUS-1D has also proved to be a successful tool for describing water content and solute concentrations when in situ measurements (Ramos et al. 2011). In this study, the HYDRUS-1D model was used to describe the movement of water and salinity in the vadose zone.

Although soil salinization and water movement in the vadose zone have been studied extensively, few studies of these properties in lowland reservoir exist, especially those comparing the impact of GW depth and GW salinity on soil salinization through specific models and data. In this study, a 2-year-long field experiment was conducted to monitor GW levels, soil water contents, and salinity relationships before simulating soil moisture and salinity migration in the area surrounding a lowland reservoir. Different GW levels and salinity scenarios were then used to evaluate the effects of GW depth and salinity on soil salinization.

The study site (44°12′N, 87°48′E) is located in northwestern China in the vicinity of Urumqi in Xinjiang province (Figure 1). The region is characterized by an arid continental monsoon climate with an average annual potential evaporation of 2,013.2 mm and an average annual precipitation of 127.6 mm. The annual mean wind speed is 2.1 m/s, with a maximum of 24.0 m/s. The annual mean temperature is 5.7 °C, with a maximum of 43.5 °C and a minimum of −42.2 °C. The annual mean relative humidity and hours of sunlight are approximately 58.6% and 2,962.8 h, respectively.

The reservoir considered here is a typical artificial diversion lowland reservoir, located in the northern part of the Tien Shan Mountains. The reservoir was built on a silty loam foundation and the vertical profile of soil textures can be divided into silt, silty clay, and silty sand from top to bottom, the thicknesses of which are 2.4–3.0, 3.0–13.0, and 13.0–30.0 m, respectively. The construction of the reservoir started in April 2002 and ended in September 2004, when the reservoir successfully began operation. The main objectives of this reservoir are to supply water to the nearby cities of Fukang and Urumqi and to protect the area's ecology and agriculture. The vegetation in this desert environment is scarce, mainly consisting of Haloxylon ammodendron and Salsola soongorica. A conceptual model of water flow transport in the structured soils of this area is shown in Figure 2.

Fifteen GW wells distributed in three typical sediment sections were designed to monitor the fluctuation of the depth to the GW level and total dissolved solids (TDS). The wells are separated by approximately 1,500 m and arranged downstream of the reservoir (Figure 1). To explore the soil water and salt variation trends, different layered monitoring profiles were distributed near each water quality monitoring well site. In order to avoid duplication, the sampling position was selected by circular wrapping near the monitoring well. The backfill was then sampled and compacted at a thickness of approximately 20 cm. The distances between each of the five monitoring wells in each section were 100, 200, 400, 800, and 1,200 m, respectively (Table 1). All monitoring well depths ranged from 11.8 to 20.5 m. Due to the CJ3-3, CJ3-4, and CJ3-5 monitoring wells being too small to collect GW samples, they were used to investigate the soil water content, salinity, and GW depth. The monitoring section of the reservoir was sampled approximately every 2 months during a 2-year period from August 2012 to November 2014.

Soil water content, salinity, and soil samples were collected at increasing depths from eight profile layers: 0–30, 30–60, 60–100, 100–150, 150–200, 200–250, 250–300, and 300–350 cm. A total of 1,703 soil samples from different times were analyzed by the Second Hydrological Engineering Geological Brigade Laboratory at the Xinjiang Bureau of Geology and Mineral Resources. The soil particle size and dry bulk density were determined in the laboratory.

The HYDRUS-1D model package (version 4.16) was used to simulate the coupled processes of water vapor and solution transport in the downstream zone of the reservoir. In the HYDRUS-1D model, vertical soil water movement was simulated using the Richards equation. The average daily soil water contents measured on 28 August 2012 were used as the initial conditions. The upper boundary condition was set as an atmospheric boundary (precipitation and evaporation), while the variable pressure head was set as the lower boundary. The variable pressure head was specified using GW depth observation values collected approximately every 2 months.

Models of solution transport in the vadose zone were based on analytical or numerical solutions of the Fickian-based convection–dispersion equation. The average daily contents of soluble salt measured on 28 August 2012 were used as initial conditions. Soil salinity at the surface of each monitoring well and the measured TDS concentration of the GW (measured approximately every 2 months at each site) were set as the top and bottom concentration boundary conditions, respectively.

HYDRUS-1D modules were used to apply two different approaches to relate solutes in the liquid and solid phases, including the standard HYDRUS-1D module and the major ion chemistry module (UNSATCHEM) (van Genuchten & Simunek 2004). The standard HYDRUS solute transport module was applied to the solution transport considered in this study.

Cumulative soluble salt concentrations and water movement were simulated by HYDRUS-1D for four different scenarios. In all four scenarios, the boundary conditions (i.e., the top and bottom boundary conditions) were uniform, and the soil hydraulic properties were the same as those at well CJ2-1. The simulation period was from 28 August 2012 to 6 November 2013 (436 days). This study was conducted to understand the mechanism and causes of soil salinization in the areas surrounding lowland reservoirs in order to develop adaptation and mitigation guidelines. HYDRUS-1D simulated the soil water contents and salinity transport under different GW depths and GW TDS concentrations.

GW depth can be the most intuitive response to GW-level changes. The three characteristic sediment sections contained 15 GW monitoring wells (Figure 1). In this study, a total of 195 samples were collected from all wells between 28 August 2012 and 6 November 2013, and a further 15 samples collected in recorded at September 2004 were used as a fixed basement datum (Table 2). GW depths in the observation wells during the month of September in 2004, 2012, 2013, and 2014 were compared (Table 2). Due to leakage, the GW level rose significantly after the reservoir was completed in 2004 in the wells that are closer to the reservoir. The GW depths in CJ1-1, CJ1-2, CJ2-1, CJ2-2, CJ2-3, CJ3-1, and CJ3-2 fluctuated quickly over large ranges, while those in CJ1-3, CJ1-4, CJ1-5, CJ2-4, CJ2-5, CJ3-3, CJ3-4, and CJ3-5 fluctuated slowly, which is consistent with their distance from the dam (Table 2). Observation data from 2004, 2012, 2013, and 2014 show that the GW elevation increased continuously, indicating that GW recharge may not have reached a stable state by the end of the monitoring period.

Fifteen GW monitoring wells were arranged downstream of the reservoir, while CJ3-3, CJ3-4, and CJ3-5 were silted up for the experiment. The mineral contents of water generated from the remaining 12 monitoring wells were investigated. A total of 144 samples were collected from all wells during the 2-year monitoring period. The data showed that, during 2 years of operation, the GW salinity in these wells fluctuated only slightly (Table 3). However, the average annual GW salinity in the 12 monitoring wells ranged from 0.32 to 38.83 g/L (Table 3). The GW average annual mineral content was lowest at CJ1-5 (0.32 g/L) and highest at CJ2-1 (38.82 g/L). The minimum GW mineral content was 0.27 g/L, which was at monitoring well CJ1-2 (26 April 2013). The maximum GW salinity was 48.07 g/L at CJ2-2 (3 November 2014).

Because the GW depth of CJ2-1 changed substantially and exhibited high GW salinity (Table 3), the soil properties of CJ2-1 were selected as a case study for the calibration and validation of the HYDRUS-1D model and other simulated scenarios. Soil cores (0–4 m) were excavated from the CJ2-1 site, and four soil horizons were observed. The vertical soil texture profile (0–4 m) was divided into silt loam and sandy clay (Table 4). The soil is composed of 38% sand, 55.25% silt, and only 6.75% clay. Undisturbed soil samples were collected from different soil layers within each soil profile at the beginning of the experiment to measure their properties. The horizon depths and soil properties are shown in Table 4.

Soil properties of site CJ2-1 (Table 4) were analyzed in the laboratory, and initial estimates of the hydraulic parameters were obtained from the soil texture and the bulk density of each horizon using the neural network analysis module of the HYDRUS-1D software (Schaap et al. 1998). Water and salinity contents measured in 786 samples collected in the period of 28 August 2012 to 29 August 2013 were used to calibrate the parameters of Richards' equation and the Fickian-based convection–dispersion equation by using the inverse option in HYDRUS-1D. Values from a further 917 samples collected from 6 November 2013 to 3 November 2014 were then used to validate the model.

Field measured values were also compared with the results of the HYDRUS-1D simulations using the root mean square error (RMSE). RMSE is a good indicator of the presence and extent of outliers or the variance in the differences between the modeled and observed values; it is given in the units of a particular variable. Thus, the CJ2-1 quantitative RMSE statistics were used to evaluate the numerical accuracy of the model. The RMSE varied from 0.000 to 0.053 cm³/cm³ during the calibration period and from 0.009 to 0.039 cm³/cm³ during the validation period (Figure 3). The RMSE varied from 0.56 to 5.66 g/L (calibration) and 1.86 to 7.52 g/L (validation) (Figure 4).

The modeled parameters of water content variation have a better fit to the observed data than those of soil salinity. These RMSE results showed that the model exhibits a good fit; thus, the calibrated soil hydraulic parameters (Table 5) reflect the soil hydraulic transport properties at the study site. Each layer was described using the nine parameters. Table 5 lists the van Genuchten–Mualem parameters for the CJ2-1 field site. Identical soil hydraulic parameters were considered for all experimental scenarios in the subsequent studies, neglecting the probable effect of the spatial variability of soil hydraulic properties on water flow and solute transport.

In order to investigate the movement of soil water and soil salt in the unsaturated zone downstream of the lowland reservoir, four scenarios with different initial conditions were established using distinctly different GW depths and GW salinity values (Table 6): (a) a low GW TDS value (GWS1; average annual value 0.35 g/L from CJ1-2, Table 3) with an average GW depth from CJ2-1 (GWD1; average annual value 3.35 m, Table 2); (b) a high GW TDS value (GWS2; average annual value of 38.82 g/L, Table 3) with an average GW depth from CJ2-1 (GWD1; 3.35 m, Table 2); (c) a low GW TDS (GWS1; 0.35 g/L, Table 3) with a GW depth from CJ2-1 in September 2004 (GWD2; approximately 10.29 m); and (d) a high GW TDS (GWS2; 38.82 g/L, Table 3) with a GW depth from CJ2-1 in September 2004 (GWD2; 10.29 m, Table 2). Soil hydraulic parameters were under the same conditions and estimated for each scenario using the CJ2-1 calibrated model. The atmospheric and variable pressure head were defined as boundary conditions at the top and bottom of each field plot, respectively. Atmospheric boundary conditions were specified using meteorological data obtained from daily weather observations collected by the weather station at the reservoir management office from 2012 to 2014.

Soil water movements under different GW depths and GW mineralization degrees were identified using HYDRUS-1D (Figures 5 and 6). Scenarios A and B present similar and almost identical trends, as do scenarios C and D, as both pairs had the same GW depths (Figure 5). This may be because mineralization had a negligible effect on soil water movement at the same GW depth. By comparing scenario A with scenario C, and scenario B with D, the transport of soil water in soil layers at depths of 30, 60, 100, and 150 cm all exhibited similar curves. However, the tracks of the other four horizons at depths of 200–350 cm were all different; the deepest layer showed a particularly large difference. These results indicate that surface VSWCs (0–60 cm depth) may be critically affected by atmospheric conditions, while VSWCs in deeper layers correlated with GW depth. Thus, it may be because soil water transport is more dependent on the GW depth.

Figure 6 presents the effects of GW depth and GW salinity on VSWCs at different depths for each depth and scenario plotted in Figure 5 using box plots of the full range of values over 436 days. The characteristic values of the box plots at 30 and 60 cm depth were similar, and at depths from 150 to 350 cm showed that VSWCs differed, showing that the effects of GW depth and GW salinity have a little effect on the VSWCs in the topsoil (0–60 cm). The median values of scenarios C and D were similar, as were the median values of scenarios A and B. The VSWCs of soil profiles with the same GW depths and different GW salinity values had similar trends. The box plot medians of scenarios C and D were lower than those of scenarios A and B. These findings may also verify that GW salinity has a little effect on soil water migration and may indicate that the GW depth has a large impact on soil water transportation (Figure 6).

All these results indicated that GW salinity has a little effect on soil water migration, surface VSWCs (0–60 cm) are critically affected by atmospheric conditions, and GW depth has a substantial impact on soil water transportation. The GW in arid regions exhibits overall high salinities; thus, these findings reveal that GW depth plays a more critical role in soil water movement than GW mineralization. As lowland reservoirs with high leakage rates are now widely used for irrigation in arid areas, rising GW levels are a major issue for salinization in the surrounding region.

Soil salinity transport under different GW depths and mineralization conditions (Figures 7 and 8) reveals a substantial difference in the distributions of soluble salt content and soil water content. Scenarios A and B showed similar soil salt migration processes at the 0–250 cm soil layer, while scenario A (Figure 7) shows that the soil–soluble salt contents in the deepest two soil layers (300 and 350 cm) exhibited almost linearly decreasing trends, while the other scenarios presented slight fluctuations. These results indicate that GW salinity has a substantial influence on the transport of soil salts close to the GW table. Scenarios C and D exhibited approximately the same trend in all eight soil profiles, and higher soil salinity values at 0–60 cm depth may indicate that GW salinity had a little impact on the soil–soluble salt transport for a greater GW depth of approximately 10 m. These results agree with previous findings of the high evaporation rates in arid lands playing an important role on superficial salinity transportation (Grunberger et al. 2008) (Figure 7).

The soil salt contents at GWD1 (3.35 m) were higher than those at GWD2 (10.29 m); the movement of soil–soluble salts with GWD1 at different GW salinity values showed different trends at a depth of 250–350 cm, while soil salinity transport was unclear at GWD2, indicating that GW depth played a positive role in soil salinity migration at lower GW levels. With increasing GW depth, the effect of GW salinity on soil salinity migration decreased (Figure 8).

Figure 8 presents the effect of different GW depths and salinity values on the soil–soluble salt content at different soil depths, with a box plot graph for the full range of values from 28 August 2012 to 6 November 2013 for each depth and scenario. For scenarios A and B, the characteristic values of soil–soluble salt content at 0–200 cm depth were equal, and the median value of soluble salt content at 250–350 cm was lower in scenario A than scenario B, i.e., there was a slight difference in soil salinity at 0–200 cm soil depth that had an obvious effect on soil–soluble salt transport at 250–350 cm under the same GW depth as scenarios A and B. For higher GW salinity values, there may have been more accumulation of soil–soluble salt at soil profiles closer to the GW depth (Figure 8). Scenarios C and D had equal eigenvalues of soil–soluble salt contents throughout the soil profile, indicating that the deeper the GW depth, the smaller the GW salinity's influence on soil salinity transport. Moreover, the median value of soil–soluble salt content at 0–250 cm depth in scenario A was substantially higher than that in scenario C, until the trend inverts at 250–350 cm depth. Figure 6 shows that soil water contents (in the 250–350 cm layer) of GWD1 were higher than those of GWD2. Additionally, the median values of all soil–soluble salt contents in scenario B were higher than those in scenario D. These findings provide evidence that soil water contents and GW depth play important roles in soil–soluble ion movement.

All scenarios indicated that GW salinity has a large influence on the transport of soil salts close to and below the GW table. However, GW depth plays a positive role in soil salinity transport. With increasing GW depth, the effect of GW salinity on soil salinity migration decreases. At higher GW salinity levels, there may be a greater accumulation of soil–soluble salt closest to the GW depth (Figure 8). Additionally, soil water contents play an important role in soil–soluble ion movement at 0–250 cm, which is consistent with soil salinity values at 0–60 cm.

These model simulations are very helpful for evaluating the mechanisms of soil salinization by saline water downstream of a lowland reservoir. Based on model predictions, soil moisture content plays an important role in soil–soluble ion movement, and the GW level also plays a key role in soil salinization (Li et al. 2018). Anti-seepage measures significantly reduce GW level rise; therefore, the improvement of dam foundations can prevent and control soil salinization in areas surrounding reservoirs (Mao et al. 2016). Another way to mitigate the adverse effects of lowland reservoirs on the downstream soil and water environment is by employing measures for seepage cutoff (Yan et al. 2015). In areas surrounding lowland reservoirs, the construction of cutoff ditches connected to the river can reduce the leakage effect and protect the surrounding land from salinization and swamping. In this study, a one-dimensional model that varies only in the vertical dimension was considered; further research should be conducted using a two-dimensional model that takes lateral variation into account.

Lowland reservoirs are uniquely important and widely used in arid areas for irrigation and freshwater supply. However, their high leakage rates are a substantial threat to soil salinization due to an increasing water table. The long-term sustainability of lowland reservoir operation requires a sound understanding of the local soil salinity distribution. To clarify the mechanism of water and salt behavior, the HYDRUS-1D numerical model was employed to successfully simulate four scenarios characterized by different GW salinities and depths.

The HYDRUS-1D simulation results for the four scenarios revealed that GW salinity has a little effect on soil water movement, and that GW depth plays a more critical role in soil water movement. Soil moisture contents and solute transport in the surface soil depend mainly on atmospheric conditions. Soil moisture content plays an important role in soil–soluble ion movement. Rising GW levels present the largest problem for salinization around lowland reservoirs. The understanding of the relationships between GW depth, salinity, soil water content, and soil–soluble ion movement has been enhanced through the simulations evaluated in this study. These findings can be used to develop safety measures to prevent water logging and salinization, as well as for adaptation and mitigation guidelines to improve reservoir design in areas surrounding lowland reservoirs. HYDRUS-1D proved to be an effective and versatile tool for simulating vertical water and solute transport dimensions.

This research was funded by the National Key R&D Program of China (2017YFA0603702) and the National Natural Science Foundation of China (41890822). This research was also supported by the Strategic Priority Research Program of the Chinese Academy of Sciences (grant no. XDA23040304). We gratefully acknowledge the Second Hydrological Engineering Geological Brigade Laboratory at the Xinjiang Bureau of Geology and Mineral Resources for providing working space and housing. The authors also thank Professor Jinlong Zhou and Dr Xianwen Li for their assistance in fieldwork and data analysis work.