Simulation of real-time variations of saline drainage water: comparing system dynamics with DRAINMOD-S

Soil salinity and lack of good-quality irrigation water are significant problems in saline and semi-arid areas, especially where the water table is shallow. In Khuzestan Province, Iran, almost a third of harvestable lands are saline due to upward flux of shallow groundwater (Jorenush & Sepaskhah 2003). One way to cultivate such lands, which typically lack appropriate drainage, is to improve subsurface drainage. However, a poorly designed and/or implemented drainage network can cause salt accumulation in the topsoil, and, accordingly, reduce soil fertility and crop productivity gradually over time. Accurate and optimized determinations of drain depth, spacing and radius, and so on, thus play vital roles in the drainage system's performance and environmental impact. For example, in Khuzestan, increasing drain depth may increase the salt concentration from shallow saline groundwater. In contrast, reducing drain depth may induce water-logging and salt accumulation near the root zone, and reduce crop yield (Nazari et al. 2008).

Drainage is very important to prevent soil water-logging and control salinity in irrigated cultivation, but can have negative environmental impacts. For instance, traverses (pollutants can reach drainage water by transport through the soil profile), and mixing of high-salinity, low-quality drainage water with surface water (e.g., a river or lake) may pollute surface and subsurface water resources severely. Nonetheless, the use of unconventional water resources – such as local low-quality and/or drainage water – seems vital for irrigation. Due to population growth and global demand for further agricultural products (Nozari & Liaghat 2014). For instance, Barnes (2014) concluded that unconventional water sources have become important worldwide and can be used to supplement irrigation water.

Drainage and salinity control projects have been implemented all round the world. The effects of drainage pipes on land reclamation have been studied by investigating the impact of drainage systems on the water table, water quality, and soil salinity (Ghumman et al. 2011, 2012). Hornbuckle et al. (2007), for example, showed that multi-depth drainage systems controlled water table levels and, thus, leaching, more efficiently than single-depth systems. More specifically, in multi-depth drainage systems, shallow drains were found to have significantly lower drainage water salinities than deep drains because of saline groundwater (Hornbuckle et al. 2007). Calculating the length and depth of subsurface drainage pipes accurately would remove salinity effectively and, consequently, increase crop yield. Optimizing subsurface drainage efficiency has a remarkable impact on soil conditions (Shao et al. 2012). Other studies also show that optimal designs leach more nitrogen from the root zone for deeper drain depth (see, e.g., Smedema 2007).

Tiwari & Goel (2015) reviewed salinity management in India and several other countries, and showed that drain depth and spacing depend on the soil's texture and classification. Wesstrom et al. (2014) used controlled drainage systems to improve drainage water quality. Such systems reduce nutrient losses from agricultural land compared to formal drainage strategies under the same conditions. In another study, Mahjoubi et al. (2013) evaluated controlled-drainage systems in sugarcane fields. They showed that such systems improved crop yield, water consumption and drainage, and reduced drainage water salinity and environmental contamination, such as river pollution.

Recently, Nozari et al. (2017) evaluated the effect of drain depth and spacing on the salinity variation of both drainage and groundwater, with an experimental model. They showed that saline groundwater significantly affects the salinity of drainage water. Shakiba et al. (2013) used a lysimeter (length 186 cm, width 20, depth 90) to investigate the impact of water table head on drainage water salinity and the depth at which irrigation and groundwater mix under tile drainage conditions (the equivalent depth). The equivalent depth is the depth below the drain which radial flow will travel to reach the drain. They demonstrated that a large portion of the drainage water came from groundwater, and that groundwater salinity affected drainage water salinity. When the water table rises, the thickness of the zone leached by infiltration flow increases (Shakiba et al. 2013).

Although field and laboratory experiments can be used to study the effect of various factors on drainage system efficiency, optimal depth, and irrigation scheduling, such observations are costly and time-consuming, restricted to a few simple scenarios in irrigation management, and limited to experimental conditions. In this paper, simulation models and their ability to model drainage and groundwater salinities, water table fluctuation, and drainage discharge are discussed.

Unlike experiments, simulation models can be used effectively to study water management systems using complex scenarios and various boundary conditions (Nozari et al. 2014) if they are calibrated suitably (Mostafazadeh-fard et al. 2009). Various models – for example, SaltMod and Nonlinear MINOS – have been applied to assess irrigation and drainage system management (see, e.g., Bonneton et al. 1985; Das 2000; Oosterbaan 2000; Idris et al. 2006). Among them, DRAINMOD (https://www.bae.ncsu.edu/agricultural-water-management/drainmod/), first developed by Skaggs (1980) to simulate the water table and subsurface drainage outlets in wet areas with shallow water tables, has been evaluated against a variety of soil and climate conditions by numerous researchers. For instance, Wahba et al. (2002) assessed DRAINMOD-S in water table simulation in Egypt, and found a good match between simulation and measurement. Wang et al. (2006) also used DRAINMOD to estimate drainage discharge in Southwestern Ohio and concluded that the model results were reasonably good. In another study, Wahba & Christen (2006) used DRAINMOD-S to simulate daily water table, drainage, and salt output in irrigated areas in Southeast Australia, and showed that various drainage plans and management strategies could be designed in semi-arid regions in Australia using DRAINMOD-S.

System dynamics (SD) models, introduced by Forrester (1961) to understand strategies in complex dynamic systems, have also been proposed for water resource management studies. For example, Luo et al. (2009) invoked the SD approach to study infiltration, evapotranspiration, surface runoff, and capillary rise under rice paddies downstream in the Yellow River basin, China. Analysis and feedback relationships between hydrologic processes indicated that the simulations agreed well with their observations. Nozari et al. (2014) applied the SD approach to simulate crop yield and estimated product benefit at the Right Abshar irrigation network, Esfahan, Iran. The SD model showed promising ability to simulate irrigation networks and cropping patterns under different scenarios. Sadeghi Khalegh Abadi et al. (2014) used system dynamics (VENSIM software) to analyze sustainable water resources in the Karkheh basin, developing a comprehensive dynamic model including population and economic activity, as well as water and water resources demand. Pourfallah Koushali et al. (2015) used the SD approach to model water resources in the Sefīd-Rūd watershed, northern Iran, under three different scenarios. Matinzadeh et al. (2017) used SD to develop a simulation model of nitrogen dynamics in a shallow aquifer. They showed that the approach was appropriate for understanding the nonlinear behavior of complex systems over time and could be used to model fertilizer management to prevent environmental pollution.

DRAINMOD-S has been used widely to model drainage and groundwater salinities, as well as water table fluctuations (Skaggs 1980). Studies have shown that the SD approach can be an excellent way to simulate water and soil systems (Saysel & Barlas 2001; Nozari et al. 2014). Although Nozari & Liaghat (2014) simulated drainage water quality via SD, to the best of the authors' knowledge, its accuracy has not been compared to that of DRAINMOD-S.

The main objectives of this study were to: simulate drainage water quality, groundwater salinity, water table fluctuation due to irrigation and precipitation, and drainage discharge using both the SD approach and DRAINMOD-S, and evaluate their accuracies by comparing the simulations with experimental data. The experimental observations are described first in the paper, then the SD and DRAINMOD-S simulations are discussed. It is noted in this context that the experiments have limitations, including high cost, inability to perform in large and complex scenarios, and limited accuracy of results depending on the region and test conditions. Simulation models can be used as a tool to solve these issues, the accuracy of these models' outcomes depending largely on the accuracy of input data, and the quality and validity of the assumptions that are built into them. They can be used to simulate a variety of irrigation management scenarios, if calibrated appropriately.

The experimental observations used were collected in 2008 from the ARC1-18 and ARC2-5 units at Amir Kabir Sugarcane Research Center, Khuzestan Province, Iran. Three rows of piezometers were installed 40, 120, and 200 m from the collector in the ARC1-18 unit, and 100, 250, and 375 m in the ARC2-5 unit, with three in each row in ARC1-18 and seven in ARC2-5. During irrigation, daily water table fluctuation, drainage discharge, irrigation water salinity, and the water salinity in each piezometer were measured, as well as drainage water salinity. Drainage discharge and water table level were recorded simultaneously one or two days after irrigation, via a manhole at the collector entrance. On the basis of the USDA soil classification (Burt 2011), the dominant soil texture class in both units is silty clay with approximately 51% clay, 46% silt, and 4% sand. To determine the representative soil moisture characteristic curve for each unit, disturbed soil samples were collected at different locations, mixed, air-dried, and ground to pass a 2 mm sieve. The average soil moisture characteristic curves for the two units are shown in Figure 1.

Table 1 is a summary of the initial soil profile salinity at different depths and the drainage system's salient characteristics are shown in Table 2. The essential soil information, as input data, is hydraulic conductivity. The auger-hole and inverse auger-hole methods were used to measure hydraulic conductivity below and above the water table, respectively, and the values reported in Table 2 are the averages, collected over several measurements, from various locations in the two units.

In both units, sugarcane rows were spaced 1 m apart. The sugarcane species is perennial, with root length between about 0 and 60 cm. To run simulations using both SD and DRAINMOD-S, the crop coefficient was assumed to vary throughout the growing season – see Table 3.

Irrigation started in April 2008 and lasted for six months using water from the Karoon River. The irrigation water's salinity varied from month to month – see Table 4. Meteorological information – hourly rainfall, and daily potential evapotranspiration and temperature – were collected from a local weather station.

System dynamics is a computer-aided means of understanding nonlinear trends in complex systems over time. It is based on feedback processes that are affected by previous behavior. One of the main processes involved is design of the causal loop diagram, illustrating the cause and effect relationships between system variables. A polarity, either positive (+) or negative (−), is assigned to each causal link, indicating its direction. A negative causal link means that the two nodes change in opposite directions – that is, if the node in which the link starts increases, the other node decreases. Equally, a positive causal link means that the two nodes change in the same direction – if the starting node decreases, the other node also decreases (see, e.g., Nozari et al. 2014).

The two basic concepts within the SD approach are stocks and flows. A stock is an accumulation of material that has built up over time, and shows system status and decisions. The system's activities are, thus, controlled by stocks. Flows comprise material, energy, or information entering or leaving a stock over time.

In the present study, the soil profile is taken as a four-layer stock, each of which has its own characteristics. Water table level and groundwater volume are also considered stock variables with distinct inputs and outputs. In particular, deep percolation is assumed to be the input flow to the water table level, while drainage water, upward flux, and penetration to groundwater are output flows.

VENSIM DSS32 version 4.0a was used to develop and analyze the dynamic feedback model (Ventana Systems 2000). This modeling tool, an object-oriented simulation environment, enables the creation of complex models with less difficulty than conventional programming languages.

The input variables for the top layer of the unsaturated zone were precipitation, irrigation, and upward flux, and its outputs evapotranspiration and deep percolation. The input water to lower layers included deep percolation from the upper layer(s) and upward groundwater flux, while their outputs were deep percolation to groundwater and upward flux to the upper layers. It was also assumed that water did not penetrate the lower layer until the upper layer was at field capacity. As the unsaturated zone has four layers all with the same hydraulic properties, the spatial variation of the hydraulic characteristics is assumed to be negligible.

Accurate determination of salt concentration in the root zone is essential since crop growth and, accordingly, agricultural productivity depend largely on dissolved salt concentrations in soils. In this study, the mass balance equation was used to determine the salt concentration in each soil layer. It is noted that the amount of water moving down from the root zone as deep percolation does not remove salt effectively, because some portions pass through soil macro pores and along preferential pathways without being mixed with soil root zone water. As a result, including leaching efficiency may enable better estimates of root zone and deep percolation salinities.

The main mechanisms generally involved in salt transport are advection, dispersion, and molecular diffusion. Although the advection-dispersion equation, also known as ADE (see, e.g., Bear 1972), is only valid in homogeneous media, it was used to model solute transport for the sake of simplicity and comparison with DRAINMOD-S in this study. The dispersivity value – 0.08 m – in ADE was determined by calibrating the SD model using the first 30% of the measured drainage water salinity data.

DRAINMOD has been used widely to simulate water management systems' performance – for example, in controlled and uncontrolled drainage, and underground irrigation, as well as combinations of them. The model is based on the water balance equation, which includes two parts relating to surface and subsurface flow. DRAINMOD-S is a sub-model, and can simulate soil and drainage water salinities simultaneously, and the governing mass balance equation to simulate the soil profile salinity is ADE.

As for SD, the inputs for DRAINMOD-S are climate data, soil physical properties, initial soil salinity, drainage system parameters, and crop information (Kandil et al. 1995). Kandil (1992) calibrated DRAINMOD-S by changing the dispersivity and comparing the results with collected soil salinity profiles, and found that 0.5 m dispersivity gave the best match between simulated and real measurements. In this study, dispersivity values between 0.03 and 0.5 m were used, and the modeled results compared to measured drainage water salinity data to calibrate DRAINMOD-S. In fact, both Kiran (2009) and Kale (2011) report that DRAINMOD-S results were not sensitive to the dispersivity value, while Brevé et al. (1997) and Wang et al. (2005) showed that DRAINMOD-N is also insensitive to the dispersivity value. Accordingly dispersivity was set at 0.08 m in this study.

The simulated water table fluctuations using DRAINMOD-S and SD, and the real ones, are shown in Figure 2. Reasonably good agreement can be seen between the real levels and the DRAINMOD-S simulations, although SD estimated the water table substantially more accurately in both units. In general, the SD simulations match the field measurements better than those of DRAINMOD-S. For ARC1-18 (Figure 2(a)), RMSE and SE were determined, respectively, as 15.3 cm and 0.08 for SD, and 21.2 cm and 0.12 for DRAINMOD-S. ARC2-5 was slightly different, with RMSE = 15.2 cm and SE = 0.09 for SD, and RMSE = 28.5 cm and SE = 0.17 for DRAINMOD-S (Figure 2(b)). The RMSE value for DRAINMOD-S found in this study is within the 6.5 to 36 cm range reported by Maurizio et al. (2000), who simulated the water table level in the Veneto region of Italy using DRAINMOD. As can be seen, the statistical parameters calculated for SD are considerably lower than those determined for DRAINMOD-S and demonstrate that SD gave better water table estimates than DRAINMOD-S.

The drainage discharge simulations using DRAINMOD-S and SD are shown in Figure 3, along with the measured levels. Again, the simulated discharges from SD are more accurate than those of DRAINMOD-S. The RMSE and SE values are respectively 1.1 L/s and 0.57 for SD, and 1.2 L/s and 0.67 for DRAINMOD-S in ARC1-18 (Figure 3(a)). For ARC2-5, RMSE = 3.1 L/s and SE = 0.61 for SD, and RMSE = 3.3 L/s and SE = 0.65 for DRAINMOD-S (Figure 3(b)) - that is, SD estimated drainage discharge slightly better than DRAINMOD-S.

Singh et al. (2006) simulated subsurface drainage systems using DRAINMOD for two experimental plots near Gilmore City, USA. They reported RMSE = 2.24, 0.98, 1.24, and 1.51 cm/d for the drainage coefficient for the four years from 1990 to 1993. Given that their average drainage coefficients were 7.25, 3.55, 3.37, and 4.40 cm/day, respectively, the approximate dimensionless SE values would be 0.31, 0.28, 0.37, and 0.34 for 1990 to 1993, around half of those found in this study.

Figure 4 shows the measured and simulated drainage water salinities reflected in EC values. The RMSE and SE values are respectively 3.36 dS/m and 0.24 for SD, and 4.76 dS/m and 0.34 for DRAINMOD-S in ARC1-18 (Figure 4(a)). In ARC2-5, RMSE = 0.73 dS/m and SE = 0.20 for SD, and RMSE = 1.57 dS/m and SE = 0.44 for DRAINMOD-S (Figure 4(b)). As shown in Figure 4, the SD drainage water salinity simulation has fairly gradual descending and ascending trends in ARC1-18 and ARC2-5, respectively, during irrigation, while the real salinity measurements fluctuate considerably more and more frequently over time. The measured drainage water salinity fluctuations arise because of the impact of variations in water table level on the mixing depth where radial flow occurs under tile drainage conditions. For example, Shakiba et al. (2013) demonstrated that, as the water table rises after irrigation, the mixing depth also rises. Because of the salt concentration below the drain, the greater the mixing depth, the higher the drainage water salinity. The effluent quality index (Rathnayake & Tanyimboh 2015) can be determined to quantify water quality but, for consistency with the DRAINMOD-S results, salinity was used in this study.

Figure 4 also shows that the general trend of drainage water salinity variations as simulated by SD is less than that simulated by DRAINMOD-S, particularly for ARC1-18. Over short time scales (e.g., Figure 4(a) around the 115 day mark), the DRAINMOD-S drainage water salinity fluctuates substantially from zero to almost 20 dS/m. Later in the simulation period; however, the DRAINMOD-S salinity simulation is more restricted, in the range between 7 and 12 dS/m (Figure 4(a)). In ARC2-5 (Figure 4(b)), the DRAINMOD-S salinity simulation, like the measured one, fluctuates in a more restricted range between 2 and 5 dS/m over almost the entire period.

The actual – that is, measured – drainage water salinity ranges from about 7 to 20 dS/m and 2 to 5 dS/m in ARC1-18 and ARC2-5, respectively. This difference arises because ARC1-18 has been operating for a shorter time than ARC2-5. The groundwater salinity is 20 dS/m in ARC1-18 and 10.6 dS/m in ARC2-5 – see Table 1, which confirms that ARC2-5 has reached equilibrium due to its long record of irrigation and leaching, for example, compared to ARC1-18.

The simulated and actual groundwater salinities are compared in Figure 5 with respect to ARC1-18, at depths of 220, 260, and 300 cm, and the values of RMSE and SE for the two models are reported in Table 5. The average values of RMSE and SE over all three depths are respectively 0.43 dS/m and 0.11 for SD, and 0.82 dS/m and 0.20 for DRAINMOD-S. The latter RMSE and SE values are nearly twice the former. As both models are based on the water balance approach and use the same input parameters, this considerable discrepancy is probably due to the insensitivity of DRAINMOD-S to dispersivity and the calibration process.

As shown in Figure 5, the SD's simulated groundwater salinity agrees well with the measured values over the entire period. In contrast, DRAINMOD-S systematically overestimates groundwater salinity over short time-scales and underestimates it substantially over long ones, particularly at shallower soil depths. At short time-scales, the measured groundwater salinity decreases, initially, to something close to the salinity of the irrigation water, that is, 2 dS/m, particularly at shallow depths (220 and 260 cm), but increases to higher values at longer time-scales. The salinity variations and equilibrium time differ from one layer to another; however, the deeper layers are more saline, and the equilibrium time longer. This is because the salinity in layers near the drain is washed away by irrigation water and radial flow more quickly than elsewhere, so that layers close to the drain reach equilibrium more quickly.

In this study, the system dynamic (SD) approach was used to simulate water table fluctuations, drainage discharge, and drainage water and groundwater salinity, and to compare with experimental measurements collected from ARC1-18 and ARC2-5, field units at Amir Kabir Sugarcane Research Center, Khuzestan Province, Iran. Comparisons were also made with output from DRAINMOD-S. For the model evaluation the root mean square error (RMSE) and standard error (SE) were calculated, with the following results:

- Water table

ARC1-18 (and ARC2-5), RMSE = 15.3 cm (15.2 cm), and SE = 0.08 (0.09) for SD, and RMSE = 21.2 cm (28.5 cm) and SE = 0.12 (0.17) for DRAINMOD-S.

- Drainage discharge

ARC1-18 (and ARC2-5), RMSE = 1.1 L/s (3.1 L/s), and SE = 0.57 (0.61) for SD, and RMSE = 1.2 L/s (3.3 L/s) and SE = 0.67 (0.65) for DRAINMOD-S.

- Drainage water salinity

ARC1-18 (and ARC2-5), RMSE = 3.36 dS/m (0.73 dS/m), and SE = 0.24 (0.20) for SD, and RMSE = 4.76 dS/m (1.57 dS/m) and SE = 0.34 (0.44) for DRAINMOD-S.

- Groundwater salinity

ARC1-18, the average RMSE = 0.43 dS/m, and (average) SE = 0.11 for SD, and average RMSE = 0.82 dS/m and (average) SE = 0.20 for DRAINMOD-S.

The results showed that the SD approach simulated the water table level, drainage salinity, and groundwater salinity considerably more accurately than DRAINMOD-S. Although SD estimated drainage discharge slightly better than DRAINMOD-S, the models' accuracy is not significantly different. DRAINMOD-S also simulated salinity typically at values below those measured, while the SD salinity estimates were closer to reality.

None.

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