Impacts of Climate Change and Different Crop Rotation Scenarios on Groundwater Nitrate Concentrations in a Sandy Aquifer

Nitrate in groundwater is a major concern in agricultural sub-watersheds. This study assessed the impacts of future climate and agricultural land use changes on groundwater nitrate concentrations in an agricultural sub-watershed (Norfolk site) in southern Ontario, Canada. A fully integrated hydrologic model (HydroGeoSphere) was used in combination with the root zone water quality model (RZWQM2) (shallow zone) to develop water flow and nitrate transport models. Three climate change models and three crop rotations (corn-soybean rotation, continuous corn, corn-soybean-winter wheat-red clover rotation) were used to evaluate the potential impact on groundwater quality (nine predictive scenarios). The selected climate change scenarios yielded less water availability in the future period than in the reference period (past conditions). The simulated nitrate nitrogen (Nitrate-N) concentrations were lower during the future period than the reference period. The continuous corn land use scenario produced higher Nitrate-N concentrations compared to the base case (corn-soybean rotation). However, the best management practices (BMP) scenario (corn-soybean-winter wheat-red clover rotation) produced significantly lower groundwater nitrate concentrations. BMPs, such as the one examined herein, should be adopted to reduce potential negative impacts of future climate change on groundwater quality, especially in vulnerable settings. These findings are important for water and land managers, to mitigate future impacts of nutrient transport on groundwater quality under a changing climate.


Introduction
Climate has always played an important role in the agriculture industry and has shaped the outlook of agriculture all over the world [1,2]. Farmers select their crops based on climatic factors such as annual rainfall, low and high temperatures, and growing season length [3,4]. In recent decades, farming activities in agricultural watersheds have become a major concern for water quantity and quality [5,6]. Conventional farming practices in some locations were to apply fertilizers uniformly at high rates without considering the crop and soil requirements [7][8][9] resulting in leaching of excess nutrients to groundwater especially after a heavy rainfall event [9]. However, more farmers are moving towards better nutrient management considering the negative impacts of high fertilizer rates on crop productivity and the environment [10][11][12]. Elevated concentrations of nitrogen and phosphorous in water bodies can cause many health and environmental problems, such as methemoglobinemia (blue baby syndrome) and eutrophication, respectively [13,14]. Therefore, it is necessary to evaluate the impact of land and crop management practices on groundwater quality under changing climate conditions. The Norfolk research site with a catchment area of 155 km 2 ( Figure 1) is located in Norfolk County in southwestern Ontario and it is part of the Lynn River sub-watershed. The average annual precipitation at for the 2010-2016 period was 1009 mm (ranging between 830 to 1349 mm) [48] . Both surface water and groundwater flows from northwest to southeast in this sub-watershed. The river outlet is located at the southeastern boundary of the modelled domain and the surface elevation varies from 188 to 265 meters above sea level (masl). The study site is dominated by cultivated crops (mainly corn (~16%-25%), soybeans (~16%-25%), and specialty crops (~15%)) except near the southwestern boundary, where forest covers large areas of the research site [49].  [50], and land use map for Norfolk site for 2012 [49]. Red symbols indicate the positions of the monitoring wells (MWs) used in this study.

Geological and Hydrogeological Framework
The study site is comprised of the Norfolk sand plain, having a thick layer of unconsolidated material above the Paleozoic bedrock. Borehole data from monitoring wells installed for a Tier 3 (water budget and local risk assessment) source water protection study [51] was used to create a conceptual model for the study site with detailed definition of geological and hydrogeological layers ( Figure 2).
The top layer, commonly known as Norfolk sand plain (Whittlesey/Warren phase), identified as a thick layer of sand laid during the end of the Huron stadial consists of the most recent lake sediments. Whittlesey/Warren phase sediments are made up of medium sand or fine sand to medium sand and is split into two sublayers based on material grading [52]. The sublayer comprised of fine sand sediments sits at the bottom of the coarse-to-medium sand sublayer [52]. All shallow monitoring wells are screened in the Norfolk coarse sand layer and have high hydraulic conductivity (K) ranging from 5.5×10 -4 m s -1 to 3.6 ×10 -4 m s -1 [51]. The Norfolk fine sand layer has relatively lower K that ranges from 4.2×10 -5 m s -1 to 3.6 ×10 -7 m s -1 [51].  [50], and land use map for Norfolk site for 2012 [49]. Red symbols indicate the positions of the monitoring wells (MWs) used in this study.

Geological and Hydrogeological Framework
The study site is comprised of the Norfolk sand plain, having a thick layer of unconsolidated material above the Paleozoic bedrock. Borehole data from monitoring wells installed for a Tier 3 (water budget and local risk assessment) source water protection study [51] was used to create a conceptual model for the study site with detailed definition of geological and hydrogeological layers ( Figure 2).

Numerical Modeling
A two-step numerical modeling approach was adopted for this study. A 1D vadose zone crop model RZWQM [26], was selected to model the detailed crop growth and Nitrate-N leaching though the root zone. The RZWQM2 model was calibrated and validated for various cropping systems using soil and leachate data from soil samples and suction lysimeters installed in the selected fields. The RZWQM2 model outputs such as daily Nitrate-N mass and daily leaf area index (LAI) were used as inputs in the integrated hydrological model. The field-scale RZWQM2 model provides comprehensive outputs in terms of plant growth and vadose zone hydrogeology that can be neglected in watershed scale models. The detailed vadose zone modeling procedure is described in Reference [57].
A fully integrated (surface-subsurface) HGS model was used to simulate groundwater flow and transport of Nitrate-N in the groundwater at the study site. HGS uses a two-dimensional (2D) wave The top layer, commonly known as Norfolk sand plain (Whittlesey/Warren phase), identified as a thick layer of sand laid during the end of the Huron stadial consists of the most recent lake sediments. Whittlesey/Warren phase sediments are made up of medium sand or fine sand to medium sand and is split into two sublayers based on material grading [52]. The sublayer comprised of fine sand sediments sits at the bottom of the coarse-to-medium sand sublayer [52]. All shallow monitoring wells are screened in the Norfolk coarse sand layer and have high hydraulic conductivity (K) ranging from 5.5 × 10 −4 m s −1 to 3.6 × 10 −4 m s −1 [51]. The Norfolk fine sand layer has relatively lower K that ranges from 4.2 × 10 −5 m s −1 to 3.6 ×10 −7 m s −1 [51].
The Norfolk fine sand layer is underlain by an Interstadial coarse sand layer (Ypsilanti layer) that represents a well-defined period of quick drops in lake levels resulting in relatively coarser sediments [52]. The Interstadial coarse sediment layer also has high K and deep wells are mostly screened in this layer with K ranging from 5.6 × 10 −6 to 1.70 × 10 −3 m s −1 [51]. The Interstadial fine sediment layer (Lake Maumee/Arkona) is dominated by a low energy deposition of very fine to fine sediments (very fine sand to clay) [52]. The thickness of the Interstadial fine sediment layer increases towards the east within the study area and it is present throughout all cross-sections.
The Port Stanley Till layer represents re-advance of the Erie lobe during the Port Burce stadial picked up the older Erie phase glaciolacustrine deposits and consists a clayey silt to sandy silt diamicton with a small percentage of coarse sediments [52]. The Erie Phase layer consists of glaciolacustrine sediments comprised of bedded sand, and silt and clay layers represent a period of ice retreat referred to as the Erie interstadial [53]. The Catfish Creek Till layer is the smallest and oldest sediment layer delineated that comprises a gravelly sandy till which correlates to the Nissouri stadial [52]. The Paleozoic bedrock is very deep (> 30 m below ground surface) and acts as an aquitard [54].
Tritium and other stable isotopes analyses on groundwater samples collected during field work in July 2014 indicated that modern groundwater is present at the research site [55]. The isotope results showed that the recharge water is only 5 to 10 years old at the research site and the position of the data relative to the global meteoric water line indicates cold climate region groundwater [55,56]. The nitrogen-15 isotope results indicated that agricultural activities are the main source of groundwater Nitrate-N concentrations at the research site, with only one well location where a septic system was identified as the main cause of groundwater Nitrate-N concentrations [55]. High infiltration rates and the sandy nature of the soil indicates that most of the contribution from the subwatershed to the Lynn river is through baseflow. The regional groundwater flows from northwest to southeast towards the Lake Erie.
In summary, the developed conceptual hydrogeological model indicates that the groundwater flow and Nitrate-N transport occur in unconsolidated shallow sand layers with high K. The groundwater present at the research site is relatively fresh and nitrate transported to the aquifer by infiltration will first reach the shallow groundwater and will flow southwards towards Lake Erie. The water quality in the aquifer can be controlled by water infiltration (recharge) and nitrate leaching from the catchment area Cross-sections were made using all Tier-3 monitoring wells and elevation points were inferred from these geologic cross sections. The geological layers were interpolated at a resolution of 1 × 1 m using an inverse distance weighted scheme and the bottom of the conceptual model was set at the top of bedrock (175 masl). A total of nine hydrogeological layers were selected for numerical modeling purposes ( Figure 2).

Numerical Modeling
A two-step numerical modeling approach was adopted for this study. A 1D vadose zone crop model RZWQM [26], was selected to model the detailed crop growth and Nitrate-N leaching though the root zone. The RZWQM2 model was calibrated and validated for various cropping systems using soil and leachate data from soil samples and suction lysimeters installed in the selected fields. The RZWQM2 model outputs such as daily Nitrate-N mass and daily leaf area index (LAI) were used as inputs in the integrated hydrological model. The field-scale RZWQM2 model provides comprehensive outputs in terms of plant growth and vadose zone hydrogeology that can be neglected in watershed scale models. The detailed vadose zone modeling procedure is described in Reference [57].
A fully integrated (surface-subsurface) HGS model was used to simulate groundwater flow and transport of Nitrate-N in the groundwater at the study site. HGS uses a two-dimensional (2D) wave approximation of the Saint Venant equation [58] for surface flow. The finite-element approach was selected for groundwater modeling and Algomesh [59] was used to create a 2D mesh of the study site. The Saint Venant equation to quantify surface flow requires Manning coefficients that are used to calculate friction slopes, and a specified rill storage height to determine how much water is stored in the surface depressions. Rill storage is an important factor in rural settings and influences the surface Sustainability 2020, 12, 1153 6 of 25 runoff generated from the study site. Another important parameter for the surface flow simulation is the coupling length between the surface and subsurface domains [60]. All the overland flow properties were assumed to be uniform throughout the study area. The model uses the Richards' equation [61] to calculate three-dimensional (3D) subsurface flow in variably saturated media. Evapotranspiration (ET) affecting both the surface and subsurface domains is modelled as a combination of plant transpiration and evaporation following the approach described in Reference [62]. The rate of transpiration from the modeling domain depends on LAI, root distribution function, and variability in available soil moisture content. Transpiration nonlinearly depends on soil moisture content and wilting point (θ wp ), field capacity (θ fc ), oxic limit (θ o ), and anoxic limit (θ an ). Maximum rooting depths and root extraction functions are also important factors in simulating evapotranspiration.

Data Collection
A hydrologically enforced digital elevation model (DEM) was provided by the Grand River Conservation Authority (GRCA) with a resolution of 1 × 1 m. Land use data for the study site was downloaded from the Agriculture and Agri-food Canada (AAFC) website [49]. At the site, groundwater elevation and Nitrate-N data was collected from selected monitoring wells on a bi-monthly basis, from June 2014 to July 2016 [55]. Outflow data was collected for 1990-2016 from a Water Survey of Canada gauge near the outlet of the watershed [63]. The climate data, including precipitation and temperature, were collected from the local weather station [48]. The solar radiation data for both sites was downloaded from NASA at a spatial resolution of 0.5 • × 0.5 • [64]. Potential evapotranspiration (PET) zones were delineated based on different land uses. Details of data types, spatial and temporal resolutions, and data sources are described in Table 1. The 3D HGS model for the study site consists of 222,363 nodes and 415,560 elements. A critical depth boundary condition was applied in HGS at the watershed outlet with all other boundaries specified as no flow. The horizontal and vertical hydraulic conductivities (K xy and K z , respectively) along with surface flow and evapotranspiration parameters for various cropping systems were calibrated during the "steady state" calibration. Some of the parameter values were adopted from the literature values.
For initial conditions, a spin up period (1000 years) was used to obtain reasonable groundwater elevation and flow values and avoid inappropriate initial conditions. The initial conditions were based on average rainfall and PET (1990-2009) and average values for Nitrate-N concentrations from different land use types generated in RZWQM2. Once the initial heads and concentrations were obtained, the model was calibrated for steady state conditions. The PEST software package [65] was used to calibrate the K values and average outlet flow for the study area. Transpiration and evaporation limiting saturations and canopy storage parameters were also calibrated to match the evapotranspiration generated by RZWQM2. Results showed a good agreement between modelled Sustainability 2020, 12, 1153 7 of 25 and observed groundwater elevations (masl) with an R 2 value of 0.96 and a root mean square error (RMSE) of 2.7 m (Figure 3a). Similarly, observed and modelled Nitrate-N concentrations in sampled monitoring wells also had a good agreement with an R 2 value of 0.91 and a RMSE of 1.83 mg L −1 (Figure 3b). The final calibrated values of all input parameters are provided in Table 2, Table 3, and  Table 4.
based on average rainfall and PET (1990PET ( -2009 and average values for Nitrate-N concentrations from different land use types generated in RZWQM2. Once the initial heads and concentrations were obtained, the model was calibrated for steady state conditions. The PEST software package [65]was used to calibrate the K values and average outlet flow for the study area. Transpiration and evaporation limiting saturations and canopy storage parameters were also calibrated to match the evapotranspiration generated by RZWQM2. Results showed a good agreement between modelled and observed groundwater elevations (masl) with an R 2 value of 0.96 and a root mean square error (RMSE) of 2.7 m (Figure 3a). Similarly, observed and modelled Nitrate-N concentrations in sampled monitoring wells also had a good agreement with an R 2 value of 0.91 and a RMSE of 1.83 mg L -1 (Figure 3b). The final calibrated values of all input parameters are provided in Table 2, Table 3, and  Table 4.     The steady state calibration was followed by transient state calibration and validation with a special focus to adjust the transport parameters, such as dispersivities and the specific storage to accurately model the Nitrate-N in the groundwater. The outflow was calibrated for 2013-2014 and validated for 2015-2016, while the groundwater elevation and groundwater Nitrate-N concentrations were spatially calibrated and validated. Automatic calibration software PEST-HP [65] was used for transient state parameters [57]. The detailed transient state calibration and validation results can be accessed in Reference [57]. For this reason, only a brief description of the data calibration and validation procedure is included in this paper. Calibrated values of the transient state parameters are provided in Table 5. The calibrated model was used for the evaluation of crop rotation scenarios and climate change projections in this study.

Sensitivity Analysis
Relative Sensitivity (S r ) of averaged groundwater elevation, outflow, and groundwater Nitrate-N concentrations was evaluated using a procedure adopted from Reference [75] under steady state conditions. The calibration parameters subjected to sensitivity analysis were K, precipitation, and PET ( Table 6). The K values selected for sensitivity analysis were based on maximum and minimum field calculated for Layers 1-3 and Layers 4-5 (Table 2). For climate parameters, maximum and minimum values between 2010 to 2016 were selected [48]. The absolute value of S r may vary from 0 to infinity with smaller values representing a more robust model (minor changes in output by changing the parameter values) while larger values represent greater parameter sensitivity. Table 6. Relative Sensitivity (S r ) of groundwater elevation, outflow, and groundwater Nitrate-N concentrations to changes in K, precipitation, and PET for steady state conditions.

Parameter
Groundwater Elevation S r Outflow S r Nitrate-N S r Sensitivity analysis results showed that both outflow and groundwater Nitrate-N concentrations were more sensitive to the modified parameters than groundwater elevation. All three model outputs were more sensitive to the changes in precipitation and PET compared to changes in K. Only groundwater Nitrate-N was sensitive to changes in K values for Layers 1-3 (Table 6). These results indicate that climate change and crop coverage can highly influence the model outputs.

Future Climate Projections
Dynamically downscaled regional climate model (RegCM) climate change projections (25 km resolution) under the representative concentration pathway (RCP) 8.5 emission scenario and driven by five global climate models (GCMs) (CanESM2, HadGEM2-ES, GFDL-ESM2M, IPSL-CM5A-LA, and MPI-ESM-MR) were downloaded from the Ontario Climate Change Data Portal (CCDP) [76]. This data corresponds to a reference period spanning from 1986 to 2005 with focus on a mid-century future period (2040-2059). Regional dynamic downscaling provides better representation of extreme and mean climate conditions than statistical downscaling [77]. The selection of climate change projections was based on future temperature and precipitation changes relative to the reference period. The criteria were to select a projection with maximum precipitation change and temperature change, minimum precipitation change and maximum temperature change, and business as usual (both temperature and precipitation change close to observed data). To assess this, the annual means of climate change projections from all five driving GCMs for the reference period were compared ( Figure 4). The MPI (MPI-ESM-MR) projection was selected for the highest precipitation changes, CanESM (CanESM2) was selected for the highest temperature and lowest precipitation change, and GFDL (GFDL-ESM2M) was initially selected as a business as usual scenario based on annual temperature and precipitation changes for the reference period.

Future Climate Projections
Dynamically downscaled regional climate model (RegCM )climate change projections (25 km resolution) under the representative concentration pathway (RCP) 8.5 emission scenario and driven by five global climate models (GCMs) (CanESM2, HadGEM2-ES, GFDL-ESM2M, IPSL-CM5A-LA, and MPI-ESM-MR) were downloaded from the Ontario Climate Change Data Portal (CCDP) [76]. This data corresponds to a reference period spanning from 1986 to 2005 with focus on a mid-century future period (2040-2059). Regional dynamic downscaling provides better representation of extreme and mean climate conditions than statistical downscaling [77]. The selection of climate change projections was based on future temperature and precipitation changes relative to the reference period. The criteria were to select a projection with maximum precipitation change and temperature change, minimum precipitation change and maximum temperature change, and business as usual (both temperature and precipitation change close to observed data). To assess this, the annual means of climate change projections from all five driving GCMs for the reference period were compared ( Figure 4). The MPI (MPI-ESM-MR) projection was selected for the highest precipitation changes, CanESM (CanESM2) was selected for the highest temperature and lowest precipitation change, and GFDL (GFDL-ESM2M) was initially selected as a business as usual scenario based on annual temperature and precipitation changes for the reference period. However, monthly mean temperature ( Figure 5a) and precipitation (Figure 5b) data showed that HadGEM (HadGEM2-ES) better represents the observed condition for the reference period and it was selected as a replacement for the GFDL scenario. The final three climate change scenarios used for this study were CanESM, HadGEM, and MPI. For the study site, local climate data was downloaded from Environment and Climate Change Canada for the reference period (1986-2005) [48]. The average observed annual and monthly precipitation for the reference period was 976.71 mm and 80.61 mm, respectively [48]. The observed temperature during the reference period ranged between -24 °C and 29 °C [48], while the climate models' minimum and maximum temperatures for the reference period were:  However, monthly mean temperature ( Figure 5a) and precipitation (Figure 5b) data showed that HadGEM (HadGEM2-ES) better represents the observed condition for the reference period and it was selected as a replacement for the GFDL scenario. The final three climate change scenarios used for this study were CanESM, HadGEM, and MPI. For the study site, local climate data was downloaded from Environment and Climate Change Canada for the reference period (1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005) [48]. The average observed annual and monthly precipitation for the reference period was 976.71 mm and 80.61 mm, respectively [48]. The observed temperature during the reference period ranged between −24 • C and 29 • C [48], while the climate models' minimum and maximum temperatures for the reference period were: CanESM: −15.7 • C and 29 • C; HadGEM: −23.7 • C and 30.3 • C; and MPI: −24 • C and Sustainability 2020, 12, 1153 12 of 25 29 • C. Several studies have suggested the application of a bias correction technique before using the RCM data in order to better represent the local subwatershed scale climatic conditions [78][79][80][81]. There are different methods to bias correct the climate data which include the delta change approach, local intensity scaling, daily bias correction, and quantile mapping [82]. For this study, a quantile mapping technique was used to remove the biases between simulated and observed temperature and precipitation variables. Details of the quantile mapping technique are described in [82]. The net recharge (Pnet), that is precipitation-PET, was calculated on a monthly basis to estimate the variation in water available for potential recharge and surface runoff throughout the year. PET was calculated using the Penman-Monteith equation [83]. Once the climate models were corrected for their bias, outputs for the average monthly temperature and precipitation for the reference period were reproduced.
Sustainability 2020, 12, x FOR PEER REVIEW 13 of 25 using the RCM data in order to better represent the local subwatershed scale climatic conditions [78][79][80][81]. There are different methods to bias correct the climate data which include the delta change approach, local intensity scaling, daily bias correction, and quantile mapping [82]. For this study, a quantile mapping technique was used to remove the biases between simulated and observed temperature and precipitation variables. Details of the quantile mapping technique are described in [82]. The net recharge (Pnet), that is precipitation-PET, was calculated on a monthly basis to estimate the variation in water available for potential recharge and surface runoff throughout the year. PET was calculated using the Penman-Monteith equation [83]. Once the climate models were corrected for their bias, outputs for the average monthly temperature and precipitation for the reference period were reproduced.  [76], model ensemble [76], and observed data at a local climate station (Delhi, Ontario).  [76], model ensemble [76], and observed data at a local climate station (Delhi, Ontario).

Crop Rotation Scenarios
Many different agricultural land use change scenarios can be devised, including those based on crop area coverage, forestation/deforestation, urban/rural area coverage, etc. In the current study, the focus was to define and quantify the transport of excess nutrients, specifically nitrogen related to cash crop modifications and variable weather, into groundwater to anticipate and mitigate potential water quality impacts. Three future crop scenarios were selected to consider potential variation and extremes in nitrogen application: (1) crop area under corn-soybean rotation will continue to be under the same corn-soybean rotation in the future (business as usual scenario or base case scenario), (2) crop area under corn-soybean rotation changes to continuous corn plantation (cash crops only) in the future (worst case scenario), and (3) crop area under corn-soybean rotation changes to corn-soybean-winter wheat-red clover rotation (BMP scenario).

Statistical Evaluation
The surface flow at the watershed outlets as well as groundwater elevation and Nitrate-N concentrations in monitoring wells were compared based on four seasons: winter (Dec, Jan, Feb), spring (Mar, Apr, May), summer (June, July, Aug), and fall (Sep, Oct, Nov). This quarterly comparison provides useful insight on how future climate change will vary seasonally and how it will impact Nitrate-N concentrations during the growing season relative to the rest of the year. Box plots were created to compare the data distribution based on percentile for groundwater elevation, surface flow, and Nitrate-N concentrations for each season using Grapher 11 (Golden Software, LLC, Golden, CO, USA). The ensemble means of surface water, groundwater elevations, and nitrate concentrations were compared using one-way analysis of variance (ANOVA) followed by means comparison using a least significant difference (LSD) for significantly different treatments (p < 0.05) in SPSS (IBM SPSS®, New York, NJ, USA).

Climate Projections
The ensemble means of future (2040-2059) simulations predicted an increase of 0.37 • C (May) to 0.19 • C (February) for monthly temperatures (Figure 6a). The envelope of uncertainty remains relatively equal throughout the year except for March, in which the MPI scenario showed a decrease of 1.3 • C in mean temperature between future and historic model data. The ensemble means of future simulations predicted an increase in monthly precipitation for every month from December to April and a decrease in monthly precipitation from the months of May to November (Figure 6b). All three models suggest high precipitation during winter months. The CanESM model showed an erratic behavior for precipitation change between future and reference periods during the summer months. The envelope of uncertainty remains high (6% to 59%) and the ensemble means predicted less rainfall in the future during the growing season (May to October) (Figure 6b). The ensemble means of calculated Pnet values had a maximum increase of 27.4 mm in April and a maximum decrease of 19.8 mm in July. The envelope of uncertainty showed a range between 1 mm (April) and 45 mm (August). Pnet variations were high between different models during the months of July to September (Figure 6c). The envelope of uncertainty for temperature change also had similar trends for these months (Figure 6a). Overall, the envelope of uncertainty for Pnet illustrated that the growing season (May to October) has less available water for potential surface runoff and groundwater recharge; however, the actual rate of availability depends on which model represents the future climate conditions more accurately (Figure 6c). Overall, the envelope of uncertainty indicated that future growing seasons will have less available water for potential surface runoff and groundwater recharge compared to the past. However, the rate of availability depends on which model will represent the future climate conditions more accurately at the study site (Figure 6c).

Outflow
Outflow at the study site was higher for the reference period compared to the future period for all three climate change scenarios and flow at the outlet ranged from 0.75 to 2.1 m 3 s −1 (Figure 7). The monthly Pnet (Figure 6c) values at the study site indicated that the amount of water available for hydrologic processes throughout the winter and spring months was close to zero for the reference and future periods, while the amount of water available for hydrologic processes for summer and fall months were less for the future period compared to the reference period. Higher flows at the outlet were more prominent during summer and fall periods for the study site (Figure 7).

Outflow
Outflow at the study site was higher for the reference period compared to the future period for all three climate change scenarios and flow at the outlet ranged from 0.75 to 2.1 m 3 s -1 (Figure 7). The monthly Pnet (Figure 6c) values at the study site indicated that the amount of water available for hydrologic processes throughout the winter and spring months was close to zero for the reference and future periods, while the amount of water available for hydrologic processes for summer and fall months were less for the future period compared to the reference period. Higher flows at the outlet were more prominent during summer and fall periods for the study site (Figure 7). The average values for model ensemble outflows at the study site for the reference and future periods were 1.12 and 1.09 m 3 s -1 , respectively ( Table 7). The results for the model ensemble indicated that the outflow will be decreased (statically significant) in the future period by 2.86% at the study site. Table 7. Comparison of outflow changes between the reference and future model ensembles at the study site.

Period
Outflow (m 3 s -1 ) % Change The average values for model ensemble outflows at the study site for the reference and future periods were 1.12 and 1.09 m 3 s −1 , respectively ( Table 7). The results for the model ensemble indicated that the outflow will be decreased (statically significant) in the future period by 2.86% at the study site.  (a versus b) are significantly different at a significance level of 0.05.

Groundwater Elevation
Modeling results of the groundwater elevations show that future groundwater elevations are generally lower, especially in the spring/summer periods (March to August) for the HadGEM and MPI climate change scenarios (Figure 8). However, the CanESM model increased water levels at the study site for future periods compared to the reference period at both monitoring wells (Figure 8). The groundwater elevation at two monitoring wells are discussed and presented in detail. The groundwater elevations for all climate change scenarios ranged between 220 masl and 223.  Figure 8b). This might be due to high rainfall for the MPI model during the reference period compared to all other climate change models at the study site (Figure 6b). The MPI scenario for the future period showed the greatest decrease in groundwater elevation (based on the median) for both monitoring wells (Figure 8).

Groundwater Elevation
Modeling results of the groundwater elevations show that future groundwater elevations are generally lower, especially in the spring/summer periods (March to August) for the HadGEM and MPI climate change scenarios (Figure 8). However, the CanESM model increased water levels at the study site for future periods compared to the reference period at both monitoring wells (Figure 8). The groundwater elevation at two monitoring wells are discussed and presented in detail. The groundwater elevations for all climate change scenarios ranged between 220 masl and 223.  Figure 8b). This might be due to high rainfall for the MPI model during the reference period compared to all other climate change models at the study site (Figure 6b). The MPI scenario for the future period showed the greatest decrease in groundwater elevation (based on the median) for both monitoring wells (Figure 8).

Groundwater Nitrate-N Concentrations
For the climate change scenarios and land use scenarios, the Nitrate-N concentration is generally lower for the future conditions at the study site ( Figure 9). The Nitrate-N concentrations at MW-04 ranged between 8 to 26 mg L −1 for the corn-soybean land use scenario (Figure 9a), between 9 to 31 mg L −1 for the continuous corn land use scenario (Figure 9b), and between 5 to 29 mg L −1 for the BMP scenario (Figure 9c). At the MW-07 monitoring well, the Nitrate-N concentrations ranged between 3 to 13 mg L −1 for the corn-soybean land use scenario (Figure 10a), between 5 to 21 mg L −1 for the continuous corn land use scenario (Figure 10b), and between 1 to 13 mg L −1 for the BMP scenario (Figure 10c). The MPI scenario had the highest values for Nitrate-N concentrations for all three land use scenarios (Figure 10). Sustainability 2020, 12, x FOR PEER REVIEW 18 of 25 L -1 for the continuous corn land use scenario (Figure 9b), and between 5 to 29 mg L -1 for the BMP scenario (Figure 9c). At the MW-07 monitoring well, the Nitrate-N concentrations ranged between 3 to 13 mg L -1 for the corn-soybean land use scenario (Figure 10a), between 5 to 21 mg L -1 for the continuous corn land use scenario (Figure 10b), and between 1 to 13 mg L -1 for the BMP scenario ( Figure 10c). The MPI scenario had the highest values for Nitrate-N concentrations for all three land use scenarios ( Figure 10). The model ensemble results at both MW-04 and MW-07 suggested a small increase of 0.02% and 0.16%, respectively, in groundwater elevation in the future period ( Table 8). The groundwater elevation results were in contrast with outflow results when the reference period was compared with the future period, suggesting that groundwater elevations are less effected by future climate changes compared to outflows in the study sub-watershed (Tables 7 to 9).  The Nitrate-N concentrations were also decreased in the future model ensemble for all three land use and climate scenarios. The smallest reduction (0.14%) in Nitrate-N concentration was reported in MW-04 for the continuous corn scenario and the highest reduction (-47.36%) was reported in MW-07 for the BMP scenario (Tables 8 and 9). The results of the comparison of different land use scenarios at the study site with corn-soybean as the base case scenario suggested that continuous corn production caused an increase of 47.56% and 79.43% during the reference period for MW-04 and MW07 respectively, while continuous corn production increased Nitrate-N concentrations by 51.21% and 89.52% in MW-04 and MW-07 respectively, for the future period (Tables 8 and 9).
The likelihood of higher Nitrate-N concentrations in the groundwater in the future compared to the reference period is greater if continuous corn is planted (Tables 8 and 9). The comparison of the base case scenario with the BMP scenario produced a reduction of 18.79% and 32.29% during the reference period for MW-04 and MW07, respectively (Tables 8 and 9). Further, BMP adoption at the study site reduced Nitrate-N concentrations by 31.49% and 47.36% in MW-04 and MW-7 respectively, compared to the base case scenario for the future period (Tables 8 and 9). The comparison of future and reference periods showed that adopting the BMP scenario (a rotation of corn-soybean-winter wheat-red clover) can have more positive impacts on groundwater quality in the future compared to the reference period (Tables 8 and 9). The Nitrate-N concentrations are almost two times higher for The model ensemble results at both MW-04 and MW-07 suggested a small increase of 0.02% and 0.16%, respectively, in groundwater elevation in the future period ( Table 8). The groundwater elevation results were in contrast with outflow results when the reference period was compared with the future period, suggesting that groundwater elevations are less effected by future climate changes compared to outflows in the study sub-watershed (Tables 7-9).
The Nitrate-N concentrations were also decreased in the future model ensemble for all three land use and climate scenarios. The smallest reduction (0.14%) in Nitrate-N concentration was reported in MW-04 for the continuous corn scenario and the highest reduction (-47.36%) was reported in MW-07 for the BMP scenario (Tables 8 and 9). The results of the comparison of different land use scenarios at the study site with corn-soybean as the base case scenario suggested that continuous corn production caused an increase of 47.56% and 79.43% during the reference period for MW-04 and MW07 respectively, while continuous corn production increased Nitrate-N concentrations by 51.21% and 89.52% in MW-04 and MW-07 respectively, for the future period (Tables 8 and 9).  The likelihood of higher Nitrate-N concentrations in the groundwater in the future compared to the reference period is greater if continuous corn is planted (Tables 8 and 9). The comparison of the base case scenario with the BMP scenario produced a reduction of 18.79% and 32.29% during the reference period for MW-04 and MW07, respectively (Tables 8 and 9). Further, BMP adoption at the study site reduced Nitrate-N concentrations by 31.49% and 47.36% in MW-04 and MW-7 respectively, compared to the base case scenario for the future period (Tables 8 and 9). The comparison of future and reference periods showed that adopting the BMP scenario (a rotation of corn-soybean-winter wheat-red clover) can have more positive impacts on groundwater quality in the future compared to the reference period (Tables 8 and 9). The Nitrate-N concentrations are almost two times higher for continuous corn compared to the BMP scenario. Also, more extreme climate changes (e.g., the MPI scenario) can increase the risk of groundwater contamination under continuous cash crop production at the study site.
The selected climate change scenarios yielded lower Pnet values in the future period, which means less water is available for hydrologic processes. The model ensemble results indicated that the lower Pnet reduced the surface water flows significantly at the study site for the future period. A previous study conducted in Prince Edward Island (PEI), Canada [27] also produced a decrease in groundwater elevation for future climate change scenarios in agricultural lands. Another modeling study conducted in the same major watershed as the current study (Grand River watershed) resulted in a reduction in annual stream flow during the future (2040-2059) period compared to the reference period   [84]. The low annual Pnet values were the main reasons for lower streamflow and lower groundwater elevation. The findings of Reference [44] contrasted with the current study. Their results suggested an increase in groundwater discharge (spring flows) during the future period in a fractured bedrock aquifer of Chateauguay River watershed in Quebec, Canada. However, the Pnet values of the selected future climate change scenarios were higher compared to the reference period [44]. The future hydrologic water balance is highly dependent on study location and associated future climate model predictions.
The reduction in future Nitrate-N concentrations is in contrast with previous studies investigating groundwater or surface water quality conducted on agricultural lands [27,37,85,86]. These studies were based on increased nitrogen-based fertilizer application in future periods compared to reference periods. In the current study, a primary assumption was that current fertilizer application rates were also applied in the future periods. Only the impact of various land use scenarios (and their associated nitrogen application rates for the specific crops), including a BMP scenario, was studied. Therefore, the primary reason for lower nitrate concentrations observed in the modelled monitoring wells is less recharge (Pnet) which may result in less nitrate leaching from the surface. The lower nitrate concentrations can also mean that the nitrate could still be present in large amounts in the soil profile and can be leached down if precipitation occurs in higher amounts in the future. More research is needed to test climate change and land use change scenarios relating to nitrate impacts on groundwater from agricultural lands. The focus of this research was the comparison of land use change scenarios with similar inputs for reference and future periods. More numerical modelling studies are required to incorporate the changes in future nitrogen demands into these land use changes scenarios and to predict the impact on groundwater quality.

Conclusions
The numerical modeling results to represent outlet flow, groundwater elevation, and Nitrate-N concentrations for climate change and land use scenarios, showed: (1) a significant reduction in surface water flow in the future period, (2) a slight increase in groundwater elevation in the future period compared to the reference period at the study site, (3) a significant reduction in Nitrate-N concentrations in the future period compared to the reference period for all three climate change and all three land use scenarios if current fertilizer application rates are maintained, (4) a significant increase in groundwater Nitrate-N concentration (up to 80%) observed in monitoring wells if all the land under corn-soybean production is converted to continuous corn production for both reference and future periods, and (5) a significant reduction in Nitrate-N concentrations if the BMP scenario (corn-soybean-winter wheat-red clover rotation) is adopted at the study site for all three climate change scenarios during reference and future periods. The overall results suggest that BMPs should be adopted for unconsolidated hydrogeological settings to reduce potential impacts of Nitrate-N leaching on groundwater quality from agricultural lands. Planting of continuous corn along with the most extreme climate change scenario (MPI) can increase the risk of groundwater contamination.
Water managers and land use managers can use the results of this study and future similar studies in other jurisdictions to inform water protection and nutrient management policies in rural contexts to protect water at the source.