Calibration and Global Sensitivity Analysis for a Salinity Model Used in Evaluating Fields Irrigated with Treated Wastewater in the Salinas Valley

: Treated wastewater irrigation began two decades ago in the Salinas Valley of California and provides a unique opportunity to evaluate the long-term effects of this strategy on soil salinization. We used data from a long-term ﬁeld experiment that included application of a range of blended water salinity on vegetables, strawberries and artichoke crops using surface and pressurized irrigation systems to calibrate and validate a root zone salinity model. We ﬁrst applied the method of Morris to screen model parameters that have negligible inﬂuence on the output (soil-water electrical conductivity (EC sw )), and then the variance-based method of Sobol to select parameter values and complete model calibration and validation. While model simulations successfully captured long-term trends in soil salinity, model predictions underestimated EC sw for high EC sw samples. The model prediction error for the validation case ranged from 2.6% to 39%. The degree of soil salinization due to continuous application of water with electrical conductivity (EC w ) of 0.57 dS/m to 1.76 dS/m depends on multiple factors; EC w and actual crop evapotranspiration had a positive effect on EC sw , while rainfall amounts and fallow had a negative effect. A 50-year simulation indicated that soil water equilibrium (EC sw ≤ 2dS/m, the initial EC sw ) was reached after 8 to 14 years for vegetable crops irrigated with EC w of 0.95 to 1.76. Annual salt output loads for the 50-year simulation with runoff was a magnitude greater (from 305 to 1028 kg/ha/year) than that in deep percolation (up to 64 kg/ha/year). However, for all sites throughout the 50-year simulation, seasonal root zone salinity (saturated paste extract) did not exceed thresholds for salt tolerance for the selected crop rotations for the range of blended applied water salinities.


Introduction
Salinization of soils, groundwater, and surface waters from irrigation is a well-known problem often associated with the decline of ancient civilizations dependent on irrigated agriculture around the world, such as Mesopotamia [1]. Today, the salinity problem associated with irrigation in low rainfall regions continues to have numerous grave economic, social, and political consequences. For example, there is a high economic cost associated with salinity; the US Bureau of Reclamation spends $32 million annually to limit salt additions to the Colorado River and the Natural Resource Natural Resource Conservation Service-US Department of Agriculture (USDA-NRCS) spends some $13 million annually to control salinity in irrigation programs across the upper Colorado River Basin [2]. Simultaneously, as competition for available freshwater resources intensifies and use of treated wastewater (recycled • Complete a model calibration and validation using parameters to which model outputs are most sensitive as a guide: and • Predict long-term (five decades) root zone salinity, salt output load with deep percolation and salt output load with surface runoff from fields using treated wastewater for irrigation.

Study Area
The study area in Castroville, overlies two main aquifers referred to as the "180-foot" and "400-foot" aquifers, respectively. These were formed from fluvial sands and gravels associated with the old Salinas River channel and possible delta conditions. Above and between the two aquifers are deposits of blue clay overlying the "180-foot" aquifer range from 8-m thick at Salinas to more than 30-m thick at Castroville [26][27][28]. The typical overlying soil profile in the study area is comprised of Pacheco sandy to clay loams as summarized in Table 1. Table 1. Average soil profile texture variations in the study region [29]. We assembled the base data for the model using the estimates for saturated hydraulic conductivity, organic matter content, soil texture, and bulk density from SoilWeb [29] an interactive webtool to access detailed soil survey data (SSURGO). We then determined saturation water content (Ts), wilting point (Twp), and field capacity (Tfc) according to soil texture using artificial neural networks techniques implemented in Hydrus-1D [30]. Meteorological records for 1983 to 2014 were taken from the local (California Irrigation Management Information System (CIMIS) station number 19 [31] and the average monthly rainfall and grass-reference ET o are shown in Figure 1. Average monthly reference ET o exceeded average monthly rainfall during April through November, annual ET o ranged from 862.3 mm to 1072.6 mm (±5.3%). Rainfall was concentrated from November to March (87% of annual rainfall) with annual rainfall depths ranging from 134.5 mm to 1026 mm (±47.1%) as shown in Figure 2.
The main crops grown in the project area are cool-season vegetables (lettuce, broccoli, cauliflower, artichoke, cabbage, spinach, celery) and strawberries. In 1995, the Monterey County Water Resources Agency (MCWRA) passed an ordinance prohibiting extraction of groundwater between sea level and −76.2 m in Salinas and Castroville. In 1998, Monterey County Water Recycling Projects (MCWRP) began delivering recycled water to 486 hectares (12,000 acres) in the northern Salinas Valley. Crop rotations and management practices at the eight sites of the study area are listed in Table A1; the crops include lettuce, broccoli, cauliflower, cabbage, celery, spinach, artichokes, and strawberries. Drip irrigation was used at the control site and vegetable crops were established with sprinklers for 20 to 30 days. At site 2, sprinkler irrigation switched to drip after plant establishment in 2002 while site 3 used sprinklers for vegetables and drip for strawberry. Sites 4 and 5 used sprinklers and drip and sites 6 and 7 used sprinklers for plant establishment and followed by furrow irrigation.
irrigation method and farming practices. At each site, soil samples were collected from depths of 0.03 to 0.30 meters, 0.30 to 0.61 meters and 0.61 to 0.91 meters at four different locations within 1 meter of a designated global positioning system (GPS) point. Sample analysis was done by an independent accredited lab (Valley Tech, Tulare, CA, USA) and included pH, soil water electrical conductivity (ECsw), extractable cations B (boron), Ca, Mg, Na and K, and extractable anions Cl, NO3 (nitrate) and SO4 (sulphate).  Irrigation water salinity varied between sites and years as recycled water was blended with groundwater (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009) or diverted Salinas River water (2010-2012) and the average applied water EC (ECw) at the different sites for these time periods are summarized in Table 2. Tertiary treated wastewater in the study had on average Sodium Adsorption Ratio (SAR) value of 5.58, containing 192.1 mg/L of Na, 246.1 mg/L of Cl, and EC of 1.4 dS/m. The rain salinity (ECp) values were taken from the National Atmospheric Deposition Program station in the Pinnacles National Park located ~322 kilometers east of the study site. The ECp varies by month ranging from 0.001 dS/m to 0.004 dS/m with May having the highest ion deposits with rain. Tertiary treated wastewater effluent from Monterey Regional Water Pollution Control Agency, (MRWPCA) was sampled on a weekly basis to determine the levels of salt present in it before blending with the supplemental well water used to meet peak irrigation demand. Monthly delivery system sampling confirmed the quality of the water received by growers after supplemental well water was Agriculture 2019, 9, 31 6 of 33 added to the recycled water. In addition, the quality of the well water delivered to the control site was sampled monthly. The water samples were analyzed for pH, EC w , Na, Mg, Cl and K (potassium) by an accredited laboratory run by MRWPCA. The one control and seven test sites were randomly distributed throughout Castroville, USA area and were chosen based on soil characteristics, drainage systems, types of crops grown (lettuce, cole crops and strawberries), irrigation method and farming practices. At each site, soil samples were collected from depths of 0.03 to 0.30 m, 0.30 to 0.61 m and 0.61 to 0.91 m at four different locations within 1 m of a designated global positioning system (GPS) point. Sample analysis was done by an independent accredited lab (Valley Tech, Tulare, CA, USA) and included pH, soil water electrical conductivity (EC sw ), extractable cations B (boron), Ca, Mg, Na and K, and extractable anions Cl, NO 3 (nitrate) and SO 4 (sulphate).
Irrigation water salinity varied between sites and years as recycled water was blended with groundwater (2000-2009) or diverted Salinas River water (2010-2012) and the average applied water EC (EC w ) at the different sites for these time periods are summarized in Table 2. Tertiary treated wastewater in the study had on average Sodium Adsorption Ratio (SAR) value of 5.58, containing 192.1 mg/L of Na, 246.1 mg/L of Cl, and EC of 1.4 dS/m. The rain salinity (EC p ) values were taken from the National Atmospheric Deposition Program station in the Pinnacles National Park located 322 km east of the study site. The EC p varies by month ranging from 0.001 dS/m to 0.004 dS/m with May having the highest ion deposits with rain.

Root Zone Salinity Model
We coupled crop growth and soil water models applied across the root-and vadose-zones to simulate both upward flow from shallow water tables as well as downward percolation to the groundwater (see Appendix B for a detailed description of the model). These were combined with a root zone salinity model and used to predict root zone soil salinity (see model configuration in Figure 3). The two driving criteria for model selection included the simplification required to describe the processes mathematically without losing the detail needed to develop realistic results, and the model reliance on readily available input data. The Isidoro and Grattan [17] daily time-step model uses a closed-form solution of Reference [32] to describe vertical unsaturated water movement in the root zone and unsaturated zone (Equations (A1) to (A4)). A number of closed-form formulas have been proposed to empirically describe the dependence of unsaturated hydraulic conductivity and water content on pressure head [33][34][35][36][37]. We used the Clapp and Hornberger equation [32] to extend vertical flow through a continuous soil profile to compute the movement of water and salt across the entire vadose zone to the account for shallow groundwater flow processes. Thus, two additional layers from the root zone to the groundwater table were added for both unsaturated and saturated zones.
The crop component of the model includes crop development stages (Table A2), root growth (Equation (A8)), root water uptake and water stress response functions (Equations (A5)-(A7)). Evapotranspiration (ET) includes a combination of two separate processes whereby water is lost from the soil surface by evaporation and from the root zone by crop transpiration. The crop ET is calculated as the product of the reference ET o and the estimated crop coefficient (K c ) that depends on crop characteristics, vegetative growth state, canopy cover, and height as well as surface-soil properties (Equations (A5)-(A7). In each layer (k), the actual crop ET can be lower than potential ETc(k) due to water stress, which depends on the soil water content and the sensitivity of the crop to low water contents, accounted for through the crop-specific parameter p: the ratio of readily available soil water (RAW) to total available water (TAW) (p = RAW/TAW) [38]. Plant water uptake was assumed to be a descending extraction pattern that depends on irrigation frequency such that greater uptake is at the top quarter of the root zone [18,[41][42][43]. Plant growth and root development parameters are summarized in Table A2 and Equation (B8) [38,[44][45][46][47] and we assumed that strawberries were planted on 1.3 m wide raised beds as is common in the region. The model was calibrated for soil salinity generation due to dissolution and a dissolution rate is used to account for these processes. Irrigation and rainwater salinity are specified by the user and the model neglects plant root uptake of salts. While preferential flow and irrigation non-uniformity may also be important features, they were beyond the scope of this model.
Each simulation extended for a period of 13 and 50 water years (1 October to 30 September) and more importantly the model simulates carry over effect from one year to the next. The surface boundary conditions of rainfall, irrigation and ETo were specified daily together with irrigation and rain water EC. The lower boundary conditions were specified as fixed water table depth and  groundwater salinity ECgw; though fixed water table depths are unlikely in the field, water table  fluctuations are assumed to be dampened by capillary rise and evaporation from the water table. The model is written in R to make it more widely accessible to water managers and possibly growers.

Calibration and Sensitivity Analyses
Parameter sensitivity analyses provide insight into those model parameters that are most critical towards approximating measured results and are often used to help focus field sampling or measurement efforts and/or refinement of modeled processes. Here, we take a different approach and first used a global sensitivity analysis that considers variations within the entire variability space of the input factors. We used the elementary effect (EE) method for screening important input factors among the 33 factors initially considered important. Finally, the variance-based "Sobol" method was used with those factors determined to be significant from the Morris screening method for factor fixing and to identify those factors which, left free to vary over their range of uncertainty, make no significant contribution to the variance of the output results of interest. We applied the modified Isidoro and Grattan model to the 13 years of soil salinity observations up to 91.4 cm below the ground The model domain consists of a one-dimensional vertical 7.7-m soil profile, representing the crop root zone and the unsaturated zone that overly a fixed saturated (water table) zone. The domain is discretized so that the clay loam root zone is divided into four quarters of equal depth to enable determination of plant water uptake fractions. The unsaturated zone below the root zone is divided into two layers: a sandy-loam layer immediately below the root zone and a silty clay-loam between the sandy-loam and the capillary fringe. Both upward flow from and downward flow to shallow groundwater is possible in the model. Surface runoff depths were calculated using the Soil Conservation Service (SCS) runoff method [39,40] see Equations (A9) and (A10). Equations (A11) and (A12) detail the soil-water balance simulations.
Salt balance calculations were performed in conjunction with the soil-water balance assuming complete mixing of water entering each layer with that already stored in that layer (Equations (A13)-(A19)). The soil-water EC is used as a salinity indicator, implicitly assuming there is a unique relationship between EC and total dissolved solids (TDS), and that EC behaves as a non-reactive (conservative) solute. Salinity of irrigation water (EC w ) and precipitation (EC p ) are input values. The mass of salts in layer k (Z(k)) is estimated from the product EC sw (k) W(k), where EC sw is the electrical conductivity of the soil water in that layer.
Plant water uptake was assumed to be a descending extraction pattern that depends on irrigation frequency such that greater uptake is at the top quarter of the root zone [18,[41][42][43]. Plant growth and root development parameters are summarized in Table A2 and Equation (A8) [38,[44][45][46][47] and we assumed that strawberries were planted on 1.3 m wide raised beds as is common in the region. The model was calibrated for soil salinity generation due to dissolution and a dissolution rate is used to account for these processes. Irrigation and rainwater salinity are specified by the user and the model neglects plant root uptake of salts. While preferential flow and irrigation non-uniformity may also be important features, they were beyond the scope of this model.
Each simulation extended for a period of 13 and 50 water years (1 October to 30 September) and more importantly the model simulates carry over effect from one year to the next. The surface boundary conditions of rainfall, irrigation and ET o were specified daily together with irrigation and rain water EC. The lower boundary conditions were specified as fixed water table depth and groundwater salinity ECgw; though fixed water table depths are unlikely in the field, water table fluctuations are assumed to be dampened by capillary rise and evaporation from the water table. The model is written in R to make it more widely accessible to water managers and possibly growers.

Calibration and Sensitivity Analyses
Parameter sensitivity analyses provide insight into those model parameters that are most critical towards approximating measured results and are often used to help focus field sampling or measurement efforts and/or refinement of modeled processes. Here, we take a different approach and first used a global sensitivity analysis that considers variations within the entire variability space of the input factors. We used the elementary effect (EE) method for screening important input factors among the 33 factors initially considered important. Finally, the variance-based "Sobol" method was used with those factors determined to be significant from the Morris screening method for factor fixing and to identify those factors which, left free to vary over their range of uncertainty, make no significant contribution to the variance of the output results of interest. We applied the modified Isidoro and Grattan model to the 13 years of soil salinity observations up to 91.4 cm below the ground surface. Measured data from the control site, irrigated only with available groundwater, was used for calibration and one of the other eight test sites using predominantly recycled water (94-98% recycled/freshwater blend) was used for model validation so as to bracket possible model predictions. The calibrated model was then used to assess the long-term (50 years) salinity outcomes of variable management strategies including cropping patterns, irrigation technology, and irrigation water quality.
We use the Morris and Sobol's methods to support model calibration as shown in Figure 4. The sensitivity analysis was used to address the following questions:

•
What input parameters or factors cause the largest variation in the output? • Are there any factors whose variability has a negligible effect on the output? • Are there interactions that amplify or dampen the variability induced by individual factors?
Agriculture 2019, 9, 31 9 of 34 The calibrated model was then used to assess the long-term (50 years) salinity outcomes of variable management strategies including cropping patterns, irrigation technology, and irrigation water quality. We use the Morris and Sobol's methods to support model calibration as shown in Figure 4. The sensitivity analysis was used to address the following questions:


What input parameters or factors cause the largest variation in the output?  Are there any factors whose variability has a negligible effect on the output?  Are there interactions that amplify or dampen the variability induced by individual factors?
Sensitivity Analysis Library (SALib version 1.1. 3) an open source library written in Python, was used for performing the sensitivity analyses [48]. The model input variables considered for the sensitivity analysis are listed in Table 3 and parameters showing the greatest sensitivity were selected for calibration.   3) an open source library written in Python, was used for performing the sensitivity analyses [48]. The model input variables considered for the sensitivity analysis are listed in Table 3 and parameters showing the greatest sensitivity were selected for calibration. Different crops have different water uptake patterns, but all take water from wherever it is most readily available within the rooting depth. The root zone water-uptake pattern depends on irrigation frequency. With infrequent irrigations, the typical extraction fractions by root zone layer is 40-30-20-10%. For frequent drip or sprinkler irrigation, the water uptake fractions are skewed towards greater uptake from the upper root zone, or a 60-30-7-3% uptake pattern [18]; this pattern is assumed in many classical analysis of saline soils [42]. Some have suggested use of an exponential model that specifies a greater proportion of uptake near the soil surface, that is, uptake fractions of 71-20-6-3% [41,43]. Ranges for the crop related data, including crop coefficients (Kc), rooting depths (Zr), and average fractions of total available water that can be depleted from the root zone before moisture stress (p), were taken from the Food and Agriculture Organization of the United Nations (FAO)-56 [38]. Estimates of strawberry crop coefficients were found in Reference [47].
All field sites considered were situated on Pacheco clay, clay-loam, and sandy loam soils with ranges of soil texture, available water content, bulk density, organic matter content, and saturated hydraulic conductivity (Ks) taken from SSURGO soil surveys as summarized previously in Table 1. Soil hydraulic properties required for model application were inferred from the soil survey information. We fitted that information to the van Genuchten model using a Multiobjective Retention Curve Estimator (MORE) based on the Multiobjective Shuffle Complex Evolution Metropolis (MOSCEM-AU) algorithm implemented in HYDRUS-1D [30]. The fraction "α" of the excess water that drained the first day (0 < α < 1) was calculated from the soil texture in the layer through an empirical relation obtained to match the results presented in [49]. Grismer [20] provides a relationship between capillary fringe heights (Hd) and saturated intrinsic permeability for different soil textures. Groundwater levels were taken from regular measurements by the MCWRA at well 13S/02E-32E05 about 5 km west of study area. Estimated range of groundwater salinity was taken from Reference [50].
The primary output variable of concern was the soil-water EC sw as determined for root zone layers up to 0.91 m deep. As the Isidoro and Grattan model is a dynamic model, the "output" term in the sensitivity analysis does not refer to the range of spatial and temporal distribution of EC sw but to a summary variable. In this case, the root mean square error (RMSE) that is obtained as a scalar function of the simulated time series output EC sw values. As such, for calibration, the objective function minimized the RMSE associated with model prediction.

Model Parameter Screening-Elementary Effects Method
The elementary effects (EE) method is an effective way of screening for important input factors contained in a model [51]. The fundamental concept of this method involves deriving measures of global sensitivity from a set of local derivatives, or elementary effects, sampled on a grid throughout the parameter space [52]. It is based on one-at-a-time (OAT) analysis, in which each parameter Xi is perturbed along a grid of size ∆i to create a trajectory through the parameter space. For a given value of X, the elementary effect of the ith input factor is defined as: where Y(X 1 , X 2 , . . . , X k ) is a prior point in the trajectory and X = X 1 , X 2 , . . . , X k is any selected value in the parameter space such that the transformed point is still in the parameter range for each index i = 1, . . . , k. The sensitivity measures µ and σ are the mean and the standard deviation of the distribution of EEi proposed by Morris. Mean parameter (µ) assesses the overall influence of the factor on the output parameter of interest; σ assesses the extent to which parameters interact. Thus, a small σ value implies that the effect of Xi is almost independent of the values taken by other factors; on the other hand, a large σ indicates that a factor is interacting with others because its sensitivity changes across the variability space. Campolongo and others [53] suggest that µ* is a good proxy of the total sensitivity index, a measure of the overall effect of a factor on the output parameter inclusive of interactions. We analyzed µ* for all input factors to screen out non-influential factors, and then performed a variance-based analysis with the remaining important factors. Once trajectories are sampled, the resulting r elementary effects per input are available, the statistics µ, σ2 and µ* for each factor are computed as:

Factor Fixing-Sobol's Variance Method
Sobol's sensitivity analysis is a global-variance based method. Sensitivity measures are based on the decomposition of the model output variance to individual parameters and the interaction between parameters [54,55]. Variance-based sensitivity analysis relies on three principles: • input factors are assumed to be stochastic variables of the model that induce a distribution in the output space; • the variance of the output distribution is a good proxy of its uncertainty; and • the contribution to the output variance from a given input factor is a measure of sensitivity.
Contribution to total output variance by individual input factors and their interaction can be written using an ANOVA high-dimensional model representation (HDMR) decomposition [51]: where V(Y) is the total or unconditional variance of the output, the conditional variance; V i is the conditional variance or first-order effect of X i on Y; V ij is the joint effect of X i and X j minus the first order-effects for the same factors. Several variance-based indices can be defined; the first order index represents the main contribution effect of each input factor to the output variance can be determined from: The total order index, S T , a measure of the overall contribution to output variance from an input factor considering its direct effect and its interactions with all other factors and is determined from: where V ∼i is the conditional variance with respect to all the factors but one, i.e., X ∼i . The condition S Ti = 0 is necessary and sufficient for X i to be a noninfluential factor on the output. That is, if S Ti ∼ = 0, then X i can be fixed at any value within it range of uncertainty without appreciably affecting the value of the output variance V(Y). Here, we calculated all of the indices to determine the factors that can be fixed in the calibration process. A recommended sampling technique uses sequences of quasi-random numbers generating n, 2k matrix of random numbers where n is called a base sample and k is the number of input factors. This scheme allows for n (k + 2) model evaluations. We evaluated up to n = 12 and found no changes occurred after n = 10 and concluded that it was sufficient.

Long-Term Salinity Indicators
Salt output loads with deep percolation were used to assess the potential for groundwater resources deterioration and salt output loads with runoff indicate the salinity threats posed by treated wastewater irrigation on the salinization of the Salinas River. Similarly, crops respond to the salinity in the root zone over the entire growing season [18]. Thus, we used seasonal-averaged root zone EC of the saturated extract or EC eS (dS/m), deep percolation S d (kg/ha), and surface runoff S r (kg/ha) salt output loads as key output state variables to describe long-term (decadal) impacts of irrigation with treated wastewater and farm management practices (i.e., applied water quality and depths, irrigation technology, and crop rotations). Annual rainfall mitigates impacts of irrigating with saline water as such accounting for rainfall leaching is important for evaluating long-term dynamics.
Daily EC of the saturated extract (EC e ) in the root zone layers (k) is calculated as: where SP(k) is the saturation percentage (water content of the saturated soil paste expressed on a dry weight basis) for layer k. Traditionally, for most mineral soils it is assumed that field capacity is half of SP, so EC et is the mean EC e of the 4 rootzone layers; daily values EC et are then averaged over the entire growing season yielding the seasonal-averaged root zone EC e (EC eS ): Days in the growing season (9) We applied the model to both meteorological series and management practices for 13 years of cropping practices in the study area. The objective was to simulate a 13-year continuous cropping and provide a multiple-year record of seasonal EC eS and related parameters. Following [15], we assumed a factor of 640 to convert EC (dS/m) into TDS (mg/L) for EC ≤ 5 dS/m and a factor of 800 for EC > 5 dS/m.
The model estimates daily lower boundary water flux (D) along with an estimate of soil-water EC in the bottom layer to determine the daily salt output load associated with deep percolation water calculated as: where S d is the daily drainage salt output load in kg/ha; EC sw (4) is the EC of soil water in the bottom root zone layer in dS/m and 6.4 is the conversion factor assuming a factor of 640 to convert EC (dS/m into TDS (mg/L) and flux in mm. Additionally, the model estimates daily runoff volumes (SR) from the soil surface along with the runoff water EC such that the daily salt output load associated with runoff is calculated as: where S r is the daily runoff salt output load in kg/ha; EC sr is the EC of surface runoff in dS/m and 6.4 is the conversion factor.

Calibration and Validation
The primary objective of model calibration was to capture the long-term soil salinity dynamics in the fields irrigated with treated wastewater in the study region by satisfactorily reproducing the 2000 to 2012 soil-water salinity (EC sw ) data set described by [11]. The study consisted of six test sites and one control site randomly distributed across the Castroville region that were chosen to provide a typical range of soil characteristics, drainage systems, types of crops grown, irrigation methods, and farming practices found in the region. Average annual water quality delivered to each site was determined as well as soil samples collected from depths of 0.03-0. Soil samples were collected following winter rains before the spring planting, during and at the end of the summer growing season prior to winter rains. Saturated paste extracts from these soil samples were analyzed for EC and solute concentrations. The control site received only well (2000-2009) or surface (2010-2012) water and this site was used for calibration, while site 3 received 94-98% recycled water and this site was used for validation (see Table 2). Annual crop rotations that included lettuce, broccoli, cauliflower, cabbage, and strawberry are shown in Figure A1. Control site vegetables were established with sprinklers for 20-30 days and drip irrigated, while at site 3, vegetables were sprinkler irrigated and then drip irrigation was used for strawberries.
Sensitivity analysis, calibration and validation were performed for soil-water EC (EC sw ) at three depth intervals. The goodness-of-fit measures of the model predictions were evaluated using a set of statistical indices including root-mean-squared-error (RMSE), Nash-Sutcliffe efficiency (NSE), coefficient of determination (R 2 ), coefficient of regression (b), and mean relative error (MRE). A perfect fit between observations and model predictions yields a RMSE = 0.0, NSE = 1.0, R 2 = 1.0, b = 1.0 and MRE = 0.0.

Sensitivity Analysis
The sensitivity analyses were performed on the 33 model parameters as indicated in Table 3. We measure the sensitivity of the root mean squared error (RMSE) metric, calculated using the sampled soil water salinity (EC sw ) at the three depth intervals to ensure that our sensitivity indices are grounded relative to the observed soil salinity. The sample sizes and corresponding number of model evaluations required for both the Elementary Effect (EE) and Sobol' methods are listed in Table 4. For the EE method, a sample size of n = 1000 with 10 optimal trajectories were used resulting in 340 model evaluations.
Since the Sobol analysis was performed with the reduced parameter space results determined from the prior EE approach, the number of samples selected ranged from 10 to 12 and iterations were discontinued at the point where the sensitivity results converged. The Latin Hypercube sampling design was used with n = 50 sampling points to generate a near-random sample of parameter values as proposed by Reference [56]. The open-source implementations were used for both methods [48]. Convergence was considered acceptable when within the 95% confidence interval. Among the 33 rootzone model parameters (Table 3), the EE method revealed that only nine strongly influenced the output soil-water EC (EC sw ) as summarized in Table 5. Table 5 indicates parameter sensitivity values µ, µ*, µ*_conf and Σ/EE from the elementary effect analysis. The mean of the distribution of EEi (Σ/EE) is a proxy for a total sensitivity index of i th input factor and µ*_conf is the confidence interval of the µ* at 95%. We chose these nine parameters based on the combination of small µ* values with corresponding large σ values (Figure 5). High σ value indicates that the EE are strongly affected by the choice of the sample points at which they are computed, therefore a non-negligible interaction with other factor values as illustrated in Figures 5 and 6.
Next, we applied Sobol's sensitivity analysis to determine factor fixing using the nine influential factors from the EE results as summarized in Table 6, we reported the first-order (Si) and total order sensitivity (S Ti ). Of the nine factors, only two were non-influential by this analysis, that is, they resulted in S T = 0. These non-influential parameters towards average root zone salinity included the capillary fringe height (Hd) and depth to groundwater (Hwt), so these were fixed at average values for the model calibration analysis. Interestingly, S T values for Tsrz for layer 2 and layer 3 indicate that saturated soil moisture content is more influential for these layers than Ks, and especially that Tsrz in layer 2 is the most influential parameter in the model. Given that we used the control site for model calibration and irrigation of vegetables for this site is managed with sprinklers for establishment and then drip for the rest of the plant development stages, we suspect that plant water uptake is mainly from the bottom layers. However, in our model plant water uptake even with drip is assumed to be mainly from the top layer: 60-30-7-3% uptake pattern. Root water uptake (RWU) although included in the sensitivity analysis was not influential. Mostly likely, upward flow of soil water from the second layer provides water required for uptake in the top layer. Soil moisture is directly related to unsaturated conductivity (a driving force for soil water from wet to dry soil layer). As such the water saturated water content in the second layer ended up being most influential in our calibration. It is important to note that with respect to the crop rooting depth, it is likely that the lettuce rootzone depth is important for the control site specifically as it was the main crop grown for the majority of the experiment period. For other sites, it may be important to include variability of the crop rooting depth.   Next, we applied Sobol's sensitivity analysis to determine factor fixing using the nine influential factors from the EE results as summarized in Table 6, we reported the first-order (Si) and total order sensitivity (STi). Of the nine factors, only two were non-influential by this analysis, that is, they resulted in S T = 0. These non-influential parameters towards average root zone salinity included the capillary fringe height (Hd) and depth to groundwater (Hwt), so these were fixed at average values for the model calibration analysis. Interestingly, ST values for Tsrz for layer 2 and layer 3 indicate that saturated soil moisture content is more influential for these layers than Ks, and especially that Tsrz in layer 2 is the most influential parameter in the model. Given that we used the control site for model calibration and irrigation of vegetables for this site is managed with sprinklers for establishment and then drip for the rest of the plant development stages, we suspect that plant water uptake is mainly from the bottom layers. However, in our model plant water uptake even with drip is assumed to be mainly from the top layer: 60-30-7-3% uptake pattern. Root water uptake (RWU) although included in the sensitivity analysis was not influential. Mostly likely, upward flow of soil water from the second layer provides water required for uptake in the top layer. Soil moisture is directly related to unsaturated conductivity (a driving force for soil water from wet to dry soil layer). As such the water saturated water content in the second layer ended up being most influential in our calibration. It is important to note that with respect to the crop rooting depth, it is likely that the lettuce rootzone depth is important for the control site specifically as it was the main crop grown for the majority of the experiment period. For other sites, it may be important to include variability of the crop rooting depth.

Calibration and Validation
Model calibration was performed allowing the parameters k1, k2, ZrL, Tsrz, Ksrz, Tfcrz, and Twprz identified as most sensitive to model determination of soil-water EC (EC sw ) to vary within their ranges and output compared to the measured soil salinity at the control site. Model validation was completed by simulating EC sw for site 3. The inclusion of the saturated soil water content and saturated hydraulic conductivity in the calibration was crucial as infiltration rates were expected to vary with changes in exchangeable sodium in the soil. O'Geen [57] provides classification of salt-affected soils based on trends of soil water EC, exchangeable sodium percentage (ESP), and SAR. We assessed the potential infiltration problems caused by irrigation water quality following Reference [6] (p. 44) and found that slight-to-moderate reduction in infiltration rates due to irrigation water salinity were expected for the control site and site 4 ( Figure A2). It is however interesting to note that blending well water and recycled water alleviated the possible adverse effects of well water on soil infiltration rates.
The latin hypercube sampling design was used with N = 50 sampling points as proposed by Reference [56]. Intervals were sampled without replacement to ensure even distribution of points with respect to each variable. We executed the model 50 times and computed the corresponding RMSE associated with model predictions. An open-source global optimization code DEoptim written in R was used to find a global minimum RMSE [58]. DEoptim implements the differential evolution algorithm for global optimization. The estimated best fits with the least RMSE values are listed in Table 7. Summarized in Table 8 are the indices associated with comparisons between EC sw measured in the field and model predicted values at different soil depths. For all depth intervals, the sum of first-order effects and sum of total order indices is greater than one, indicating that there are interactions among model factors. Moreover, for factors that have total indices greater than their first-order values, other factors are taking part in the interaction such that throughout the soil profile the root zone hydraulic parameters, rooting depth and dissolution rate are taking part in determination of the soil-water EC, with the dissolution rate accounting for the largest fraction of output variance. This observation provides insight about how well water flow and salinity is modeled.   For all calibrated model runs, regressions of the predicted vs. observed values resulted in a non-zero intercept, b of nearly 1 dS/m. Taken together with the low R 2 -values, we conclude that the model persistently underestimates EC sw , especially those observed EC sw values greater thañ 2 dS/m. Based on the relatively small MRE values that provide an indication of the magnitude of the error relative to observed values without considering the error direction, the model captures salinity dynamics for all layers in the root zone. However, the Figure 8 plot of residuals vs. predicted EC sw for the calibration and validation results exhibits heteroscedasticity, that is, residuals grow as the predicted EC sw values increase. Overall, this latter observation suggests that although the NSE, RMSE, and MRE statistics show that the model has some predictive capacity, it does not capture some processes apparently involved in the soil salinity dynamics. The possible explanation for the differences in measured and predicted soil salinity is that the model does not account for fertilizer and soil amendment management, or plant root uptake of solutes and fertilizer. Generally, transformation (e.g., dissolution) in the soil of different chemicals added during fertigation will increase soil salinity. For example, urea is converted to ammonium that is adsorbed in the soil depending on sol temperature; ammonium is converted to nitrate by nitrification that depends on soil temperature, soil moisture, pH and oxygen content, and nitrate is highly mobile but ammonium, potassium, and phosphorus remain relatively immobile in the root zone. Although the salt index (SI) based on equivalent units of sodium nitrate (developed in 1943 to evaluate the salt hazard of fertilizers) alone cannot be used to evaluate the effect of increased soil salinity from fertilizer applications, it can be used as indicator for the long-term effects on soil salinity. The most commonly used fertilizers in the study region (from California Department of Food and Agriculture annual reports) were nitrogen fertilizers that included urea ammonium nitrate solution (SI = 95), ammonium nitrate (SI =102) and calcium nitrate (SI = 53); phosphorus fertilizer (SI = 7.8-29), potassium sulfate (SI = 46), gypsum, and lime. Sodium nitrate was arbitrarily set at 100, where EC of 0.5 to 40 mass percentage of sodium nitrate is 0.54 to 17.8 dS/m and for a mixture of materials it is reasonable to assume EC is additive for horticulture. As such, the lower the index value the smaller The possible explanation for the differences in measured and predicted soil salinity is that the model does not account for fertilizer and soil amendment management, or plant root uptake of solutes and fertilizer. Generally, transformation (e.g., dissolution) in the soil of different chemicals added during fertigation will increase soil salinity. For example, urea is converted to ammonium that is adsorbed in the soil depending on sol temperature; ammonium is converted to nitrate by nitrification that depends on soil temperature, soil moisture, pH and oxygen content, and nitrate is highly mobile but ammonium, potassium, and phosphorus remain relatively immobile in the root zone. Although the salt index (SI) based on equivalent units of sodium nitrate (developed in 1943 to evaluate the salt hazard of fertilizers) alone cannot be used to evaluate the effect of increased soil salinity from fertilizer applications, it can be used as indicator for the long-term effects on soil salinity. The most commonly used fertilizers in the study region (from California Department of Food and Agriculture annual reports) were nitrogen fertilizers that included urea ammonium nitrate solution (SI = 95), ammonium nitrate (SI =102) and calcium nitrate (SI = 53); phosphorus fertilizer (SI = 7.8-29), potassium sulfate (SI = 46), gypsum, and lime. Sodium nitrate was arbitrarily set at 100, where EC of 0.5 to 40 mass percentage of sodium nitrate is 0.54 to 17.8 dS/m and for a mixture of materials it is reasonable to assume EC is additive for horticulture. As such, the lower the index value the smaller the contribution the fertilizer makes to the level of soluble salts. Thus, fertilizer applications likely add to soil salinities exacerbating the problem over time.
Although the model underestimates EC sw , it adequately captures salinity trends in the leaching of salt during winter months and an increase in salinity water applications and ET during the crop season.
In an effort to evaluate performance of transient vs. steady state models, Reference [12] concluded that the transient models better predict the dynamics of the chemical-physical-biological interactions in an agricultural system. However, since we account for irrigation water salinity, rainfall salinity and dissolution of salts in the soil and exclude additions of fertilizer and soil amendment our simulated EC sw values can be viewed as a likely lower bound of soil salinity associated with the irrigation and farm management practices considered in the model description.
Another complexity possible affecting model prediction is the spatial distribution of salinity with drip irrigation as noted by Reference [59]. They used the transient Hydrus-2D model to compare results between field experiments having both drip and sprinkler irrigated processing tomatoes under shallow water table conditions for a wide range of irrigation water salinities. Both field and model results showed that soil-wetting patterns occurring under drip irrigation caused localized leaching which was concentrated near the drip line. In addition, a high-salinity soil volume was found near the soil surface that increased with increasing applied water EC. Overall, localized leaching occurred near the drip line while soil salinity increased with increasing distance from the emitter and with increasing soil depth. Such localized non-uniformities in leaching are not captured in the one-dimensional model but may have affected field soil sampling in the drip-irrigated fields of our study region. That is, soil samples collected some distance from where dripline emitters were previously operating would likely have greater salinities than would otherwise occur under the uniform leaching and dissolution conditions assumed in the model. Nonetheless, it is important to note that this model is user-friendly and less data intensive and it can be very useful for setting reference benchmarks of long-term salinity impacts of using saline water for irrigation.

Long-Term Salinity with Treated Wastewater Irrigation
We used the calibrated and validated model to simulate the long-term (50-year periods) soil salinity in the fields irrigated with varying fractions of treated wastewater, that is, we applied the model to control site and sites 2 to 7 (Tables 2 and 3). Fifty-year simulations assumed randomly selected rainfall and ET o data from historical records (1983 to 2014) from 2013 to 2049 and the 13-year cropping patterns and irrigation management. In the model calibration, the groundwater table height (Hwt) and groundwater salinity (ECgw) were found to be non-influential parameters with respect to the measured soil water EC (EC sw ). As such average values of measured Hwt and ECgw from the monitoring well located~5 km west of the study site were used, that is, 0.95 m and 0.97 dS/m, respectively. For each simulated case, the three output variables of interest were averaged root zone salinity over the growing season expressed as EC eS (dS/m), annual drainage salt output load as S d (kg/ha/year), and annual runoff salt output load S r (kg/ha/year).
Values for the annual average water-balance terms over the 13-water year simulation at each site are summarized in Table 9a and for the 50-water year simulation in Table 9b. Somewhat greater actual crop water uptake (ET) is achieved under drip and sprinkler irrigated fields with similar crop rotations (sites 2 and 3 as compared to 6 and 7). Leaching fractions were generally very low for all sites at less than 2%. The greatest leaching occurred in the vegetable-strawberry rotation that was first sprinkler than drip irrigated, while the other irrigation and cropping practices yielded similar leaching volumes with the exception of no leaching from drip-irrigated artichokes. Similarly, surface runoff is smaller from drip or sprinkler irrigated fields as compared to sprinkler-furrow irrigated fields. Simulated average (50-year) annual root zone soil water salinity for all sites is shown in Figure 9. Sites managed with sprinkler or drip irrigation and with higher salinity water (EC w ) had higher estimated annual average EC sw , for example, at sites 3, 4, and 5, EC sw was 2.94 dS/m, 3.15 dS/m, and 3.44 dS/m, respectively. For sites 6 and 7 that used sprinklers for germination then furrow irrigation, the latter site received water with higher salinity than the former 1.37 dS/m compared to 1.1 dS/m on vegetable and strawberry rotation, but the resulting average annual soil water salinity differed little, that is, 2.05 versus 2.89 dS/m, respectively. The control site irrigated with well water had the least annual average rootzone soil-water salinity. Furthermore, soil-water salinity equilibrium EC sw ≤ 2.0 dS/m was reached throughout the 50-year horizon for the control site irrigated with well water and after 12 years of irrigated with blended wastewater for sites 2, 6, and 7, whereas for sites 3, 4, and 5 soil water EC increased above 2.0 dS/m in the simulation period. Using Mann-Kendall analysis we found that actual ET had a positive and significant association whereas irrigation amounts had a negative and significant association with EC sw ≤ 2.0 dS/m (Tau = 0.321, p-value = 0.016 and Tau = −0.268, p-value = 0.046 respectively). The Mann-Kendall Tau values indicate the strength and direction of monotonic trends, with −1 and 1 representing perfectly negative and positive monotonic trends, respectively, while the p-value indicates relative significance [60]. Crops are generally assumed to respond to seasonal-averaged root zone salinity of the saturated paste (ECeS) and yield loss thresholds and rates of decline with increasing salinity have been determined based on salinity thresholds in [61]. We calculated ECeS for all the sites and plot the range of ECeS in Figure 10. The maximum seasonal-averaged saturated paste EC for each site was 0.19, 0.27, 0.20, 0.23, 0.24, and 0.26 dS/m for sites 2, 3, 4, 5, 6, and 7, respectively. These values are less than half that of the lowest Mass-Hoffman threshold value of EC e * =1.0 dS/m associated the most sensitive crop (strawberry) in the rotations considered. As such, it is unlikely that long-term irrigation with treated wastewater in the region will adversely affect crop yields significantly. As Platts and Grismer [11] found that salt leaching to deeper soil layers occurred during the rainy season (October-March), while during the growing season soil-water EC increases in soil layers near the surface due to evapo-concentration, on the other hand, applied water salinity causes soil-water EC spikes during the growing season. Rainfall was important towards salt leaching from the root zone at all sites as evident during the wet years 2010, 2016, 2025, 2026, 2038, and 2046 when average annual soil water EC decreased. With the exception of sites 3, 4, and 5, there were no upward trends in soil water salinity over the 50-year period. Overall, relatively constant soil-water EC after 50 years simulation of EC sw < initial EC sw of 2.19 dS/m for all sites except site 5 suggest that there was adequate soil leaching in the region for sustained use of the treated wastewater for irrigation. However, the question remained as to what level of soil salinity would be acceptable especially for annual strawberry production.
Crops are generally assumed to respond to seasonal-averaged root zone salinity of the saturated paste (EC eS ) and yield loss thresholds and rates of decline with increasing salinity have been determined based on salinity thresholds in [61]. We calculated EC eS for all the sites and plot the range of EC eS in Figure 10. The maximum seasonal-averaged saturated paste EC for each site was 0.19, 0.27, 0.20, 0.23, 0.24, and 0.26 dS/m for sites 2, 3, 4, 5, 6, and 7, respectively. These values are less than half that of the lowest Mass-Hoffman threshold value of EC * e = 1.0 dS/m associated the most sensitive crop (strawberry) in the rotations considered. As such, it is unlikely that long-term irrigation with treated wastewater in the region will adversely affect crop yields significantly. In terms of possible adverse environmental effects associated with salinization of surface and ground waters in the region, we determined the cumulative salt output load with deep percolation (Sd) and salt output load with runoff (Sr) during the 50-water year simulation period for the different sites as shown in Figure 11. Salts accompanying surface runoff pose a larger threat in the watershed as these are an order-of-magnitude greater than the cumulative salt output loads to groundwater. Salt loading with deep percolation for the 50-year simulation range up to 3,377 kg/ha with the greatest loads from site 3 and minimal loading for sites 4 and 5 (40 kg/ha and 49 kg/ha respectively) on which annual artichoke crops were grown. Cumulative salt loads accompanying runoff ranged from 19,918 to 59,552 kg/ha, with the greatest loading emanating from site 3 and least from site 4. In comparison to the control site irrigated with well water cumulative salt output loading with deep percolation was 1325 kg/ha and salt output load with runoff 21,505 kg/ha.
To clarify what factors were key to affecting adverse environmental salinization within the region, we tested the effect of applied water EC (ECw) and depths, rainfall depths, actual crop ET and potential crop ET, and number of days fallow on soil water EC (ECsw) and salt output loads with runoff (Sr) or deep percolation (Sd) using the non-parametric Mann-Kendall trend analysis (the 'Kendall' test in the R package [60]. With respect to the Mann-Kendall Tau values (Table 10), the applied water EC (ECw) and rainfall depths had positive and significant effects on annual average soil water EC (ECsw) and annual salt output loads with runoff (Sr) and deep percolation (Sd). Calculated actual crop ET had a positive and significant effect on ECsw; this is expected as water uptake by the plant and evaporation leave salts behind. On the other hand, actual crop ET had a negative and significant effect on Sd, and had no significant effect on Sr. Applied water depths had positive and significant effect on Sd and Sr. The number of days fields were fallowed had a negative effect on ECsw and a positive effect on Sd and Sr. Overall, these observations were consistent with that expected from the field observations and described above. In terms of possible adverse environmental effects associated with salinization of surface and ground waters in the region, we determined the cumulative salt output load with deep percolation (S d ) and salt output load with runoff (S r ) during the 50-water year simulation period for the different sites as shown in Figure 11. Salts accompanying surface runoff pose a larger threat in the watershed as these are an order-of-magnitude greater than the cumulative salt output loads to groundwater. Salt loading with deep percolation for the 50-year simulation range up to 3,377 kg/ha with the greatest loads from site 3 and minimal loading for sites 4 and 5 (40 kg/ha and 49 kg/ha respectively) on which annual artichoke crops were grown. Cumulative salt loads accompanying runoff ranged from 19,918 to 59,552 kg/ha, with the greatest loading emanating from site 3 and least from site 4. In comparison to the control site irrigated with well water cumulative salt output loading with deep percolation was 1325 kg/ha and salt output load with runoff 21,505 kg/ha.
To clarify what factors were key to affecting adverse environmental salinization within the region, we tested the effect of applied water EC (EC w ) and depths, rainfall depths, actual crop ET and potential crop ET, and number of days fallow on soil water EC (EC sw ) and salt output loads with runoff (S r ) or deep percolation (S d ) using the non-parametric Mann-Kendall trend analysis (the 'Kendall' test in the R package [60]. With respect to the Mann-Kendall Tau values (Table 10), the applied water EC (EC w ) and rainfall depths had positive and significant effects on annual average soil water EC (EC sw ) and annual salt output loads with runoff (S r ) and deep percolation (S d ). Calculated actual crop ET had a positive and significant effect on EC sw ; this is expected as water uptake by the plant and evaporation leave salts behind. On the other hand, actual crop ET had a negative and significant effect on S d , and had no significant effect on S r . Applied water depths had positive and significant effect on S d and S r . The number of days fields were fallowed had a negative effect on EC sw and a positive effect on S d and S r . Overall, these observations were consistent with that expected from the field observations and described above.

Discussion
We calibrated and validated a modified root zone salinity model originally developed by Isidoro and Grattan [17], which was then applied to estimate long-term soil salinity in fields irrigated with treated wastewater. We conducted a global sensitivity analysis using the elementary effect/Morris and Sobol's methods to first reduce the number of influential model parameters important to calibration and that need to be acquired from the field. Seven of the thirty-three model parameters were found to be critical to root zone soil salinity dynamics. These were parameters accounting for salt dissolution in the soil, root zone hydraulic parameters, and crop rooting depth. Model calibration resulted in a satisfactory fit to the observed field data; however, the model underestimated soil water  and that need to be acquired from the field. Seven of the thirty-three model parameters were found to be critical to root zone soil salinity dynamics. These were parameters accounting for salt dissolution in the soil, root zone hydraulic parameters, and crop rooting depth. Model calibration resulted in a satisfactory fit to the observed field data; however, the model underestimated soil water salinity (EC sw ), especially for large EC sw > 2 dS/m measured during the growing season. We attributed this error to the model's failure to account for fertilizer and soil amendment applications and transformation thereof (e.g., gypsum dissolution) in the soil. In addition, drip irrigation leads to very localized variations in soil salinity that depend on the distance from the emitter that are not considered in the model and may have affected field soil sampling between plantings and harvests. Nonetheless, the model adequately captured soil-water EC trends that were congruent with observed data. Sites irrigated with greater salinity water (EC w ) combined with sprinkler or drip had greater estimated annual average soil-water salinity (EC sw ). Sites that combined sprinkler irrigation for germination with furrow for the remaining development stages resulted in lower annual average EC sw in the root zone even when more saline irrigation water was applied. Rainfall played an important role in the leaching of salts from the root zone as during wet years average annual soil water EC decreased at all sites. Moreover, rainfall had a negative and significant effect on annual average root zone EC sw . We found that for all sites use of treated wastewater for irrigation over the 50-year period does not affect strawberry yields, the most salt sensitive of the crops in the rotations encountered. Overall irrigation water EC (EC w ), rainfall amounts, actual calculated crop ET and the number of days fields were fallowed had significant effects on annual average soil water EC (EC sw ). EC w , rainfall and actual crop ET effects were positive and fallowing decreased average root zone EC sw . On the other hand, irrigation amounts and number of days fallowed had positive and significant effects on salt output loads associated with runoff and deep percolation. Moreover, soil water salinity equilibrium EC sw ≤ 2.0 dS/m is reached throughout the 50-year horizon for the control site irrigated with well water and after 8, 9 and 14 years of irrigated with blended wastewater for sites 2, 6, and 7 respectively. For sites 3, 4, and 5 soil water EC increased above 2.0 dS/m in the simulation period. Actual ET had a positive and significant association whereas irrigation amounts had a negative and significant association with ECsw ≤ 2.0 dS/m.
While we believe that the modeling results can inform recommendations about irrigation management practices and for estimating salt output loading resulting from use of saline waters for irrigation, difficulties in linking field observations of soil salinity and model predictions remain troubling. However, since we account for irrigation water salinity, rainfall salinity, and dissolution of salts in the soil and exclude additions of fertilizer and soil amendment our simulated EC sw values are likely a lower bound of soil salinity associated with the irrigation and farm management practices considered in the modeling. Nonetheless, it is important to note that this model is user-friendly and less data intensive and it can be very useful for setting reference benchmarks of long-term salinity impacts of using saline water for irrigation.

Unsaturated Soil Water Movement
The soil matric potential (ψ) was related to the volumetric water content (θ) by means of Equation (A1) [32]: Where θ s was the volumetric water content at saturation, ψ s is the water entry potential or "saturation" water potential and b is the slope of the water retention curve on a logarithmic plot. For each soil type, b and ψ s were calculated from the volumetric water content at field capacity and wilting point and their respective potentials in absolute value (ψ FC = 316 cm and ψ WP = 15,849 cm; so that pF (FC) = 2.5 and pF (WP) = 4.2). Taking logarithms, the expression of the potentials for FC and WP become linear equations: log(ϕ FC ) = 2.5 = logϕ s − b. log( θ FC θ s ) log(ϕ WP ) = 4.2 = logϕ s − b. log( θ WP θ s ) (A2) from which, b and ψ s are estimated. The unsaturated hydraulic conductivity for a given θ was given by: where K s is the saturated hydraulic conductivity. Thus, unsaturated flow between layers (U) can be calculated as: where ∆Z is the center to center simulation distance selected between layers and neglecting the gravitational gradient.

Crop Water Uptake
Non-stressed crop ET is calculated as: where K c is the crop coefficient and varies with the crop development stages (Table A2) and ET o is the reference ET. Between cropping seasons, all ET or evaporation E was assumed to take place from the upper layer. For this period Kc was calculated from the mean interval between precipitation events of each month and the mean precipitation event in each month and ET o [38].
In each layer (k), the actual crop ET can be lower than ET c (k) due to water stress, which depends on the soil water content and the sensitivity of the crop to low water contents, accounted for through the crop-specific parameter p: the ratio of readily available soil water (RAW) to total available water (TAW) (p = RAW/TAW) [38]. When the soil water content (W(k)) in a layer fell below We(k) = WP + (1 − p) TAW, the ET from that layer actual crop ET(k) dropped below the ETc(k), and the actual ET of the layer was calculated as: actual crop ET = K s ×ET c (A6) where K s is a stress coefficient [38]: when one layer was stressed during the growing season (W(k) < W e (k)), the model allowed increase in the extraction coefficient of the lower layer to supply the ET demand of the day. The root zone water uptake pattern depends on irrigation frequency. Root uptake patterns were taken from [18,42,43]. Root length increase a function time is calculated as [45]: where L z is the rooting depth at time t, L o is the starting root depth, L max is the maximum root length, t L max time after planting when L max is reached and t o is time to reach 90% crop emergence. This is a linear root expansion; the method assumes that once half of the time required for crop emergence is passed by t o 2 , the rooting depth starts to increase from an initial depth L o till L max is reached.

Water Balance
Surface runoff for winter rainfall and a fraction of applied water is modelled using the SCS method. We define curve number (CN) associated with row crop cover for the growing season and bare soil for non-growing season from the SCS tables and calculate precipitation runoff as: where P is runoff producing precipitation and S is the potential maximum retention after runoff begins related to CN by: Daily water balance for the 4-layered root zone and 2-layered vadose zone is performed. To account for the slow water movement between layers for low water content below field capacity, a slow upward or downward flow U is calculated dependent upon the difference in matric potential between soil layers (Equations B4). In the first quarter of the root zone inflows and outflows include applied water (I) and rainfall (P), the drainage above field capacity (D (1)) to layer 2, actual crop ET and U. For the underlying root zone layers, inflows and outflows include drainage (D) from the overlying layer, U and actual crop ET and finally for the unsaturated layers below the root zone inflows and outflows include D and U.
When the soil water content in layer "k" is above field capacity, the excess water drains to the lower layer over a two-day period, the higher flow in the first day than the second. The fraction α of the excess water that drains the first day is calculated from the soil texture in the layer through an empirical relation obtained to match results presented by approximately 0.9 for sand, 0.85 for loam and 0.7 for clay [49]. Two arbitrary water contents were defined from field capacity to saturation for each layer W a and W b defined as: Drainage (D) is calculated as: