Dynamic assessment of the impacts of global warming on nitrate losses from a subsurface-drained rainfed-canola field

The impact of global warming on water and nitrate losses from a rainfed-canola cropping system under various artificial drainage systems was assessed using an integrated field-modeling approach. Four subsurface drainage systems with different drain depths (Dx) and spacings (Ly), including D0.90L30, D0.65L30, D0.65L15, and Bilevel (with a drain spacing of 15 m and alternate drain depths of 0.65 and 0.90 m), were considered. The HYDRUS (2D/3D) model was first calibrated and validated using data collected for all drainage systems during the 2015–2016 and 2016–2017 canola cropping cycles, respectively, and then applied to simulate water/nitrate losses for different drainage systems under meteorological conditions predicted assuming expected future global warming. Future weather data were downscaled from 20 general circulation models and four RCP scenarios for the mid 21st century (for 2041–2070). The model capability of representing experimental field data was evaluated using the mean bias error (MBE), the normalized root mean square error (nRMSE), and the model efficiency (EF). The HYDRUS (2D/3D) model provided reliable description of soil water contents (MBE=-0.5 % to 0.2 %, nRMSE = 0.005−0.034%, and EF = 0.73−0.99), drainage fluxes (MBE= -21.7 × 10 to 24.9 × 10 mm d, nRMSE = 0.23−0.37%, and EF = 0.69−0.85), soil nitrate concentrations (MBE= -0.002 to 1.00 mg cm, nRMSE = 0.08−0.18%, and EF = 0.51−0.88), and nitrate fluxes (MBE= -0.97 to 0.72 mg cm d, nRMSE = 0.35−0.57%, and EF = 0.77−0.87). The modeling results indicate that climate change will cause an increase of up to 148 % in average daily drainage fluxes and up to 125 % in average daily nitrate fluxes compared to the base case. This will result in an increase of 4–125 % in seasonal nitrate losses from various drainage systems, with the lowest and highest projections for the D0.65L15 and D0.65L30 systems, respectively. The HYDRUS-simulated results indicate that the D0.65L15 system is environmentally safer than the other evaluated drainage systems for predicted global warming conditions concerning water/nitrate losses.


Introduction
The need to feed the growing global population on the one hand, and increasing global scarcity of blue (fresh) water (WWAP, 2009;Hoekstra et al., 2012) on the other hand, indicate that it is essential to expand rainfed agriculture in the world. In this regard, humid regions may be of higher importance since they receive a sufficient quantity of precipitation to fully supply crop's water requirements by green (from rainfall) water (Shahsavari et al., 2019). However, expanding dry-land farming in these regions may be restricted by waterlogging problems, which occasionally occur after heavy rainfalls and, particularly, in heavy soils (Darzi-Naftchali et al., 2017). Waterlogging threatens the year-round cropping, resulting in considerable areas either going out of production or experiencing reduced yields (Darzi-Naftchali et al., 2013). Under such circumstances, installing subsurface drainage systems may help in providing suitable conditions for winter-crops-based cropping systems and consequently improving the annual productivity of these lands. Subsurface drainage systems can speed up the water table drawdown and provide better aeration during the growing season.
While being beneficial in terms of the water table control, installed drainage systems may cause serious environmental challenges due to their negative impacts on nutrient leaching. Previous researchers have demonstrated that drainage systems may accelerate the leaching process of soil nutrients, and in particular of nitrogen (N) (Kalita et al., 2006;Furukawa et al., 2008;Zhang et al., 2007;Darzi-Naftchali et al., 2017), which is well known to be an essential crop nutrient affecting crop growth and yield (Wienhold et al., 1995;Jia et al., 2014). The design parameters of drainage systems, including the type of the systems and the drain depth and spacing, may affect water or N losses from agricultural lands (Christen and Skehan, 2001;Wahba and Christen, 2006;Jafari-Talukolaee et al., 2015Darzi-Naftchali et al., 2017). For certain drainage systems, the rate of water and N losses additionally also depends on field management practices. Improving field management practices may increase the water/nutrient holding capacity of the soil and thus reduce water and N losses (Azooz and Arshad, 1996;De Vita et al., 2007;Triplett and Dick, 2008;Udayasoorian et al., 2009;Constantin et al., 2010;Mitchell et al., 2012;Liu et al., 2013;Phogat et al., 2013;Chukalla et al., 2017).
Drainage fluxes represent a dominant factor for water and N losses from drained areas. Among various influencing factors, climatic variables, including precipitation (P) and potential evapotranspiration (ET 0 ), play key roles in determining the rate of drainage fluxes. Increased surplus precipitation, defined as P-ET c , may accelerate the rate of water and N losses due to enhanced drainage fluxes. Water and N losses from drained dry-farming lands may become a more important concern in the future when predicted global warming, which affects these climatic variables, takes place. Numerous researchers have demonstrated various ranges of both negative as well as positive projections of climate change into P and ET 0 in different parts of the world (Abbaspour et al., 2009;Dastoorani and Poormohammadi, 2012;Terink et al., 2013;Agarwal et al., 2014;Kazemi-Rad and Mohammadi, 2015;Karandish and Mousavi, 2018;Adham et al., 2019;Darzi-Naftchali and Karandish, 2019). Through such investigations, earlier research has mainly focused on determining probable economic consequences of global warming on the agricultural sector due to the climate change impact on crop yields and productivity under different future greenhouse gas emission scenarios Harmsen et al., 2009;Massah Bavani and Morid, 2005).
Predicting probable N losses under climate change conditions is essential because the byproduct of these losses is the pollution of water resources (Karandish and Šimůnek, 2017) due to N leaching/drainage from agricultural lands (Zhu et al., 2005;Thompson et al., 2007;Dudley et al., 2008;Burow et al., 2010;Dahan et al., 2014;Darzi-Naftchali et al., 2017). This issue has higher importance in humid regions, where a low recovery of N fertilizers by crops produces excessive N losses after heavy rains (Karandish and Šimůnek, 2017). Although a few attempts have been made to investigate the effects of climate change on water and N losses in agricultural systems (e.g., Singh et al., 2009;Wang et al., 2015;He et al., 2018;Jiang et al., 2020), a literature review reveals that no such study has been carried out for subsurface-drained paddy fields under winterrainfed cropping. As such systems are in their infancy in northern Iran paddy fields, a comprehensive study of their behavior would enable policymakers to incorporate the effects of inherent complexities of climate change on drainage activities to protect the local environment.
Therefore, we carried out a field experiment to investigate water and N losses from a drained field cultivated with rainfed-canola and with various subsurface drainage systems. Our research is novel since it is the first study in which environmental hazards in rainfed agriculture due to global warming are evaluated in terms of N losses from the field. To the best of our knowledge, our study assesses for the first time such hazards in fields with subsurface drainage systems. The N dynamics under various climate change projections were assessed using the HYDRUS (2D/3D) model (from now on referred to simply as HYDRUS), which was first calibrated and validated using the field-collected data. We used the HYDRUS model since its capability in capturing drainage and N fluxes has been previously confirmed by many researchers (e.g., Öztekin, 2002;Mirjat et al., 2014;Filipovic et al., 2014;Li et al., 2014;Mguidiche et al., 2015;Karandish and Šimůnek, 2016;Mekala and Nambi, 2016;Matteau et al., 2019). The main objectives of this research thus were (i) to calibrate and validate HYDRUS using experimental drainage and nitrate fluxes collected in a subsurface-drained field cultivated with rainfed canola, and (ii) to investigate the impact of climate change projections on key climatic variables, drainage, and nitrate fluxes.

Materials and methods
2.1. Field trial 2.1.1. Study area and drainage system layout The field research was carried out during two canola growing seasons (2015−16 and 2016−17) in the subsurface drainage pilot of the Sari Agricultural Sciences and Natural Resources University in the   Haghnazari, et al. Agricultural Water Management 242 (2020) 106420 Mazandaran province, Iran (SANRU: 36.3• N, 53.04• E; 15 m below sea level) (Fig. 1). In 2011, this pilot was designed and implemented in a 4.5-ha consolidated paddy field (Darzi-Naftchali et al., 2013Darzi-Naftchali and Ritzema, 2018), consisting of four subsurface drainage systems with different drain depths (D x ) and spacings (L y ), including D 0.90 L 30 , D 0.65 L 30 , D 0.65 L 15 , and Bilevel which has a drain spacing of 15 m and alternate drain depths of 0.65 and 0.9 m (Fig. 1).
The field consisted of consolidated paddy plots 100 m long and 30 m wide. Considering the dimensions of the plots and the field size, the drainage systems were installed at a scale much larger than the laboratory scale. Drainage installation costs at such a scale were very high and required a large area. Different drainage systems were not replicated due to limited resources and land. The drainage pipes were installed so that the last drain line in each system acted as the first drain line in the adjacent system. Accordingly, at least three drain lines were considered in each drainage system. In addition, the monitoring lines (lines 2, 4, 5, 7, and 9 in Fig. 1) for measuring drain discharge and water quality were selected to fully represent each drainage system by removing the buffering lines. A detailed description of the drainage systems can be found in Darzi-Naftchali et al. (2013. Long term averages of annual precipitation, mean temperature, as well as minimum and maximum temperatures recorded in the study area, are 616 mm and 17.3°C, −6°C, and 38.9°C, respectively . The variation of climatic variables during the 2015-2016 (the calibration period) and 2016-2017 (the validation period) growing seasons is shown in Fig. 2.

Cropping system and data collection
The area was under rice-canola cropping during 2011-2017, with rice as a major crop and canola as a winter crop. The data from the 2015−16 and 2016−17 canola growing seasons were used in the present study. Canola was cultivated in the two cropping cycles on October 3 in 2015 and on September 28 in 2016. Crops were harvested on May 4, 2016, and May 20, 2017. All agricultural practices were the same as the conventional practices of local farmers in the study area. In the growing season of 2015-2016, 50 kg ha −1 of triple superphosphate was applied before cultivation, and 50 kg ha −1 of urea (i.e., CO(NH 2 ) 2 with 46 % N) was applied 52 days after sowing (DAS). In the second growing season, urea was applied at rates of 85 kg ha -1 , 85 kg ha -1 , and 115 kg ha -1 at 24, 96, and 116 DASs, respectively. The growers in the region adjust fertilization based on their activities in the crop rotation. For example, remaining residues after the previous crop are considered by the experienced growers (Darzi-Naftchali et al., 2017).
Before crop cultivation, soil samples were collected every 30 cm to a depth of 200 cm to determine soil physical and hydraulic properties and the soil N-content. Soil water characteristic curves of the soil samples were determined by measuring soil water contents at 14 different pressure heads (varied in the range of 0-15 bars). After that, the RETC model was applied to fit the van Genuchten-Mualem model (van Genuchten, 1980) to the observed retention data. For each treatment, two suction samplers were buried at depths of 30 cm and 60 cm beneath the surface before crop cultivation and remained continuously at 30 kPa. The samplers had porous ceramic caps of 5 cm in diameter. Using the vacuum pump, the leachate was then collected during the growing seasons (once every two weeks). Nitrate concentrations of water samples were determined by spectrophotometer (DR-4000 HACH).
To determine daily variations in the soil water content (SWC) during the cropping cycles, a 100-cm long TDR probe (Trime FM; IMKO; Germany) was installed midway between drains in each subsurface drainage system (i.e., 4 TDR probes were installed in the study area; 1 probe * 4 drainage systems). TDR probes were used once a day to measure SWC at a 15-cm interval during both growing seasons. The accuracy of all TDR sensors was evaluated by comparing TDR-measured SWCs with corresponding values measured using the gravimetric method. TDR-measured SWCs agreed well with gravimetrically-measured SWCs with a coefficient of determination of 95 %.
Water table depths were measured daily in observation wells that were dug out midway between adjacent drains in each drainage system. Free drainage management was adopted during canola growing seasons. Subsurface drainage discharges (q) were measured daily for all drainage systems as long as the flow was observed at the outlet of drains. Drain fluxes were measured by using partial flumes at the outlet of representative drains during drainage periods. The total water discharge for a particular day is the sum of discharges over 24 h. The above values were converted to water depths by dividing them with the plot area. Drainage water samples were taken every 15 days and several consecutive days after fertilization. Collected water samples were then analyzed to determine their nitrate concentrations.

Climate data under global warning
The projections of future climate data were downscaled for four RCP scenarios (the Representative Concentration Pathways reported in the 5th IPCC report, IPCC, 2013): the RCP2.6 (a low greenhouse gases emission scenario), RCP4.5 (an intermediate one), RCP6.5 (a high one), and RCP8.5 (a very high greenhouse gases emission scenario). For each scenario, the projections of air temperature (minimum and maximum),  precipitation, solar net radiation, wind speed, atmosphere pressure, and relative humidity were generated using 20 different GCM models (general circulation models, Table 1) for the future period of 2040-2070, and compared with the corresponding values obtained for the base period . The statistical downscaling approach was carried out based on the widely used change factor method (Jones et al., 1997). For all climatic variables, the change factor coefficients, which determine the relationship between the current and future climatic variables, were calculated based on Eq. 1 (Jones et al., 1997), except for air temperature, for which Eq. 2 was applied (Trzaska and Schnarr, 2014).
wheres P Δ i is a dimensionless monthly change factor for a particular climatic variable (p), P GCM fut i , , and P GCM base i , , are monthly averages of future (simulated) and historical values of a considered climatic variable, respectively, T Δ i is a dimensionless monthly change factor for air temperature, and T GCM fut i , , and T GCM base i , , are monthly averages of future and historical values of air temperature, respectively. These change factors were calculated at a monthly time scale (i = 1-12).
Using generated climatic variables, reference evapotranspiration (ET 0 ) was calculated using the Penman-Monteith FAO-56 equation (PMF) (Allen et al., 1998) due to its global acceptability . ET 0 was then used to estimate crop evapotranspiration as explained below.

N dynamics under current and future conditions
2.3.1. The HYDRUS model 2.3.1.1. Model description and governing equations. Soil water and N dynamics were simulated using the HYDRUS (2D/3D) model (Šimůnek et al., 2008, 2016). The HYDRUS program numerically solves the Richards equation for saturated-unsaturated water flow: where Θ is the volumetric soil water content (SWC) [L 3 L −3 ], K is the unsaturated hydraulic conductivity [LT -1 ], h is the soil water pressure head [L], x is the lateral coordinate [L], z is the vertical coordinate   Table 2 Soil properties, and soil hydraulic and solute transport parameters for the study area (i.e., the calibrated and validated data for soil water and solute transports are obtained through HYDRUS modeling).
where γ(h) is the soil water stress response function (dimensionless) of Feddes et al. (1978), RDF is the normalized spatial root water uptake distribution [L −2 ], T pot is the potential transpiration rate [LT -1 ], and W is the width of the soil surface [L] associated with the transpiration process. The stress response function was obtained from the HYDRUS database for maize (Šimůnek et al., 2011). Uniform root distribution was assumed in each soil layer. The van Genuchten-Mualem constitutive relationships were used to model soil hydraulic properties (van Genuchten, 1980). The transport of nitrate (NO 3 − ) was modeled using the convectiondispersion equation embedded in the HYDRUS model. This equation is reported to be suitable for conservative solutes such as NO 3 where c is the NO 3 − concentration in the liquid phase (ML -3 ), q x and q z are the components of the volumetric flux density (LT -1 ), D xx , D zz , and D xz are the components of the dispersion tensor (L 2 T -1 ) (Bear, 1972), S c is a sink term, which generally includes local NO 3 − uptake (through a passive process), mineralization, microbial immobilization, and denitrification (ML -3 T -1 ). The first term on the right side of Eq. 3 describes the solute flux due to dispersion, the second term describes the solute flux due to convection with flowing water, and the third term describes nutrient uptake by roots. Similarly, as in other earlier studies (e.g., Ajdary et al., 2007; Wang et al., 2010; Tafteh and Sepaskhah, 2012), mineralization gains were neglected due to the following reasons. First, the lack of organic matter (OM) lmits the mineralization process in mineral soils (i.e., soils with OM < 3% in the upper horizon, Huang et al., 2009 cronologically.), as in our study, the lack of OM limits the mineralization process (Li et al., 2003;Deenik, 2006;Wijanarko, 2015). Second, mineralization mainly occurs in coarse-textured soils with low clay content, while it decreases considerably as soil clay content increases (Karandish and Šimůnek, 2017). Third, the abundance of micropores in high clay soils causes a physical protection of OMs from being microbially decomposed and mineralized (Deenik, 2006). Our frequent measurements also showed that NO 3 − concentration in the soil were much higher than NH 4 concentrations (Darzi-Naftchali et al., 2016). Only the NO 3 − transport in the soil was thus considered in this research, assuming that the input of N-fertilizer in the form of urea was instantly nitrified into NO 3 -. This is similar to the assumptions made by many other researchers (e.g., Ajdary et al., 2007;Tafteh and Sepaskhah, 2012;Karandish and Šimůnek, 2017), who assumed that nitrification is faster than the other processes, and nitrifying urea into NO 3 -takes only a few days (Havlin et al., 2006). In addition, nitrification is reported to be stimulated in less-organic (OM < 3%), finely-textured soils with high clay content, all of which are the cases in the current research. We also considered that root N uptake was strictly passive, which has also been assumed in many other studies (e.g., Hanson et al., 2006;Tafteh and Sepaskhah, 2012). According to this assumption, N uptake can be calculated by multiplying the local soil N concentration and root water uptake.
2.3.1.2. Flow domain, and initial and boundary conditions. The impermeable layer in the study area was located at a soil depth of 200 cm (Fig. 3). Hence, the 2D transport domain was defined as a rectangle with a depth of 200 cm and a width of 30 m for D 0.9 L 30 and D 0.65 L 30 drainage systems and 15 m for the D 0.65 L 15 and Bi-level  drainage systems. The drains are considered to be on the sides of the transport domain. An unstructured triangular finite element mesh (FEM) was used to discretize the defined transport domain. A nonuniform FEM was generated by HYDRUS with finite element sizes gradually increasing with distance from the drains. The entire transport domain was divided into six soil layers, i.e., 0−30 cm, 30−60 cm, 60−90 cm, 90−120 cm, 120−150 cm, and 150−200 cm soil depths, for which different soil properties were defined. To represent the high permeability of the backfilled drain trench (gravel), an additional soil layer was defined 10 cm around, 10 cm below, and 30 cm above the drains. Measured soil water contents in different soil layers were considered to represent initial conditions for water flow simulations. The saturated soil water content was specified as an initial condition in soil layers below the water table to differentiate the saturated and unsaturated zones. The atmospheric boundary condition, defined using potential evaporation (E P ), potential transpiration (T P ), and precipitation (P), was applied at the top of the transport domain. The dual-crop coefficient approach (Allen et al., 1998) was adopted to separate crop evapotranspiration (ET C ) into E P and T P as follows: where K C is the crop coefficient, K Cb is the basal crop coefficient, and K e is the evaporation coefficient. The values of K Cb for different crowing seasons were taken from Allen et al. (1998): 0.3, 1.05, and 0.25 for the initial, mid-, and late-growing stages of canola, respectively. The values of K e were then estimated as 0.05, 0.1, and 0.1, for then mentioned growing stages, respectively. A seepage face boundary condition was used to represent the drains. All other remaining boundaries were assigned a no-flow boundary condition. The initial condition for solute transport simulations was defined using the measured soil NO 3 − content in different soil layers. A thirdtype Cauchy boundary condition was used to describe the concentration flux at the top boundary and the drains. A Cauchy boundary condition is automatically converted into a second-type Neumann boundary condition during periods of drain outflow.
2.3.1.3. Calibration and validation. The experimental data were first used to calibrate soil hydraulic and solute transport parameters and then to validate the HYDRUS model. These parameters were optimized using the inverse solution option of HYDRUS. In this process, soil hydraulic parameters, including the saturated hydraulic conductivity (K s ), the saturated soil water content (θ s ), the residual soil water content (θ r ), and solute transport parameters, including the longitudinal (D L , L) and transverse (D T , L) dispersivities, were optimized. Measured drain fluxes (q) and NO 3 − concentrations in the 2015-2016 growing season were used to calibrate the HYDRUS model. The calibrated parameters were then used to validate the model using the same data collected in the 2016-2017 growing season. The molecular diffusion coefficient was always set equal to zero since molecular diffusion in soils can usually be neglected (Radcliffe and Šimůnek, 2010;Karandish and Šimůnek, 2017 where P i and O i are simulated and observed data, respectively, O and P are the averages of observed and simulated data, respectively, and n is the number of observations.

Scenario analysis
The calibrated and validated HYDRUS model was then applied to estimate the probable consequences of climate change projections on drainage water (q) and NO 3 − fluxes for different drainage systems. The climate projections were first generated for the base period (i.e. , 1975-2005), and then for all 20 GCMs for four considered RCP scenarios over the 2040-2070 period. To obtain consistent comparisons, the initial and boundary conditions, as well as all agricultural practices, were considered to be the same as during the validation period (i.e., the 2016-2017 growing season). A comparative analysis was then carried out between the simulated projections for the 2041-2070 period and the base period. In this regard, a relative change in the considered parameter (i.e., climatic variables, drainage flux, NO 3 flux, or seasonal N loss) was calculated as follows: where RC is a relative change in the considered parameter (%), X f and X b are values of the considered parameter in the future and base periods, respectively.

The HYDRUS(2D/3D) model efficiency
The calibrated soil hydraulic and solute transport parameters are summarized in Table 2. Fig. 4 shows temporal variations of the observed and HYDRUS-simulated drain fluxes (q) for different drainage   Table 3, also indicate the strong predictive capability of the model. MBE, nRMSE, and EF ranged from -21.7 × 10 −3 to 24.9 × 10 −3 mm d -1 , from 23 % to 37 %, and from 0.69 to 0.85, respectively, among different drainage systems and different growing seasons. The efficiency of the HYDRUS model to predict drainage fluxes q is also supported by other researchers (e.g., Öztekin, 2002;Mirjat et al., 2014;Filipovic et al., 2014;Li et al., 2014;Matteau et al., 2019).
The quantitative assessment of the model performance, as well as the visual inspection presented in Fig. 5, also shows good agreement between the observed and HYDRUS-simulated soil water contents (SWCs) in different soil layers during the calibration and validation periods; with MBE, nRMSE, and EF ranging from -0.005 to 0.002 cm 3 cm −3 , from 5 to -3.4 %, and from 0.73 to 0.99, respectively. These differences may be a consequence of the fact that the HYDRUS-simulated SWCs are compared with the measured SWCs, which are averaged over a certain soil volume, in which the SWC gradient caused by irrigation/precipitation/drainage may not be linear (Karandish and Šimůnek, 2016;Mguidiche et al., 2015). The capability of the HYDRUS model to predict SWC variations is also supported by other researchers (Ramos et al., 2012;Kandelous et al., 2012). Table 3 also shows good correspondence between the observed and HYDRUS-simulated NO 3 concentrations in different soil layers during both growing seasons (MBE=(-0.002)-(+0.001) mg cm −1 d −1 , nRMSE = 7.9-18.4 %, EF = 0.51−0.88), which indicates that the calibrated model is well suited to describe the flow and transport processes observed in the experimental field.
To ensure the capability of the HYDRUS model to simulate the solute transport adequately, we also compared temporal variations of the observed and simulated NO 3 fluxes in different drainage systems during the calibration and validation periods (Fig. 6). While the HYDRUS model either slightly underpredicted or overpredicted NO 3 fluxes during different periods, the simulated fluxes overall closely matched the observed data for all drainage systems during both the calibration and validation periods, with correlation efficiencies of 0.77−0.9 (data not shown). The model performance criteria reported in Table 3 also indicate the high potential of the HYDRUS model in capturing the NO 3 transport in different drainage systems. MBE, nRMSE, and EF varied in the range of (-97)-(0.72) mg cm −1 d −1 , 35.5-57 %, and 0.77−0.87, respectively, across different drainage systems and during both growing seasons. The model capability of simulating N dynamics under different drainage conditions is also supported by others (e.g., Mekala and Nambi, 2016;Matteau et al., 2019). Fig. 7 shows the range of relative changes in different climatic variables for different RCPs, along with the range of GCMs projections for each scenario. Regardless of the type of the climatic variable, Fig. 7 shows a different behavior of projected variables for different GCMs and RCP scenarios, which could be attributed to different resolutions of the ocean models for different GCMs. Climate processes fully depend on how GCMs simulate the extent of sea ice, sea surface temperature, surface heat, ocean heat transfer, and momentum fluxes (CCSP, 2008;Karandish and Mousavi, 2018). Different projections by different GCMs could also be attributed to many other factors, which control the simulation process, such as the prognostic variable for cloud characterization and the compatibility between the heat and water budgets of the atmospheric and ocean models (Randall et al., 2007). However, these differences among RCP scenarios can be attributed to the embedded assumptions in the socio-economic and environmental models for each scenario. Other researchers have also indicated that future projections of climate variables depend on either the RCP scenarios or the selected GCM (e.g., Adham et al., 2019;Haj-Amor et al., 2020). Climate change projections showed both positive (increases) and negative (decreases) changes in monthly average values of T max and ET 0 , while monthly precipitation P is more likely to increase in the future. Except for RCP2.6 and a few GCMs for RCP4.5, climate change may cause an increase in T min in the study area. Differences among the projections of different GCMs/RCP scenarios indicate the uncertainty in the projections of future climatic variables. Quantifying these ranges of uncertainties is essential since they provide more accurate insights for developing rational plans to cope with plausible consequences of global warming. Fig. 7 shows that the range of uncertainty arising from projections of different GCMs for different RCP scenarios is lower for T min and T max , which indicates more similar predictions of temperatures compared to water fluxes, such as P and ET 0 . Such results are in agreement with those reported by . Fig. 7 also demonstrates non-uniform temporal projections of climatic variables by different RCP scenarios, which is also confirmed by others (e.g., Zickfeld et al., 2005;Girvetz et al., 2009;Harmsen et al., 2009;Agarwal et al., 2014;Adham et al., 2019). Precipitation projections seem to be less uncertain over the October-January period, while smaller uncertainties in the projections of other variables were observed in May. Fig. 7 shows that compared to the base period, and based on the average value obtained from all GCMs under four RCPs, T min may increase by 0.15-1.25°C, with the lowest and highest increases in April   Fig. 7. The range of relative changes in minimum (T min ) and maximum (T max ) temperatures, monthly precipitation, and monthly ET 0 in the study area for different future RCP scenarios (i.e., 2041-2070) compared to the base period (i.e., 1975-2005). The lower and upper ends indicate the 5% and 95 % intervals of the uncertainty ranges, and dots show the 50 % value. and October, respectively. T max may increase by 0.45-2.18°C, except in April and May, when T max decreases by 0.08°C and 0.15°C, respectively. Others have reported a general increase in both T min and T max under climate change in the study area (e.g., Abbaspour et al., 2009;Dastoorani and Poormohammadi, 2012;Kazemi-Rad and Mohammadi, 2015;. They also indicated a non-uniform pattern of temporal projections of global warming into air temperature, which is in agreement with results in our research. An increase in T min in autumn and winter may be beneficial since it may provide favorable conditions for early cultivation and, consequently, may reduce the duration of the cropping cycle since the crop's thermal energy requirement may be supplied during a shorter period . Such projections may further lead to lower crop water requirements if the positive effect of the shortened cropping cycle goes beyond the other probable negative effects of global warming on crop growth.

Climatic variables
During the 2041-2070 period, monthly P in the study area will always increase with the lowest increase occurring over the November-January period, and the highest one over the March-May period (Fig. 7). The RCP8.5 scenario led to the highest increases in monthly precipitation.  carried out a high-resolution assessment for entire Iran and reported that global warming might cause the highest increase in P along the coast of the Caspian sea, where the study area is located. They believed that such an increase might be attributed to an increase in air temperature, which consequently causes a considerable increase in the atmospheric water-holding capacity. An increase in P during wet seasons, which is in agreement with the findings of other researchers (e.g., Agarwal    installing costly drainage systems to mitigate further problems. An increase in P during seasons may also provide favorable conditions for the growth of weeds and pests and may enhance soil erosion, along with a considerable change in the soil available water (Enete and Amusa, 2010). Such consequences may reduce economic benefits by restraining the proper crop growth. Fig. 7 shows that monthly ET 0 may slightly increase by 2.4-10.3 % throughout October-March. The highest increase can be observed in autumn during the October-December period. However, monthly ET 0 is likely to be reduced by 0.6 % and 2.3 % in April and May, respectively. Such a decrease can be a result of a small increase in monthly average T max in the same period. The highest increase in ET 0 can be observed in autumn. Having a key role in the hydrological cycle, an increase in ET 0 may threaten agricultural sustainability since it may lead to a significant increase in the agricultural demand for water beyond its sustainable availability in a region. Such an increase may cause a significant reduction in water availability due to the overexploitation of surface and groundwater resources (Terink et al., 2013). Other researchers also support a projected increase in ET 0 in Iran (e.g., Terink et al., 2013;Karandish and Mousavi, 2018;Darzi-Naftchali and Karandish, 2019).

Drainage flux
Daily HYDRUS-simulated drainage fluxes for different drainage systems under future global warming projections (the ensemble average of 20 GCMs for each RCP scenario) are displayed in Fig. 8. While daily drain fluxes due to climate change projections both decreased and increased in different periods compared to the base scenario, they are more likely to increase, especially during the October-January period. Drainage fluxes are highly dependent on climatic variables, and particularly on P and ET 0 . While an increase in ET 0 may lead to a likely reduction in drainage fluxes due to an increase in crop water uptake, this positive effect may be reversed by a negative impact of increased precipitation on q.
The results in Fig. 8 indicate that the monthly green water surplus (GWS), defined as GWS=P-ET 0 for a particular month (Karandish and Mousavi, 2018), will significantly increase due to global warming and a corresponding increase in monthly precipitation, which will consequently enhance drainage fluxes. Simulated drainage fluxes for different drainage systems seem to be more uncertain early in the simulation rather than later.
Based on the ensemble average of 20 GCMs for each RCP scenario, average drainage fluxes during the growing season were calculated for each drainage system, and the results are summarized in Table 4. Compared to the base period, the average drainage fluxes over the cropping cycle may increase by 112.3 % (80.6-148.3 %), 66.9 % (62.0-74.9 %), 87.5 % (61.8-110.2 %), and 25.1 % (12.0 %-42.1 %) for the Bilevel, D 0.65 L 15 , D 0.65 L 30 , and D0.90L30 drainage systems, respectively. For all drainage systems, the RCP6.5 scenario leads to the highest increase in the average drainage flux over the cropping cycle, while the RCP2.6 scenario leads to the lowest increase.

Nitrate flux
Daily HYDRUS-simulated NO 3 fluxes for different drainage systems under future global warming projections (the ensemble average of 20 GCMs for each RCP scenario) are displayed in Fig. 9. The visual inspection shows a dramatic increase in NO 3 fluxes during the October-January period when monthly precipitations are also relatively high. This pattern is also a consequence of increased drainage fluxes in these future climate scenarios (Fig. 9). Fig. 9 also shows that during the October-January period, the range of uncertainties arising from the projections of different RCPs is also higher than during the other periods. While the Bilevel and D 0.65 L 30 drainage systems registered increased NO 3 fluxes during the entire growing cycle, the D 0.65 L 15 and D 0.90 L 30 drainage systems experienced a considerable decrease in NO 3 fluxes during the February-March period. More similar predicted NO 3 fluxes were obtained for all drainage systems during the February-March period, indicating less uncertainty during this time. Such results may be attributed to lower NO 3 soil concentrations during this period due to high crop N uptake during previous months (Darzi-Naftchali et al., 2017;DeDatta, 1981).
Based on the ensemble average of 20 GCMs for each RCP scenario, the minimum, maximum and average NO 3 fluxes during the growing season were calculated for each drainage system, and results are summarized in Table 4. The average NO 3 flux over the cropping cycle may increase by 53. 3-74.8 %, 18.6-23.9 %, 71.3-124.9 %, and 4-32.6 % in the Bilevel, D 0.65 L 15 , D 0.65 L 30 , and D 0.90 L 30 drainage systems, respectively. For all drainage systems, the RCP4.5 and RCP6.5 scenarios may result in the highest increase in the average NO 3 flux over the cropping cycle.

Seasonal nitrate losses
The seasonal NO 3 losses for different drainage systems were computed based on the ensemble average of 20 GCM projections for each RCP scenario (Table 4). Regardless of the time (i.e., either during the base or future periods), the lowest NO 3 loss always occurred for the D 0.65 L 15 drainage system, where drains are installed at a shallower depth and at a smaller distance. On the other hand, the D 0.90 L 30 drainage system always had the highest NO 3 losses during the growing seasons, indicating that NO 3 losses increase when drains are installed at deeper depths and larger distances (Cooke et al., 2002;Burchell, 2003;Yoon et al., 2006). These results indicate the importance of selecting a proper drainage system when sustainable agriculture and a safe environment are the primary concern.
Table 4 shows that climate change may result in higher seasonal NO 3 losses for all drainage systems. Compared to the base period, these increases will be 66.8 % (53.3-74.8 %), 20.5 % (18.6-23.9 %), 96.8 % (71.3-124.9 %), and 18.2 % (4-32.6 %) for the Bilevel, D 0.65 L 15 , D 0.65 L 30 , and D 0.90 L 30 drainage systems, respectively. The lowest and the highest increases in NO 3 losses are expected to occur in the D 0.90 L 30 and D 0.65 L 30 drainage systems, respectively.
Increased NO 3 losses due to climate change may be attributed either to increased drain fluxes or the field management practices and, in particular, to the N-fertilization management. Drain fluxes can be controlled by various factors, including climatic variables, especially P and ET 0 , soil hydraulic properties, the design of the drainage system, and crop management. Hence, increased seasonal NO 3 losses may be partially attributed to increased drainage fluxes resulting from a profound increase in surplus precipitation (i.e., P-ET 0 ) due to global warming (Table 4). While future changes in climate variables may be inevitable, negative consequences of climate change projections on drain fluxes may be alleviated by improving the soil water holding capacity. Such improvements may be achieved by applying various organic amendments, such as manure or gypsum (Udayasoorian et al., 2009), organic mulches (Chukalla et al., 2015), or proper tillage practices (Liu et al., 2013). Organic mulching is also reported to have a positive effect on reducing nutrient leaching in agricultural lands (De Vita et al., 2007;Mitchell et al., 2012), and is commonly proposed as an effective measure to prevent water pollution under intensive agriculture (Azooz and Arshad, 1996;Triplett and Dick, 2008;Constantin et al., 2010;Chukalla et al., 2017). In particular, organic mulching may restrict N-leaching from the soil surface layer, where sufficiently high  concentrations of N enable continued uptake by roots (Phogat et al., 2013). Drain depths and spacings also affect drain fluxes (Wahba and Christen, 2006). Jafari-Talukolaee et al. (2015) reported that increasing drain spacing might reduce drainage volumes in the subsurface-drained paddy field. The impact of installing either deep or shallow drainage systems on drainage water quality/quantity in Southern Australian irrigated lands was investigated by Christen and Skehan (2001). Their results demonstrated that lower drainage volumes occur from drainage systems with reduced drain depths and spacings. The long-term field investigation in the paddy field in northern Iran indicated that shallow drains were more effective than deeper ones in controlling the water table depth during the first three years after the installation of subsurface drainage while the reverse trend was observed in the fourth year (Jafari-Talukolaee et al., 2016). However, shallow drains may produce serious environmental problems since they enhance leaching and, consequently, may increase seasonal NO 3 losses (Darzi-Naftchali et al., 2017). Hence, drain depths and spacing need to be optimized depending on the main agricultural or environmental concerns and future global warming projections.
N-fertilization management is among the most crucial field practices, which significantly affect NO 3 losses. Numerous researchers demonstrated the close relationship between the N-fertilizer application rate/timing and N-leaching in agricultural lands (Mosier et al., 2002;Gasser et al., 2002;Haverkort et al., 2003;Daudén and Quilez, 2004;Alva et al., 2006;Barton and Colmer, 2006;Hutton et al., 2008;Wei et al., 2009;Jia et al., 2014;Karandish and Šimůnek, 2017). Besides, it is worth noting that the type of N-fertilizer can also affect the Nleaching rate (Jia et al., 2014). On the other hand, N is an essential crop nutrient, that profoundly affects crop growth and yield (Wienhold et al., 1995;Jia et al., 2014). An insufficient N supply in soil may result in severe yield and economic losses (Wienhold et al., 1995;Jia et a1., 2014). Therefore, further research should be carried out to find out the environmentally/economically safer applications of N-fertilizers under future global warming predictions.

Conclusion
Field experiments and modeling analyses involving a subsurfacedrained field of rained-canola, as a winter crop, were carried out to evaluate the integrated influence of subsurface drainage and global warming on probable future water and N losses. Our quantitative assessment successfully evaluated the capability of the HYDRUS model to predict soil water contents, soil nitrate concentrations, and drainage/ nitrate fluxes for various drainage systems. While the downscaled weather data from 20 GCMs and four RCP scenarios projected both decreases and increases in different future climatic variables, the ensemble averages predicted an increase in air temperatures, precipitation P, and potential evapotranspiration ET 0 for the 2041-2070 period. Surplus precipitation (defines as P-ET 0 ) is also likely to increase under climate change, which consequently resulted in a considerable increase in average values of drainage/nitrate fluxes during the growing season. Such increases may lead to an increase of 4-125 % in seasonal N losses for various drainage systems in the coming future. Our results demonstrate that negative environmental consequences of global warming, in terms of its effect on water and nitrate losses, may be alleviated when drains are installed at a shallow depth of 65 cm and at a spacing of 15 m (D 0.65 L 15 ). Such a system should be more beneficial regarding either the water table control or reducing seasonal N losses from rainfed cropping systems. Based on our results, it can be concluded that achieving sustainable rainfed agriculture under global warming requires further serious attempts to identify the best field management practices aiming at diminishing environmental consequences.
Additionally, the HYDRUS model could be a proper alternative to the costly and labor/time-consuming field investigations to determine the optimal management scenarios in rainfed lands under global warming conditions. While simulating the effects of climate change projections on crop yield was not among the objectives of this research, such effects may be of high importance when economic interests are considered. Further research is thus needed to compare both economic and environmental consequences to obtain a full picture of the coming future.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.