|
Predicted Contribution of Snowmelt to Subsurface Drainage Discharge in Two Subsurface-Drained Fields in Southern Quebec and Ontario
Ziwei Li1, Zhiming Qi1,*, Chandra A. Madramootoo1, Tiequan Zhang2
Published in Journal of the ASABE 67(1): 99-114 (doi: 10.13031/ja.15532). Copyright 2024 American Society of Agricultural and Biological Engineers.
1Department of Bioresource Engineering, McGill University, Sainte-Anne-de-Bellevue, Quebec, Canada.
2Harrow Research and Development Centre, Agriculture and Agri Food Canada, Harrow, Ontario, Canada.
*Correspondence: zhiming.qi@mcgill.ca
Submitted for review on 12 January 2023 as manuscript number NRES 15532; approved for publication as a Research Article and as part of the “Advances in Drainage: Selected Works from the 11th International Drainage Symposium” Collection by Community Editor Dr. Mohamed Youssef of the Natural Resources & Environmental Systems Community of ASABE on 27 July 2023.
Highlights
Abstract. Late winter/early spring has been recognized as a critical period for subsurface drainage discharge and associated nutrient losses due to snowmelt and rainfall in croplands under cold, humid climates in North America and Europe. Although the detachment and transport processes of soil particles under snowmelt and rainfall are known to be different, studies quantifying the contribution of snowmelt to winter drainage discharge are limited. This study aims to investigate the contribution of snowmelt to subsurface drainage discharge at two winterized experimental sites in southern Quebec and Ontario, Canada, using observed data and the Root Zone Water Quality Model-Simultaneous Heat and Water model (RZ-SHAW model). The RZ-SHAW model was calibrated and validated against the measured snow depth and year-round subsurface drainage discharge data from 2021 to 2023 at St. Emmanuel, Quebec, and from 2000 to 2003 at Harrow, Ontario. RZ-SHAW demonstrated satisfactory Nash-Sutcliffe model efficiency values (NSE = 0.50) for simulating snow depth, soil temperature, and winter subsurface drainage discharge. The calibrated RZ-SHAW model was used to simulate the contribution of snowmelt to subsurface drainage for the two sites from 1990 to 2022. The snowmelt was predicted to contribute 55% and 44% of the winter and annual subsurface drainage discharge for the cropland in Southern Quebec. In contrast, for the southern Ontario site, the contribution of snowmelt to winter and annual subsurface drainage discharge was 29% and 18%, respectively. Considering that snowmelt in Southern Quebec contributed a significant fraction of winter and annual subsurface drainage discharge, more attention should be devoted to adapting and developing new winter management for nutrient loss in future research in the region.
Late winter/early spring has been recognized as a critical period for surface runoff, subsurface drainage discharge, and the associated phosphorus (P) loss in croplands in cold climate regions in North America ("Dfa" and "Dfb" in the Köppen–Geiger climate classification system, Peel et al., 2007) due to the melt of snow accumulated across the winter as well as rainfall events occurring in the winter/spring period (Jamieson et al., 2003; Morrison et al., 2014; King et al., 2015; Ahmed et al., 2022; Ren et al., 2022; Ruggiero et al., 2022). Lam et al. (2016) reported that dissolved reactive P (DRP) loss in tile flow during the non-growing season (NGS) induced by snowmelt accounted for 50% of the annual DRP for cropland located in southern Ontario. Similarly, in Ontario, Plach et al. (2019) observed that DRP and TP losses during NGS account for 76% and 78% of the annual DRP and total P (TP) losses for three farms in southern Ontario. Coelho et al. (2012) reported that NGS delivered a more significant DRP load than GS for two fields south of Ontario. In Quebec, the snowmelt period contributed 99.8% of annual surface runoff and 96.7% of surface P load, and the subsurface flow during snowmelt accounted for 52% of the annual total runoff in 2000-2001 (Jamieson et al., 2003). In southern Quebec's Pike River watershed, 82% of total P discharge was linked to post-harvest rainfall events, while the TP loss from subsurface drainage was 37% (98 g/ha) of the total P load during snowmelt (Adhikari et al., 2010). Furthermore, subsurface drainage has been reported to be the dominant hydrologic pathway (42% to 97% of the total runoff) in multiple studies in the more temperate regions in Canada (Great Lakes–St. Lawrence Lowlands region) (Enright and Madramootoo, 2004; Eastman et al., 2010; Tan and Zhang, 2011; Zhang et al., 2017; Kokulan et al., 2019; Plach et al., 2019).
Surface runoff generated by rainfall and snowmelt conditions shares different characteristics in the detachment and transport of P (Pensuka et al., 2008). Runoff caused by snowmelt could be highly selective in removing fine-textured and low-specific materials (Spomer and Hjelmfelt, 1983; Chanasyk and Woytowich, 1987). On the contrary, runoff events generated by rainfall contain higher energy, facilitating an overall greater soil erosion rate than the runoff generated by snowmelt, but are less selective during particle detachment and transport. For two agricultural watersheds located in eastern Canada, Su et al. (2011) reported that rainfall-runoff events had been observed to lead to higher soil erosion and particulate phosphorus (PP) loss than snowmelt-runoff events during the snowmelt period. Panuska et al. (2008) reported snowmelt-induced surface runoff to reveal a higher concentration of dissolved P compared to the rainfall runoff for a terraced hillslope in Wisconsin. Hoffman et al. (2019) reported that snowmelt and the mixed precipitation-driven edge of field runoff mobilized higher dissolved P compared to rainfall runoff for a farm in Wisconsin. Nonetheless, under a changing climate, precipitation as snow may be reduced, and increased rainfall may occur (Henry, 2007; Qian et al., 2011; Bouchard and Qi, 2016; ECCC, 2018a); the volume of subsurface drainage and the corresponding P loss dynamics (Van Esbroeck et al., 2016). Currently, limited research exists in cold regions quantifying snowmelt and rainfall to subsurface drainage during winter and spring, where subsurface drainage has been recognized as the dominant runoff pathway for subsurface-drained croplands (Enright and Madramootoo, 2004; Eastman et al., 2010; Kokulan et al., 2019; Plach et al., 2019).
It is technically challenging to measure the infiltration of melted snow into a frozen soil profile covered by snow and ice; furthermore, precipitation and snowmelt could coincide (rain on snow events, ROS). On the other hand, the hydrological subroutine of calibrated and validated agricultural system models provides an alternative way to quantify hydrologic cycle components efficiently (Ma et al., 2005), which are difficult to measure. The Root Zone Water Quality Model 2 (RZWQM2) is an agricultural system model that considers hydrology, plant nutrition and growth, pesticide transport, and agricultural practices. RZWQM2 has been verified and tested adequately for simulating hydrological balance, greenhouse gas emissions, water quality, and crop growth under various management practices (Ma et al., 2012).
Accurate subsurface drainage simulation is crucial since it is directly related to the simulation of the associated nutrient loss in the subsurface drainage (Qi and Qi, 2017; Jiang et al., 2018; Sadhukhan et al., 2019; Boico et al., 2022; Shokrana et al., 2022). In simulating the hydrologic cycle under cold climates, RZWQM2 provides two options: a generic snowmelt subroutine: Precipitation Runoff Modeling System (which does not consider frozen soil conditions; Leavesly et al., 1983) and the SHAW freeze-thaw subroutine (RZ-SHAW) (Flerchinger et al., 2000). Studies have reported that the RZWQM generic subroutine does not work well in simulating subsurface drainage during snowmelt. Jiang et al. (2018) reported that the RZWQM2 generic subroutine vastly overestimated the subsurface drainage flow in March by 70 mm (two years in total) for a corn-soybean field in southern Quebec from 2008 to 2009. The subsurface drainage was underestimated by 20.8 mm by RZWQM in April 2009. It was suspected that incorrect modeling of snowmelt dynamics in the generic freeze-thaw subroutine caused poor simulation performance. A similar conclusion was reached by Ahmed et al. (2007) that the RZWQM2 generic subroutine was not accurately predicting spring snowmelt hydrology.
Limited research exists evaluating RZ-SHAW's performance in simulating subsurface drainage during winter periods. Jiang et al. (2020) used the RZ-SHAW model to simulate 3-year subsurface drainage at a site in southern Ontario under free drainage (FD) and controlled drainage (CD) treatments. Although RZ-SHAW generally performed satisfactorily in simulating subsurface drainage across the 3-year period (Nash-Sutcliffe efficiency, NSE, of 0.62), it was reported to underestimate subsurface drainage in the winter under both FD and CD treatments. From October 2008 to February 2009, RZ-SHAW underestimated approximately 160 mm (89%) under FD and 80 mm (57%) under CD. The underestimation was suspected to be primarily a result of underestimated infiltration over frozen soil. However, in Jiang et al. (2020), RZ-SHAW was calibrated and validated only using drainage flow and surface runoff across the years (2008-2011) and was not validated with additional measured winter data such as snow depth and soil temperature. As a result, it is unclear whether the underestimation of winter drainage flow by RZ-SHAW resulted from an inaccurate simulation of infiltration over frozen soil or the incorrect simulation of snow accumulation and the soil freeze-thaw dynamics. On the other hand, although calibration for the soil hydraulic properties was established against the observed subsurface drainage in the study, the lateral saturated hydraulic conductivity (LKsat) was not calibrated. Nonetheless, LKsat has been identified as one of the critical parameters in delivering accurate simulations of subsurface drainage in RZWQM (Ma et al., 2012).
Accordingly, the objective of the current project is twofold:
1. Evaluate the performance of RZ-SHAW in simulating winter subsurface drainage discharge under cold climates using subsurface flow and surface snow cover data from two winterized sites in Quebec and Ontario, Canada.
2. Predicting the 30-year snowmelt contribution to the winter and spring snowmelt season subsurface drainage discharge for the two winterized sites in Quebec and Ontario, Canada, through the calibrated RZ-SHAW scenarios.
Methodology
RZ-SHAW Model Overview
RZWQM2 is a one-dimensional agricultural system model capable of simulating the movement of water and nutrients and pesticide fates in agricultural soils under alternative management practices (Ahuja et al., 2000; Ma et al., 2012). Flerchinger et al. (1999) integrated simple routines for snow accumulation and melting into RZWQM. Later, to allow RZWQM to address winter hydrology and soil freeze-thaw dynamics fully, Flerchinger et al. (2000) incorporated the Simultaneous Heat and Water (SHAW) model into RZWQM. The current version of RZWQM2 includes SHAW as an option (noted as RZ-SHAW) to simulate the energy balance in the snowpack (RZ-SHAW defined, approximately 2.5 cm each layer) and soil layers (user-defined), soil freezing and thawing, snow depth, and soil temperature. The soil water movement routine is simulated by RZWQM2 using the Green-Ampt equation for the infiltration process and the Richards equation for soil liquid water distribution (Ahuja et al., 2000). RZWQM2 provides the SHAW model with the evaporation flux, soil liquid water content, and plant information at each timestep. On the other hand, SHAW provides RZWQM2 with the soil temperature, adjusted liquid water content based on the freeze-and-thaw calculation, and energy balance components (total net radiation, ground heat flux, latent heat, and sensible heat) (Li et al., 2012; Ma et al., 2012). SHAW routines have hourly or sub-hourly timesteps during the soil water redistribution. The liquid water content in frozen soil is computed through the Clausius-Claperyon equation, where the total water potential is related to the soil temperature (Flerchinger et al., 2012). Any soil water content greater than the computed liquid water content is considered as soil ice content.
In the RZ-SHAW model, subsurface drainage is calculated from Hooghoudt's steady-state equation as applied by Skaggs (1978) for the DRAINMOD. Precipitation is assumed to be snow if the air temperature is below 0°C. During snowmelt, the excess liquid water is routed through the snowpack through lag and attenuation coefficients, and the current RZ-SHAW does not consider the permeability of the snow layer, but a set of lag attenuation empirical equations was adopted to determine the timing for the snowmelt water outflow. The SHAW model computes frozen soil temperature by solving the vertical temperature distribution in the soil, considering the convective heat transport by liquid and latent heat transfer by vapor (Flerchinger et al., 2000). Details of incorporating the RZ-SHAW and SHAW models into RZWQM2 are described by Flerchinger et al. (2000).
The Experimental Sites
The current modeling study used data collected from two winterized experimental sites (table 1). The first site is a 4.2-hectare subsurface-drained field on the St. Emmanuel Road near Coteau-du-Luc, QC, in the St. Lawrence lowlands region (referred to as the St. Emmanuel site) (Madramootoo et al., 1994; Tait et al., 1995). The second site was the Eugene F. Whalen Experimental Farm, established by Agriculture, and Agri-food Canada, located near Harrow in Southwestern Ontario (referred to as Harrow site) (Tan et al., 1993). Subsurface drainage pipes are extended into walk-in half-basement style buildings with room temperature controlled by electric heaters in winter (fig. 1). At the St. Emmanuel site after, flow was measured and water samples were collected, water was drained out through a sink on the floor to an external pump station and disposed of in an open ditch nearby through underground pipes. In the current study, the winter season was defined as from the first day of November to the last day of April to cover the time span with the existence of snow for the two sites (Li et al., 2022) instead of from the winter solstice to the spring equinox. The spring snowmelt season was defined as the first day when significant snowmelt occurred during February or March (when snow depth peaked and commenced decreasing) towards the last day when the snowpack was melted entirely, plus ten days as a grace period for the snowmelt water to be fully discharged through drains. Ten days were used as the grace period as the drainage outflow ceased after nine days when the snowpack melted entirely at the St. Emmanuel site in 2021.
| Table 1. General information about the two study sites. | ||||
| Site/Variable | St. Emmanuel Site, Quebec | Harrow Site, Ontario | ||
| Coordinates | (74°11'W, 45°21'N) | (82°44'W, 42°13'N) | ||
| Average annual temp. | 6.7 | 8.9 | ||
| Long-term average annual precipitation (1990-2022) (mm) | 845 | 880 | ||
| Field experiment duration: | 2020-2022 | 1999-2003 | ||
| Number of research plots | 24 plots (15m * 75m each) | 16 plots (67.1 m × 15.2 m) | ||
| Plots used in the current study | Plot 1-1, 1-2, 1-3, 1-8, 2-1, 2-2, 2-3, 2-4, 3-5, 3-6, 3-7, 3-8 | Plot 5, 9 | ||
| Drain pipes' depth, spacing, and diameter | Depth: 1.0 m Spacing:15 m Diameter: 76 mm | Depth: 0.6m Spacing: 7.5 m Diameter: 102 mm | ||
| Soil classification | Soulanges very fine sandy loam | Brookston clay loam soil | ||
| Texture of top 25 cm soil (%) | 49 sand 19 silt 32 clay | 28 sand 35 silt 37 clay | ||
| SOC (%) | 5 | 3.6 | ||
| PH. | 7.7 | 7 | ||
| Bulk density of top 25 cm soil (g/cm3) | 1.6 | 1.3 | ||
| Depth to the restrictive layer (cm) | 150 | 300 | ||
Experimental data for the St. Emmanuel site was collected from November 2020 to April 2023. The soil in the St. Emmanuel site was described as a stone-free Soulanges very fine sandy loam (Tait et al., 1995). Tait et al. (1995) and Madramotoo et al. (2001) described each layer's soil physical properties for the site. The field in the St. Emmanuel site was overall flat, with an average slope of less than 0.5% (Kaluli et al., 1999). The field has been artificially subsurface-drained since 1992 (Tait et al., 1995). During the field experiment for the current study, corn was planted each year at a row spacing of 0.76 m. The farmer seeded at a density of 89,000 seeds ha-1 each year. Corn was planted on 28 April 2020, 27 April 2021, and 8 May 2022. The fertilizer application rate was identical for nitrogen and potassium in 2020 and 2021, 200 kg/ha and 45 kg/ha. In 2022, the applicate rate for nitrogen and potassium was 220 kg/ha and 40 kg/ha, respectively. The fertilizer application rate for phosphorus was 42 kg/ha in 2020 and 25 kg/ha in 2021 and 2022. Two fertilizer applications each year were made: (1) broadcasted along with the tillage in spring, and (2) band application during seeding. The field was tilled with a chisel plow after tooth harrow before seeding in spring.
The St. Emmanuel site contains 24 subsurface-drained plots (fig. 1), grouped into three blocks,each housing two water table management practices. The two water table management practices were: free drainage (FD) and controlled drainage (CD) (12 plots for FD, 12 plots for CD) (fig. 1). In the current study, only FD will be considered for simulation. The subsurface drains were installed at a depth of 1.0 m at a slope of 0.5% in the center of each plot with a spacing of 15 m, running the plot's length (Madramootoo et al., 2001). The diameter of the drainpipe was 76 mm (3-inch), and they were wrapped in geotextile to prevent clogging (knitted sock). The drainage water from each drain was discharged into one of the monitoring water houses located between the blocks. Each block was separated by a 30-m wide strip (Tait et al., 1995).
 |
| Figure 1. Experimental layout and winterized facilities at the St. Emmanuel site (adapted from Tait et al., 1995). Upper right picture shows the winterized house during winter. Bottom right picture shows the drainage discharge and water sample collection facility. Only drainage data from the free drainage plots were used in the current study (marked in red text). |
The field experiment at the Harrow site in the north watershed of Lake Erie was conducted from 1999 to 2003. The field has a slope between 0.05% and 0.1% and a Brooklyn clay loam soil (Tan et al., 2002, table 1). There were sixteen plots (67.1 m × 15.2 m each) with various treatment combinations, including water table management (FD vs. CD), compost amendments (yard waste leaf compost and liquid pig manure compost), and cover crops (no cover crops vs. winter wheat as a cover crop, as shown in fig. 2). Two tiles (drain spacing of 7.5 m) were located at 0.6 m depth in each plot. The tiles ran parallel to the length of each plot. Each plot had a 30 cm high berm and a 1.2 m deep plastic barrier installed on three sides (excluding the lower end with the catch basin) to prevent surface and subsurface water movement between plots. In the current study, only drainage data from plots 5 and 9, which housed free drainage and no cover crop treatments, were averaged and used to evaluate the RZ-SHAW model. Corn was planted in 2000 and 2002, while soybean was planted in 2001 and 2003. Detailed information on the experimental design and agronomic management information has been described by Zhang et al. (2017).
Meteorological Data and Field Data Retrieval
The meteorological data, including the snow depth, relative humidity, wind speed, and air temperature, were retrieved from the Environment Canada (EC) meteorological stations closest to each site. The EC Harrow AUTO CDA weather station (station ID 6133362) was located 0.5 km from the Harrow site, whereas the EC Coteau-du-Luc station (station ID 7011947) was located 0.7 km from the St. Emmanuel site (ECCC, 2021). Total precipitation for each site was retrieved from the Adjusted Daily Rainfall and Snowfall Data for Canada dataset provided by Environment Canada for the corresponding EC weather stations (Wang et al., 2017). The total precipitation was adjusted to address gauge under catch and evaporation due to the wind effect for gauge-specific wetting loss. The solar radiation data for each site was retrieved from the NASA-MERRA platform (Gelaro et al., 2017).
Soil temperature for the St. Emmanuel site was measured at four depths (2 cm, 5 cm, 10 cm, and 20 cm) using temperature probes (temperature probes 107, Campbell) for three winter seasons at a 15-minute resolution. The soil temperature data was recorded using the Campbell 1000x datalogger. For the first winter season, soil temperature was measured from December 2020 to April 2021. The soil temperature measured for the second and third winter seasons was from November 2021 to April 2022 and November 2022 to April 2023. No soil temperature data were measured for the Harrow site.
Subsurface drainage from each plot in the St. Emmanuel sites was directed to a winterized house. It was recorded for every 1 L of flow using pre-calibrated tipping buckets (Tait et al., 1995). Drainage outflow for the twelve plots of the twenty-four total plots housing free drainage (FD) in the St. Emmanuel site was used in the current study. The average of the four plots with the outflow closest to the mean daily discharge of all twelve plots was used to evaluate the RZ-SHAW model. From late fall 2020 to winter 2021, the data logger for the first eight FD plots was not functioning, and only the data from the last four plots (plots 3-5, 3-6, 3-7, and 3-8) were used (from the second water house). The two warehouses were flooded from 19 to 22 March under snowmelt and heavy rainfall. The data for the flooding period was removed during the validation.
 |
| Figure 2. Experimental layout and winterized facilities of the Harrow site. Plots in red dashed line are the plots used in the current study. |
For the Harrow site, the subsurface drainage water was routed from the plots to a winterized monitoring building using non-perforated tiles, which contained catch basins for both surface runoff and subsurface drainage for each plot (Zhang et al., 2017). A multi-channel data logger continuously monitored and recorded drainage flow rates. In addition, a backup generator was present and automatically engaged in the event of electric power interruptions to ensure uninterrupted year-round data and water sample collection. Additional detailed settings for the data collection have been previously described by Tan et al. (1993). In the current study, surface runoff and subsurface drainage data from only plots 5 and 9, which housed free drainage and no cover crop, were averaged and used to evaluate the RZ-SHAW model (fig. 2). Surface runoff and subsurface drainage water samples were collected using auto-samplers (CALPSO 2000S, Buhler Gmbh & Co., Uzwil, Switzerland) after every 500 to 3000 L of flow, depending on the time of year, on a year-round continuous basis (fig. 2, table 2).
| Table 2. Drainage flow data period for the Hallow site. | |
| Period | Duration |
| 1 | 18 May 2000 – 31 May 2000 |
| 2 | 31 May 2000 – 28 June 2000 |
| 3 | 28 June 2000 – 24 Aug. 2000 |
| 4 | 24 Aug. 2000 – 20 Sep. 2000 |
| 5 | 20 Sep. 2000 – 03 Oct. 2000 |
| 6[a] | 03 Oct. 2000 – 29 Nov. 2000 |
| 7[a] | 29 Nov. 2000 – 13 Feb. 2001 |
| 8[a] | 13 Feb. 2001 – 11 June 2001 |
| 9 | 13 June 2001 – 16 Aug. 2001 |
| 10 | 16 Aug. 2001 – 21 Nov. 2001 |
| 11[a] | 20 Nov. 2001 – 21 Nov. 2002 |
| 12[a] | 21 Nov. 2001 – 16 Jan. 2002 |
| 13[a] | 16 Jan. 2002 – 06 Mar. 2002 |
| 14[a] | 06 Mar. 2002 – 16 Apr. 2002 |
| 15 | 16 Apr. 2002 – 22 May 2002 |
| 16 | 22 May 2002 – 27 June 2002 |
| 17 | 27 June 2002 – 05 Nov. 2002 |
| 18[a] | 05 Nov. 2002 – 11 May 2003 |
| 19 | 11 May 2003 – 18 June 2003 |
| 20 | 18 June 2003 – 08 Oct. 2003 |
[a]Indicates the period is in the winter season. | |
Model Simulation
An RZ-SHAW modeling scenario was created for each site. The RZ-SHAW scenarios were initiated with the weather input (wind speed, precipitation, minimum air temperature, maximum air temperature, and solar radiation) in daily resolution. Each site's simulation period was based on the available observed winter subsurface drainage with a three-year "spinup". The soil hydraulic properties for the two sites for each soil horizon were adapted from the calibrated values reported by Jiang et al. (2019) for the St. Emmanuel site and Jiang et al. (2020) for the Harrow site. For the St. Emmanuel site, 27 November 2020, to 21 April 2021, was treated as the calibration period, and 15 November 2021, to 17 April 2023) was treated as the validation period. For the Harrow site, 2001 was treated as the calibration period for the winter subsurface drainage (periods 6-8), while 2002 to 2003 was treated as the validation period (periods 11, 12, 13, 14, and 18). The effective drain radius for the St. Emmanuel site and Harrow sites was 4.4 cm and 0.5 cm based on previously reported values and recommendations (Singh et al., 1994; Ghane, 2022).
Li et al. (2022) described detailed auto-calibration procedures for the RZ-SHAW model. Calibration of initial snow density, maximum average snowpack density, dry soil albedo, wet soil albedo, surface resistance, and lateral hydraulic conductivity was conducted against observed snow depth, soil temperature, and subsurface drainage discharge through a customized auto-calibration script. The script would generate a Sobol sequence for each calibrating variable mentioned above based on the given maximum iterations and range of the targeted parameters. RZ-SHAW was then operated under all combinations in the Sobol-generated sequences to discover the values that deliver the best statistical performance based on the Nash–Sutcliffe model efficiency coefficient (NSE) value (Sobol, 1967; PyTorch, 2019). The auto-calibration script continued to iterate unless the maximum number of iterations was achieved or the NSE value was greater than or equal to 0.65. Model auto-calibration was carried out in a distributed parallel system consisting of three computer workstations with 40 cores to reduce computation time.
Multiple studies have revealed that the correct simulation of the snow cover dynamics was the premise of an accurate simulation of the soil freezing dynamics (Zhang et al., 2000; Jégo et al., 2014; Qi et al., 2016). Thus, the scenario for each site was first calibrated and validated against the observed snow depth to ensure a correct simulation of snowmelt and insulation effects on soil through calibrating the snow density parameters: new-fallen snow density and average maximum snowpack density. Snow depth also had the most significant impact in delivering an accurate simulation for soil temperature in the SHAW model (Flerchinger et al., 2012). The newly fallen snow density and average maximum snowpack density were initialized with the default value provided with RZWQM of 100 kg/m3 and 600 kg/m3. The range of the calibration for average snowpack (100-1000 kg/m3) and initial new snow densities (40-100 kg/m3) was based on the observed minimum and maximum values from Ontario (Whiteley, 2004; ECCC, 2018a). The maximum iteration for the snow depth auto-calibration was set to 4000 times, and parameter values with the highest NSE were saved as the final calibration results (table 3).
After calibration of the snow parameters, the soil's dry albedo, wet albedo, and surface soil resistance would be calibrated against the observed surface layer winter soil temperature. However, this calibration step only applied to the St. Emmanuel site, as no observed soil temperature was available for the Harrow site. On the other hand, considering that winter was the focused period in the current study, calibration in soil albedo has little influence on the winter soil temperature simulation (the field was covered by snow during the winter). The albedo of snow cover was computed through the RZ-SHAW model using an empirical equation based on surface layer snow density. The range of the calibration for dry and wet soil albedo (0.05- 0.45) and soil surface resistance (150-250 s/m-1) was based on the suggestion from Ma et al. (2011). The maximum iteration for the soil temperature auto-calibration was set to 10,000 times.
Once the model has been calibrated for both the snow depth and soil temperature, calibration for the lateral saturated hydraulic conductivity (LKsat) against the observed winter subsurface drainage has been established (VKsat for the two sites has been calibrated in previous literature: Jiang et al. (2019) and Jiang et al. (2020)). The LKsat was initialized with two times the vertical soil saturated hydraulic conductivity (VKsat) as suggested by Ma et al. (2012) for each soil horizon. The maximum value of the LKsat was set to triple the VKsat, and the minimum value for the LKsat was set to the VKsat, within a given layer. The maximum number of iterations was set to 50000 times. Nonetheless, depending on the testing value, RZ-SHAW was computationally intensive and required 10 to 20 minutes to perform a one-time simulation for a 5-year model scenario. Comparatively, the RZWQM model with the default generic snowmelt subroutine (SHAW model turned off) was five to ten times faster than the RZ-SHAW model. Thus, the 50000 iterations were performed based on the RZWQM model first to reduce computation time. Then 2000 sets of values with the best NSE results from the RZWQM's 50000 calibrations were simulated again using the RZ-SHAW model to discover the best values for the soil lateral hydraulic conductivity (table 3).
| Table 3. Calibrated lateral hydraulic conductivity for the St. Emmanuel site and the Harrow site. | ||||||||||
| St. Emmanuel | Harrow | |||||||||
| Soil type | Depth | LKsat (cm h-1) | Initial snow density (kg/m3) | 60 | Soil type | Depth | LKsat (cm h-1) | Initial snow density (kg/m3) | 86 | |
| Sandy clay loam | 0.00-0.05 | 2.6 | Average maximum snowpack density (kg/m3) | 300 | Loam | 0.0-0.25 | 1.7 | Average maximum snowpack density (kg/m3) | 258 | |
| Loam | 0.05-0.25 | 4.2 | Soil dry albedo | 32 | Clay loam | 0.25-0.45 | 5.5 | Soil dry albedo | 0.3[a] | |
| Sandy loam | 0.25-0.45 | 3.2 | Soil wet albedo | 0.23 | Clay loam | 0.45-0.80 | 6.0 | Soil wet albedo | 0.2 | |
| Clay loam | 0.45-0.80 | 3.7 | Surface soil resistance | 210 | Clay loam | 0.80-1.20 | 3.7 | Surface soil resistance | 200[a] | |
| Clay loam | 0.80-1.2 | 4.8 | Clay loam | 1.20-3.00 | 1.5 | |||||
| Clay loam | 1.20-1.50 | 1.5 | ||||||||
| Clay loam | 1.50-2.0 | 2.0 | ||||||||
[a]Indicates the value was not calibrated. For the Harrow site, no observed soil temperature was available, and default values for dry soil albedo, wet soil albedo, and surface soil resistance were provided by the RZWQM model were used. | ||||||||||
Computation of 30-Year Snowmelt and Rainfall Contribution to Drainage
The calibrated and validated RZ-SHAW scenarios were used to investigate the contribution of snowmelt and rainfall towards drainage during the winter or spring snowmelt season. The calibrated scenarios were extended to a 30-year simulation (1990-2022) to investigate the long-term snowmelt contribution to drainage. The meteorological input for the St. Emmanuel site was retrieved through the same Environment Canada weather station (Coteau du Lac, 7011947). For the Harrow site, the meteorological input was retrieved from the same station, but the name and code changed across the years (1990-1993: Harrow CDA 6133360; 1993-2000: Harrow Automatic Climate Station, 613CC60; 2000-2022: Harrow Auto CDA, 6133363).
To estimate the contribution of snowmelt towards winter drainage, the calibrated RZ-SHAW scenario for each site was simulated again with winter precipitation falling as rain removed ("no-winter-rainfall" scenario). The winter drainage simulated under the "no-winter-rainfall" scenario was considered snowmelt-contributed drainage. The removal of winter rain was established based on the current RZ-SHAW's mechanic in differentiating snowfall and rain, which was based on whether the air temperature was below 0°C. However, in reality, rainfall and snowmelt could both contribute to evaporation and surface runoff during rain-on-snow (ROS) events. Thus, for "no-winter-rainfall" scenarios, the snowmelt may bear all the evaporation and surface runoff during ROS events, and the snowmelt-contributed drainage could be underestimated. The current version of the RZ-SHAW model could not differentiate between rainfall and snowmelt contributions to surface runoff and evaporation during ROS events. It was assumed that the contribution of snowmelt and rainfall during ROS events in evaporation and surface runoff was identical. Half of the evaporation and surface runoff during ROS events was added to the snowmelt-contributed drainage to compensate for the potential overestimation of snowmelts' contribution to evaporation and surface runoff during ROS events under "no-winter-rainfall" scenarios. The difference in the drainage discharge between the snowmelt-contributed drainage and the normal scenario (with winter rainfall) was considered the rainfall-contributed subsurface drainage discharge.
Model Performance Evaluation
Four statistical methods were used to evaluate the performance of RZ-SHAW in simulating snow depth, soil temperature, and subsurface drainage relative to the observed values: the Nash-Sutcliffe model efficiency (NSE), mean bias error (MBE), kling-gupta efficiency (KGE), and index of agreement (IOA), given as:
(1)
(2)
(3)
(4)
where
n = number of paired observed and simulated values
OBSi = ith observed value
= mean observed value
sOBS = standard deviation of the observed value
SIMi = ith predicted value
= mean predicted value
SSIM = standard deviation of the predicted value.
An NSE value of 1.0 is considered to be a perfect fit, an NSE value greater than 0.75 is considered to be a very good fit, an NSE value ranging from 0.64 to 0.74 is a good fit, and an NSE value ranging from 0.5 to 0.64 is a satisfactory fit (Moriasi et al., 2015). The value of KGE varies from -8 to 1, with KGE = 1 being a perfect fit between the observed and simulated values (Gupta et al., 2009). A KGE value greater than 0.5 can be considered acceptable simulation performance (Jeantet et al., 2021). The IOA determines the degree of model prediction error, ranging from 0 to 1, with 1 representing the optimal match (Jiang et al., 2019). An IOA greater than 0.7 can be considered a satisfactory model performance (Ma et al., 2011).
Results
Winter Snow Depth, Soil Temperature, and Subsurface Drainage Discharge Simulation
In general, the model performance was satisfactory during both calibration and validation in simulating snow depth, winter soil temperature, and winter subsurface drainage discharge with NSE = 0.5 (except for: snow depth during validation for the St. Emmanuel site; snow depth and winter subsurface drainage discharge during calibration for the Harrow site), KGE > 0.5 (except for: winter subsurface drainage discharge and 20 cm winter soil temperature during validation for the St. Emmanuel site; snow depth during calibration for both sites), and IOA > 0.7 (except for: snow depth during calibration for the Harrow sites) (table 4). The NSE, MBE, KGE, and IOA values for the three simulating variables during calibration and validation were presented in table 4. The simulated and observed daily winter soil temperature, snow depth, and subsurface drainage discharge for the St. Emmanuel site are shown in figure 3. The simulated and observed daily winter snow depth and periodic subsurface drainage discharge for the Harrow site are shown in figure 4.
| Table 4. RZ-SHAW's statistical performance in simulating snow depth, winter soil temperature, and winter subsurface drainage discharge during validation for the St. Emmanuel site and the Harrow site. | ||||||||||
| Site | Variable Type | Calibration | Validation | |||||||
| NSE | MBE | KGE | IOA | NSE | MBE | KGE | IOA | |||
| St. Emmanuel | Daily winter subsurface drainage discharge | 0.64 | -1.2 mm | 0.58 | 0.76 | 0.51 | 0.6 mm | 0.46 | 0.81 | |
| Daily snow depth | 0.61 | -9.3 cm | 0.46 | 0.82 | 0.49 | -6.5 cm | 0.54 | 0.76 | ||
| Daily soil temperature 2 cm | 0.85 | -0.1 °C | 0.61 | 0.97 | 0.88 | 0.1 °C | 0.57 | 0.94 | ||
| Daily soil temperature 5 cm | 0.71 | -0.3 °C | 0.64 | 0.98 | 0.86 | -0.2 °C | 0.75 | 0.96 | ||
| Daily soil temperature 10 cm | 0.70 | -0.3 °C | 0.68 | 0.97 | 0.83 | 0.06 °C | 0.67 | 0.95 | ||
| Daily soil temperature 20 cm | 0.55 | -0.5 °C | 0.51 | 0.93 | 0.61 | 0.4 °C | 0.42 | 0.92 | ||
| Harrow | Periodic winter subsurface drainage discharge | 0.23 | -2.0 mm | 0.56 | 0.74 | 0.61 | -7.4 mm | 0.94 | 0.77 | |
| Daily snow depth | 0.17 | -3.2 cm | 0.22 | 0.63 | 0.50 | -2.8 cm | 0.56 | 0.76 | ||
During calibration of the snow depth, the NSE value for the St. Emmanuel site reached satisfactory (0.61), and the NSE value was unsatisfactory (0.17) for the Harrow site. For the KGE value, the simulation performance during calibration was inadequate for both sites (0.46 for the St. Emmanuel site and 0.22 for the Harrow site). The IOA value for the St. Emmanuel site reached satisfactory (0.82) but was unsatisfactory for the Harrow site (0.63). The unsatisfactory simulation for the Harrow site was primarily a result of the significantly underestimated snow depth in the year 2000 (fig. 4). The observed snow depth was 12.2 cm, while the simulated average was 5.6 cm from January to March 2000. In 2001, an accurate simulation was delivered before 19 January, but RZ-SHAW overestimated the snow depth by 13.4 cm from 20 to 31 January. RZ-SHAW tended to simulate faster snow depth decline (which may indicate simulated faster snowmelt) than the observed snow depth for the St. Emmanuel site. In 2021, significant snow depth reduction occurred on 6 March in the observed data, whereas RZ-SHAW simulated the snowmelt to occur on 2 March. Similarly, in 2022, major snowmelt in the observed snow depth occurred on 4 March, but the RZ-SHAW simulated the snowmelt to appear on 25 February. Aside from the major snowmelt events, RZ-SHAW simulated the snow depth collapse immediately after snowfall events, when no subtle changes were observed in the measured snow depth. For instance, from 6 to 15 February in 2021, RZ-SHAW simulated the snow depth to decrease progressively, whereas there was little change in the observed snow depth. The simulated faster snowmelt also led to the failure to simulate the peak snow depth on 19 February. During validation, NSE values for simulating the snow depth for the St. Emmanuel and Harrow sites were 0.49 and 0.50. The KGE and IOA values for the snow depth were satisfactory during validation for both sites.
All model performance indexes were satisfactory during the calibration period for the winter soil temperature simulation at the St. Emmanuel site. During validation, NSE values of 0.88, 0.86, 0.83, and 0.61 at depths of 2 cm, 5 cm, 10 cm, and 20 cm were achieved during validation (fig. 3). For the KGE values, only the KGE value for the winter soil temperature at 20 cm was unsatisfactory (0.42). Winter soil temperature simulations across all depths achieved IOA values greater than 0.7.
For the winter drainage discharge simulation during the calibration in the St. Emmanuel site, the NSE, KGE, and IOA values were satisfactory. RZ-SHAW successfully captured the peak drainage during the snowmelt season at the end of March 2021. Still, it overestimated subsurface drainage in April 2021 (fig. 3). The observed peak drainage discharge (defined as daily observed drainage greater than 10 mm day-1) lasted ten days, while the RZ-SHAW simulated peak drainage discharge lasted only four days. From 18 March to 4 April, the observed total drainage discharge was 185 mm, while RZ-SHAW simulated 125 mm. On the other hand, RZ-SHAW overestimated the mid-April drainage discharge during the winter season of 2021, when there was little (less than 0.1 mm day-1) drainage in the observed data. RZ-SHAW simulated the drainage to occur ten days earlier than the observed drainage in the spring of 2021. In total, the observed subsurface drainage discharge during the winter season of 2021 was 218 mm, and the simulated subsurface drainage discharge was 162 mm.
The statistical performance of the winter subsurface drainage during the validation phase was worse than the calibration phase for the NSE and KGE values (NSE: 0.64 to 0.51, KGE: 0.58 to 0.46). The date when daily drainage discharge was more significant than 1 mm occurred was 22 February in the simulated result and 15 March in the observed data (21 days apart) in 2022. The difference in the occurrence of 1 mm of subsurface drainage was 13 days from the end of December 2022 to the beginning of January 2023 (simulated: 23 December 2023, observed: 23 March 2023). From 20 February to 20 March 2022, the observed total drainage discharge was 25 mm, whereas the RZ-SHAW simulated 75 mm. From 22 December 2022, to 22 January 2023, the observed total drainage discharge was 31 mm, whereas the RZ-SHAW simulated 87 mm. Overall, RZ-SHAW overestimated the drainage discharge by 23 mm and 59 mm during winter 2022 and 2023 (observed: 178 mm in 2022 and 152 mm in 2023, simulated: 203 mm in 2022 and 212 mm in 2023).
For the Harrow site, only the NSE value was unsatisfactory across all the model performance evaluation indexes during calibration (0.23). All model performance evaluation indexes were improved during validation compared to the calibration period (table 4). An MBE of -7.4 mm for the periodic winter drainage was observed during the validation phase. The underestimation in the subsurface drainage discharge predominantly occurred during the colder periods (January and February). In contrast, the overestimation occurred in the winter periods, which covered warmer times during the winter (November and April) (fig. 4). The most significant underestimation of winter drainage by the model at the Harrow site occurred in period 7 (underestimated by 32 mm compared to the measured value) of the calibration
 |
| Figure 3. RZ-SHAW simulated daily subsurface drainage discharge, winter daily snow depth, and winter daily soil temperature vs. observed values for the St. Emmanuel site. (a) Simulated winter daily subsurface drainage discharge vs. observed subsurface drainage discharge. (b) Observed rainfall and snowfall. (c) Simulated daily snow depth vs. observed snow depth. (d) Simulated winter daily soil temperature at 2 cm vs. observed soil temperature at 2 cm. (e) Simulated winter daily soil temperature at 10 cm vs. observed soil temperature at 10 cm. The drainage discharge data during the 2022 flooding event were removed from the statistical analysis for the St. Emmanuel site. Calibration and validation periods were separated by the black vertical line with the arrow. |
period, which covered drainage from November 2000 to February 2001. Periods 10, 12, 13, and 18 all demonstrated underestimation compared to the observed drainage (underestimated by 5 mm to 20 mm), which covered January to March in 2002 and the winter season of 2003. On the other hand, RZ-SHAW overestimated the subsurface drainage discharge in periods 6 and 8 (October to November in 2000 and February to June in 2001) by 8 and 13 mm. Similarly, RZ-SHAW overestimated the subsurface drainage discharge by 2 mm and 3 mm in periods 10 and 14 (November 2001 and April to May 2002). RZ-SHAW overestimated the observed drainage by 8 and 13 mm.
 |
| Figure 4. RZ-SHAW simulated winter daily snow depth and subsurface drainage discharge in periods vs. observed values for the Harrow site. (a) Simulated daily snow depth vs. observed daily snow depth. (b) Simulated periodic subsurface drainage vs. observed periodic subsurface drainage. Bars with * on the top indicate the period is in the winter season. Calibration and validation periods were separated by the black vertical line with the arrow. |
Contribution of Snowmelt to Subsurface Drainage Discharge
The 30-year (1990-2022) contribution of snowmelt and rainfall towards subsurface drainage discharge for the St. Emmanuel site and the Harrow site during the winter and spring snowmelt season were computed based on the calibrated and validated scenarios and shown in figure 5. Snowmelt was the dominant contributor in both winter and spring snowmelt drainage for the St. Emmanuel site, whereas rainfall was the dominant contributor for the Harrow site. The 30-year average snowmelt drainage in the St. Emmanuel site contributed to 52% of the winter drainage discharge and 69% of the spring-snowmelt period drainage discharge. Snowmelt drainage contributed 44% to the annual drainage discharge at the St. Emmanuel site. Comparatively, in the Harrow site, snowmelt drainage contributed 22% of the winter drainage and 48% of the spring snowmelt drainage. Snowmelt drainage contributed 9% of the annual drainage at the Harrow site.
The estimated 30-year average rainfall and snowmelt's contribution towards drainage during the winter and spring snowmelt period aligned with the observed snowfall and rainfall for the two sites (table 5). On average, the St. Emmanuel site experienced 475 mm of total winter precipitation, and the Harrow site received 451 mm of total winter precipitation. Although the observed winter total precipitation was similar, under the common conversion metric where the water equivalent of newly fallen snow is 10 to 1 (ECCC, 2018b), St. Emmanuel experienced higher snowfall as precipitation (49%) compared to the Harrow site (29%). A similar pattern could also be witnessed in the observed number of days of rainfall, snowfall, and "Rain-on-snow" days (table 5). The St. Emmanuel site also experienced a more extended spring snowmelt season than the Harrow site (35 days vs. 28 days). The difference in the partition between snowfall and rainfall was primarily the result of the winter air temperature. The Harrow site witnessed a higher average winter air temperature of 2°C compared to -3°C.
Discussion
Potential Limitations of RZ-SHAW and Future Improvement
Although RZ-SHAW delivered a satisfactory simulation for the winter subsurface drainage discharge, it tended to overestimate the subsurface drainage discharge in early spring and underestimate the peak drainage discharge (during the period when peak occurred in the observed data),
 |
| Figure 5. RZ-SHAW simulated 30-year rainfall and snowmelt contribution towards winter subsurface drainage discharge. (a) St. Emmanuel site's winter rainfall and snowmelt contribution towards drainage, (b) St. Emmanuel site's spring snowmelt season rainfall and snowmelt contribution towards drainage, (c) Harrow site's winter rainfall and snowmelt contribution towards drainage, (d) Harrow site's spring snowmelt contribution towards drainage, (e) St. Emmanuel site's winter drainage/spring snowmelt season drainage's contribution towards annual drainage, (f) Harrow site's winter drainage/spring snowmelt season drainage's contribution towards annual drainage. |
usually during late winter and early spring. Slight overestimation occurred in the Harrow site in periods 6 and 11, which covered November to December in 2000 and 2001, respectively. On the other hand, RZ-SHAW simulated spring drainage to occur earlier in spring 2022, December 2023, and March 2023 compared to the observed drainage at the St. Emmanuel site. RZ-SHAW underestimated the observed peak drainage at the St. Emmanuel site in both winter 2021 and winter 2022. Multiple potential issues, including the frozen soil hydraulic conductivity estimation, the simulation of the ice layer within the snow cover, and the snow cover simulation, could all potentially contribute to the errors mentioned above.
| Table 5. Observed 30-year winter climate patterns and simulated 30-year winter croplands hydrology. | ||||||||||||||||||||
| Observed 30-Year Average Meteorological Pattern | ||||||||||||||||||||
| Site | Annual Precipitation (mm) | Winter Precipitation (mm) | Winter Rainfall (mm) | Winter Snowfall (mm) | Number of Snowfall Days | Number of Rainfall Days | Number of ROSb Days | Winter Average Air Temperature (°C) | ||||||||||||
| St. Emmanuel | 845 (112[a]) | 475 (86) | 244 (71) | 232 (62) | 51 (11) | 44 (10) | 11 (5) | -3 (1) | ||||||||||||
| Harrow | 880 (116) | 451 (92) | 327 (80) | 131 (44) | 43 (8) | 59 (9) | 7 (3) | 2 (1) | ||||||||||||
Simulated 30-Year Average Cropland Hydrology | ||||||||||||||||||||
| Site | Annual Drainage (mm) | Winter Drainage (mm) | Spring Snowmelt Drainage (mm) | Annual Runoff (mm) | Winter Runoff (mm) | Spring Runoff (mm) | Annual Evapo-transpiration (mm) | Winter Evaporation (mm) | Spring Snowmelt Evaporation (mm) | Spring Snowmelt Starting Date (Julian Day) | Spring Snowmelt Ending Date (Julian Day) | |||||||||
| St. Emmanuel | 272 (91) | 216 (78) | 142 (48) | 63 (37) | 18 (14) | 12 (7) | 536 (38) | 89 (28) | 41 (11) | 64 (18) | 99 (15) | |||||||||
| Harrow | 246 (53) | 135 (34) | 53 (31) | 95 (50) | 38(26) | 17 (11) | 529 (47) | 90 (8) | 46 (6) | 66 (14) | 94 (17) | |||||||||
Note [a]: the number in the parenthesis indicates the standard deviation; b: ROS was the abbreviation for the "Rain-on-snow". | ||||||||||||||||||||
Firstly, the hydraulic conductivity of frozen soil may not be simulated correctly. In the current version of RZ-SHAW, the hydraulic conductivity for frozen soil was linearly reduced with the simulated soil ice content in each layer. When the ice is presented in the soil layer, the maximum pore space available for unfrozen water movement for each layer is calculated by subtracting the saturated water content from the ice content. The hydraulic conductivity for the frozen soil was then linearly reduced based on the hydraulic conductivity for the unfrozen soil with the changing porosity. The hydraulic conductivity for unfrozen soil was computed based on the suction head and water content from Brooks-Corey's curve (RZ-SHAW assumes frozen soil's water retention curve was the same as the dry portion curve) (Ma et al., 2012). If the computed maximum pore space for liquid water movement was less than 0.13, the hydraulic conductivity would be 0 (Flerchinger et al., 2012). Currently, no calibratable impedance factor was implemented in the current RZ-SHAW model.
Comparatively, two other well-known soil hydrology models, HYDRUS-1D and DRAINMOD, either implemented a calibratable impedance factor or a calibratable critical ice content to prevent the water from infiltrating by adjusting the hydraulic conductivity for the frozen soil (Luo et al., 2000; Hansson et al., 2004). Better simulation performance in simulating winter subsurface drainage discharge (NSE of 0.6) for an agricultural field in Quebec using the DRAINMOD model was reported by Morrison et al. (2014). Therefore, it may be worth considering adapting an impedance factor for the RZ-SHAW model in simulating hydraulic conductivity for frozen soil. On the other hand, soil ice content was also the main factor in estimating the hydraulic conductivity of the frozen soil in the current RZ-SHAW model. An underestimation of soil ice content could also contribute to the overestimated hydraulic conductivity of frozen soil. Nonetheless, no observed soil ice content was available for the current study. Similarly, data was rarely available across field studies, and RZ-SHAW's ability to simulate ice content had not been evaluated. Incorrect simulation of the soil ice content may affect the frozen soil hydraulic conductivity estimation, which may partially explain the tendency to simulate early drainage for the St. Emmanuel site in the current study. Thus, it would be necessary for future field experiments to include such measurements to improve the understanding of the hydraulic conductivity of frozen soil.
Aside from the potentially overestimated soil hydraulic conductivity, another potential limitation of the current RZ-SHAW model is the lack of a snow cover permeability module. The ice layer inside the snow could be developed throughout the winter due to ROS and repeated freeze-thaw cycles, which could impede the meltwater in the snow layer from infiltrating into the soil (Wever et al., 2016). The development of the ice layer at the bottom of the snowpack has been shown to be tightly coupled with the preferential flow in the snow cover (Wever et al., 2016). In the current version of the RZ-SHAW model, the permeability of the snowpack was not considered. Instead, a set of lagged and attenuation empirical equations was employed to compute the timing for the excess liquid water from the snow cover to flow out after the snowpack's water-holding capacity was satisfied (Flerchinger et al., 2012). The maximum allowable lag for the excess water from the snowpack was hardcoded inside the model to be 10 hours in the current RZ-SHAW model based on the report from Anderson (1976). Nonetheless, the actual time for the water to percolate through the snowpack and infiltrate into the soil may be longer due to the potential development of the ice layer. Thus, the lack of description of the impedance of the ice layer towards water flow may potentially be one of the contributors to the simulated early drainage, as shown in figure 3.
Simulation of the snow depth dynamics could also be a critical factor for predicting winter subsurface drainage when snow cover presents. For the Harrow site in the winter of 2001, a considerable overestimation of the snow depth occurred at the end of January. Significant underestimation in the drainage was also witnessed in period 7 (November 2000 to January 2001). RZ-SHAW simulated more snowfall as precipitation than rainfall. As it would take a longer time for the snow to melt, such an overestimation in snowfall could shift the portion of the infiltration input that should occur in January to a later date, which could also explain the overestimation in drainage discharge in period 8. Currently, RZ-SHAW hardcoded the air temperature of 0°C as the threshold value in differentiating whether the precipitation was falling as rain or snow. Nonetheless, a recent study by Jennings et al. (2018) revealed that the rain-snow thresholds for different regions in the Northern Hemisphere ranged from -0.4°C to 2.4°C, with an average of 1.0°C. Furthermore, models such as HYDRUS-1D and DRAINMOD implement a range of threshold values to consider rain and snowfall that may coincide on the same day. Thus, RZ-SHAW's snowfall simulation may be improved by implementing a calibratable snow-rain air temperature threshold range based on the observed snow-rain air temperature threshold.
A potential source of error in the snow depth simulation in the current study could be the lack of measured in-situ snow water equivalent data to allow correction for the site-specific precipitation input due to gauge under catch during winter for the two sites. The calibrated initial snow density of 60 kg/m3 for the St. Emmanuel site may be compensating for the gauge under catch. As a result, the rapid decline of snow depth and early snowmelt could be observed compared to the observed snow depth. The simulated early snowmelt due to compensated initial snow density during calibration could also be one of the contributors to the simulated early drainage.
Potential Error in Estimating Rainfall and Snowmelt's Contribution to Drainage
In the current study, "Rain-on-snow" (ROS) events that often occur during early spring could lead to errors in estimating the contribution from the snowmelt and rainfall for subsurface drainage discharge. The current study assumed that snowmelt and rainfall contributed equally towards evaporation and surface runoff during ROS events. On a 30-year average, compensation from evaporation and runoff to snowmelt drainage during ROS during the simulations were 9 mm and 12 mm for the St. Emmanuel site and 15 mm and 16 mm for the Harrow site, respectively. However, infiltration and evaporation during ROS were a dynamic process. Rain that fell on the snow cover could be refrozen and stays in the snowpack. On the other hand, relatively warm rainfall could also lead to the faster melting of the snow cover, which may lead to flooding (Pomeroy et al., 2016). Similar flooding due to ROS during snowmelt occurred in the St. Emmanuel site in March 2022, where the ditch was flooded, and the drainage water could not be drained out. Currently, it is difficult to distinguish the contribution of rainfall and snowmelt to surface runoff during ROS events (Li et al., 2019). Similarly, it is challenging to partition the evaporation from rainfall or snowmelt water sources.
Another potential source of error in estimating the rainfall/snowmelt contribution towards winter drainage was the definition of the winter period. For the winter period, it was defined as the first day in November towards the last day of April. Such a definition can be found in multiple winter hydrology research for Canada, as the snow could fall as early as October and as late as June in Canada (Newton et al., 2014; Van der Kamp et al., 2003; Zhang et al., 2001). However, heavy rainfall could continue to occur after the snow has been completely melted for the two sites (83 mm of rainfall in April for the St. Emmanuel site and 73 mm of rainfall in April for the Harrow site) in the current study. As a result, the contribution of rainfall towards winter drainage could be overestimated; therefore, the snowmelt contribution can be otherwise underestimated.
Conclusion
In the current study, winter subsurface drainage discharge data from two winterized drainage experimental sites in eastern Canada, combined with the RZ-SHAW model, was used to evaluate the RZ-SHAW model and assess the 30-year (1990-2022) contribution of snowmelt to subsurface drainage discharge under a cold climate. RZ-SHAW generally delivered satisfactory results in simulating winter drainage discharge, soil temperature, and snow depth during the model evaluation (2021 to 2023 for the St. Emmanuel site and 2000-2003 for the Harrow site). The prediction of 30-year snowmelt and rainfall's contribution towards the subsurface drainage discharge during the winter and spring snowmelt seasons were successfully predicted based on the calibrated and validated RZ-SHAW scenarios. From 1990 to 2022, snowmelt was either the dominant contributor (St. Emmanuel site, 69%) or contributed significant drainage discharge (Harrow site, 48%) for the spring snowmelt season drainage at both sites. For the winter season drainage, snowmelt was still the dominant contributor in the St. Emmanuel site (55%), while its contribution in the Harrow site was 29%. The snowmelt drainage contributed 44% and 9% of the annual drainage for the St. Emmanuel and Harrow sites, respectively. Overall, snowmelt contributes significantly towards not only winter but also annual drainage. It was discovered that RZ-SHAW tends to simulate early drainage during spring while underestimating peak drainage during spring compared to the observed values. Considering less precipitation has been projected to fall as snow under the warming climate in the temperate regions of Canada, increased attention should be given to adapting and developing new winter management for crop planting windows and nutrient loss in future research.
Acknowledgments
We acknowledge the Natural Sciences and Engineering Research Council of Canada (NSERC) for sponsoring this research. We thank the Chinese Scholarship Council and the Angus F. MacKenzie Graduate Fellowship in Sustainable Agriculture for providing the scholarship and financial support for the doctoral student.
Adhikari, B. K., Madramootoo, C. A., & Sarangi, A. (2010). Temporal variability of phosphorus flux from Pike River watershed to the Missisquoi Bay of Quebec. Curr. Sci., 98(1), 58-64. Retrieved from https://www.jstor.org/stable/24111610
Ahmed, I., Rudra, R., McKague, K., Gharabaghi, B., & Ogilvie, J. (2007). Evaluation of the Root Zone Water Quality Model (RZWQM) for southern Ontario: Part I. Sensitivity analysis, calibration, and validation. Water Qual. Res. J., 42(3), 202-218. doi:https://doi.org/10.2166/wqrj.2007.024
Ahmed, S. I., Rudra, R., Goel, P., Amili, A., Dickinson, T., Singh, K., & Khan, A. (2022). Change in winter precipitation regime across Ontario, Canada. Hydrology, 9(5), 81. doi:https://doi.org/10.3390/hydrology9050081
Ahuja, L., Rojas, K. W., Hanson, J. D., Shaffer, M. J., & Ma, L. (2000). Root zone water quality model: Modelling management effects on water quality and crop production. Colorado: Water Resources Publication.
Anderson, E. A. (1976). A point energy and mass balance model of a snow cover. 19. NOAA technical report NWS. Retrieved from https://repository.library.noaa.gov/view/noaa/6392
Askar, M. H., Youssef, M. A., Chescheir, G. M., Negm, L. M., King, K. W., Hesterberg, D. L., . . . Skaggs, R. W. (2020). DRAINMOD Simulation of macropore flow at subsurface drained agricultural fields: Model modification and field testing. Agric. Water Manag., 242, 106401. doi:https://doi.org/10.1016/j.agwat.2020.106401
Boico, V. F., Therrien, R., Delottier, H., Young, N. L., & Højberg, A. L. (2022). Comparing alternative conceptual models for tile drains and soil heterogeneity for the simulation of tile drainage in agricultural catchments. J. Hydrol., 612, 128120. doi:https://doi.org/10.1016/j.jhydrol.2022.128120
Bouchard, E., & Qi, Z. (2016). Long-term trends of climate change and its impact on crop growing season on Montreal Island. J. Water Clim. Change, 8(1), 78-88. doi:https://doi.org/10.2166/wcc.2016.139
Chanasyk, D. S., & Woytowich., C. P. (1987). Sediment yield as a result of snowmelt runoff in the Peace River region. Can. Agric. Eng., 29(1), 1-6. Retrieved from https://library.csbe-scgab.ca/all-publications/761:sediment-yield-as-a-result-of-snowmelt-runoff-in-the-peace-river-region
Coelho, B. B., Murray, R., Lapen, D., Topp, E., & Bruin, A. (2012). Phosphorus and sediment loading to surface waters from liquid swine manure application under different drainage and tillage practices. Agric. Water Manag., 104, 51-61. doi:https://doi.org/10.1016/j.agwat.2011.10.020
Crézé, C. (2015). Greenhouse gas emissions from an intensively cropped field under various water and fertilizer management practices. MS thesis. Montreal, Quebec: McGill University, Department of Bioresource Engineering.
Ding, R., Kang, S., Zhang, Y., Hao, X., Tong, L., & Du, T. (2013). Partitioning evapotranspiration into soil evaporation and transpiration using a modified dual crop coefficient model in irrigated maize field with ground-mulching. Agric. Water Manag., 127, 85-96. doi:https://doi.org/10.1016/j.agwat.2013.05.018
Eastman, M., Gollamudi, A., Stämpfli, N., Madramootoo, C. A., & Sarangi, A. (2010). Comparative evaluation of phosphorus losses from subsurface and naturally drained agricultural fields in the Pike River watershed of Quebec, Canada. Agric. Water Manag., 97(5), 596-604. doi:https://doi.org/10.1016/j.agwat.2009.11.010
Enright, P., & Madramootoo, C. A. (2004). Phosphorus losses in surface runoff and subsurface drainage waters on two agricultural fields in Quebec. Drainage VIII Proc. 8th Int. Symp. St. Joseph, MI: ASABE. doi:https://doi.org/10.13031/2013.15722
Environment and Climate Change Canda (ECCC). (2018a). Canadian historical snow survey data. Retrieved from https://data-donnees.ec.gc.ca/data/climate/systems/canadian-historical-snow-survey-data
Environment and Climate Change Canda (ECCC). (2018b). Weather tools: Frequently asked questions. Retrieved from https://www.canada.ca/en/environment-climate-change/services/weather-general-tools-resources/frequently-asked-questions.html
Environment and Climate Change Canda (ECCC). (2020). Temperature change in Canada. Retrieved from https://www.canada.ca/en/environment-climate-change/services/climate-change/canadian-centre-climate-services/basics/trends-projections/changes-snow.html
Environment and Climate Change Canda (ECCC). (2021). Historical Climate Data [Dataset]. Retrieved from https://climate.weather.gc.ca/historical_data/search_historic_data_e.html
Flerchinger, G. N., Aiken, R. M., Rojas, K. W., & Ahuja, L. R. (2000). Development of the Root Zone Water Quality Model (RZWQM) for over-winter conditions. Trans. ASAE, 43(1), 59-68. doi:https://doi.org/10.13031/2013.2688
Flerchinger, G. N., Caldwell, T. G., Cho, J., & Hardegree, S. P. (2012). Simultaneous Heat and Water (SHAW) Model: Model use, calibration, and validation. Trans. ASABE, 55(4), 1395-1411. doi:https://doi.org/10.13031/2013.42250
Flerchinger, G., Aiken, R., Rojas, K., Ahuja, L., Johnsen, K., Alonso, C., & Hanson, J. D. (1999). Soil heat transport, soil freezing and snowpack conditions. In Root Zone Water Quality Model: Modelling management effects on water quality and crop production (pp. 281-314). Colorado: Water Resources Publications.
Gelaro, R., McCarty, W., Suárez, M. J., Todling, R., Molod, A., Takacs, L., . . . Zhao, B. (2017). The modern-era retrospective analysis for research and applications, Version 2 (MERRA-2). J. Climate, 30(14), 5419-5454. doi:https://doi.org/10.1175/JCLI-D-16-0758.1
Ghane, E. (2022). Choice of pipe material influences drain spacing and system cost in subsurface drainage design. Appl. Eng. Agric., 38(4), 685-695. doi:https://doi.org/10.13031/aea.15053
Gupta, H. V., Kling, H., Yilmaz, K. K., & Martinez, G. F. (2009). Decomposition of the mean squared error and NSE performance criteria: Implications for improving hydrological modelling. J. Hydrol., 377(1), 80-91. doi:https://doi.org/10.1016/j.jhydrol.2009.08.003
Hansson, K., Šimunek, J., Mizoguchi, M., Lundin, L.-C., & van Genuchten, M. T. (2004). Water flow and heat transport in frozen soil: Numerical solution and freeze-thaw applications. Vadose Zone J., 3(2), 693-704. doi:https://doi.org/10.2136/vzj2004.0693
He, J., Diao, Z., Zheng, Z., Su, D., & Lyu, S. (2020). Laboratory investigation of phosphorus loss with snowmelt and rainfall runoff from a Steppe wetland catchment. Chemosphere, 241, 125137. doi:https://doi.org/10.1016/j.chemosphere.2019.125137
Henry, H. A. (2007). Soil freeze-thaw cycle experiments: Trends, methodological weaknesses and suggested improvements. Soil Biol. Biochem., 39(5), 977-986. doi:https://doi.org/10.1016/j.soilbio.2006.11.017
Hoffman, A. R., Polebitski, A. S., Penn, M. R., & Busch, D. L. (2019). Long-term variation in agricultural edge-of-field phosphorus transport during snowmelt, rain, and mixed runoff events. J. Environ. Qual., 48(4), 931-940. doi:https://doi.org/10.2134/jeq2018.11.0420
Jamieson, A. (2001). Evaluating phosphorus losses in surface and subsurface runoff from two agricultural fields in Quebec. MS thesis. Montreal, Quebec: McGill University, Department of Bioresource Engineering.
Jamieson, A., Madramootoo, C. A., & Enright, P. (2003). Phosphorus losses in surface and subsurface runoff from a snowmelt event on an agricultural field in Quebec. Can. Biosyst. Eng., 45. Retrieved from https://library.csbe-scgab.ca/docs/journal/45/c0208.pdf
Jeantet, A., Henine, H., Chaumont, C., Collet, L., Thirel, G., & Tournebize, J. (2021). Robustness of a parsimonious subsurface drainage model at the French national scale. Hydrol. Earth Syst. Sci., 25(10), 5447-5471. doi:https://doi.org/10.5194/hess-25-5447-2021
Jégo, G., Chantigny, M., Pattey, E., Bélanger, G., Rochette, P., Vanasse, A., & Goyer, C. (2014). Improved snow-cover model for multi-annual simulations with the STICS crop model under cold, humid continental climates. Agric. For. Meteorol., 195-196, 38-51. doi:https://doi.org/10.1016/j.agrformet.2014.05.002
Jennings, K. S., Winchell, T. S., Livneh, B., & Molotch, N. P. (2018). Spatial variation of the rain-snow temperature threshold across the Northern Hemisphere. Nat. Commun., 9(1), 1148. doi:https://doi.org/10.1038/s41467-018-03629-7
Jiang, Q., Qi, Z., Lu, C., Tan, C. S., Zhang, T., & Prasher, S. O. (2020). Evaluating RZ-SHAW model for simulating surface runoff and subsurface tile drainage under regular and controlled drainage with subirrigation in southern Ontario. Agric. Water Manag., 237, 106179. doi:https://doi.org/10.1016/j.agwat.2020.106179
Jiang, Q., Qi, Z., Madramootoo, C. A., & Crézé, C. (2019). Mitigating greenhouse gas emissions in subsurface-drained field using RZWQM2. Sci. Total Environ., 646, 377-389. doi:https://doi.org/10.1016/j.scitotenv.2018.07.285
Jiang, Q., Qi, Z., Madramootoo, C. A., & Singh, A. K. (2018). Simulating hydrologic cycle and crop production in a subsurface drained and sub-irrigated field in Southern Quebec using RZWQM2. Comput. Electron. Agric., 146, 31-42. doi:https://doi.org/10.1016/j.compag.2018.01.021
Kaluli, J. W., Madramootoo, C. A., Zhou, X., MacKenzie, A. F., & Smith, D. L. (1999). Subirrigation systems to minimize nitrate leaching. J. Irrig. Drain. Eng., 125(2), 52-58. doi:https://doi.org/10.1061/(ASCE)0733-9437(1999)125:2(52)
King, K. W., Williams, M. R., Macrae, M. L., Fausey, N. R., Frankenberger, J., Smith, D. R., . . . Brown, L. C. (2015). Phosphorus transport in agricultural subsurface drainage: A review. J. Environ. Qual., 44(2), 467-485. doi:https://doi.org/10.2134/jeq2014.04.0163
Kokulan, V., Macrae, M. L., Lobb, D. A., & Ali, G. A. (2019). Contribution of overland and tile flow to runoff and nutrient losses from vertisols in Manitoba, Canada. J. Environ. Qual., 48(4), 959-965. doi:https://doi.org/10.2134/jeq2019.03.0103
Lam, W. V., Macrae, M. L., English, M. C., O'Halloran, I. P., Plach, J. M., & Wang, Y. (2016). Seasonal and event-based drivers of runoff and phosphorus export through agricultural tile drains under sandy loam soil in a cool temperate region. Hydrol. Process., 30(15), 2644-2656. doi:https://doi.org/10.1002/hyp.10871
Leavesley, G. H., Lichty, R. W., Troutman, B. M., & Saindon, L. G. (1983). Precipitation-runoff modeling system; user's manual. Water-Resources Investigations Report 83-4238. USGS, Water Resources Division. doi:https://doi.org/10.3133/wri834238
Li, D., Lettenmaier, D. P., Margulis, S. A., & Andreadis, K. (2019). The role of rain-on-snow in flooding over the conterminous United States. Water Resour. Res., 55(11), 8492-8513. doi:https://doi.org/10.1029/2019WR024950
Li, Z., Ma, L., Flerchinger, G. N., Ahuja, L. R., Wang, H., & Li, Z. (2012). Simulation of overwinter soil water and soil temperature with SHAW and RZ-SHAW. Soil Sci. Soc. Am. J., 76(5), 1548-1563. doi:https://doi.org/10.2136/sssaj2011.0434
Li, Z., Qi, Z., Smith, W., Pattey, E., & Qian, B. (2022). Long-term simulation of snow cover and its potential impacts on seasonal frost dynamics in croplands across Southern Canada. Water Resour. Res., 58(8), e2021WR031674. doi:https://doi.org/10.1029/2021WR031674
Luo, W., Skaggs, R. W., & Chescheir, G. M. (2000). Drainmod modifications for cold conditions. Trans. ASAE, 43(6), 1569-1582. doi:https://doi.org/10.13031/2013.3057
Ma, L., Ahuja, L. R., Nolan, B. T., Malone, R. W., Trout, T. J., & Qi, Z. (2012). Root Zone Water Quality Model (RZWQM2): Model use, calibration, and validation. Trans. ASABE, 55(4), 1425-1446. doi:https://doi.org/10.13031/2013.42252
Ma, L., Ahuja, L. R., Saseendran, S. A., Malone, R. W., Green, T. R., Nolan, B. T., . . . Hoogenboom, G. (2011). A protocol for parameterization and calibration of RZWQM2 in field research. In L. R. Ahuja, & L. Ma (Eds.), Methods of introducing system models into agricultural research, Volume 2 (pp. 1-64). American Society of Agronomy. doi:https://doi.org/10.2134/advagricsystmodel2.c1
Ma, L., Hoogenboom, G., Ahuja, L. R., Nielsen, D. C., & Ascough II, J. C. (2005). Development and evaluation of the RZWQM-CROPGRO hybrid model for soybean production. Agron. J., 97(4), 1172-1182. doi:https://doi.org/10.2134/agronj2003.0314
Ma, L., Hoogenboom, G., Saseendran, S. A., Bartling, P. N., Ahuja, L. R., & Green, T. R. (2009). Effects of estimating soil hydraulic properties and root growth factor on soil water balance and crop production. Agron. J., 101(3), 572-583. doi:https://doi.org/10.2134/agronj2008.0206x
Ma, Q. L., Hook, J. E., & Wauchope, R. D. (1999). Evapotranspiration predictions: A comparison among GLEAMS, Opus, PRZM-2, and RZWQM models in a humid and thermic climate. Agric. Syst., 59(1), 41-55. doi:https://doi.org/10.1016/S0308-521X(98)00081-X
Madramootoo, C. A., Broughton, R. S., & Tait, R. K. (1994). Technical notes: Trenchless installation of vertical plastic curtains. Trans. ASAE, 37(5), 1525-1527. doi:https://doi.org/10.13031/2013.28236
Madramootoo, C. A., Helwig, T. G., & Dodds, G. T. (2001). Managing water tables to improve drainage water quality in Quebec, Canada. Trans. ASAE, 44(6), 1511. doi:https://doi.org/10.13031/2013.7034
Moriasi, D. N., Gitau, M. W., Pai, N., & Daggupati, P. (2015). Hydrologic and water quality models: Performance measures and evaluation criteria. Trans. ASABE, 58(6), 1763-1785. doi:https://doi.org/10.13031/trans.58.10715
Morrison, J., Madramootoo, C. A., & Chikhaoui, M. (2014). Modeling agricultural land drainage under spring snowmelt conditions with DRAINMOD. Can. J. Civ. Eng., 41(4), 275-284. doi:https://doi.org/10.1139/cjce-2013-0416
Newton, B. W., Prowse, T. D., & Bonsal, B. R. (2014). Evaluating the distribution of water resources in western Canada using synoptic climatology and selected teleconnections. Part 1: Winter season. Hydrol. Process., 28(14), 4219-4234. doi:https://doi.org/10.1002/hyp.10233
Panuska, J. C., Karthikeyan, K. G., & Norman, J. M. (2008). Sediment and phosphorus losses in snowmelt and rainfall runoff from three corn management systems. Trans. ASABE, 51(1), 95-105. doi:https://doi.org/10.13031/2013.24230
Peel, M. C., Finlayson, B. L., & McMahon, T. A. (2007). Updated world map of the Köppen-Geiger climate classification. Hydrol. Earth Syst. Sci., 11(5), 1633-1644. doi:https://doi.org/10.5194/hess-11-1633-2007
Plach, J., Pluer, W., Macrae, M., Kompanizare, M., McKague, K., Carlow, R., & Brunke, R. (2019). Agricultural edge-of-field phosphorus losses in Ontario, Canada: Importance of the nongrowing season in cold regions. J. Environ. Qual., 48(4), 813-821. doi:https://doi.org/10.2134/jeq2018.11.0418
Pomeroy, J. W., Fang, X., & Marks, D. G. (2016). The cold rain-on-snow event of June 2013 in the Canadian Rockies — characteristics and diagnosis. Hydrol. Process., 30(17), 2899-2914. doi:https://doi.org/10.1002/hyp.10905
PyTorch. (2019). SOBOLENGINE. Retrieved from https://pytorch.org/docs/stable/generated/torch.quasirandom.SobolEngine.html
Qi, H., & Qi, Z. (2017). Simulating phosphorus loss to subsurface tile drainage flow: A review. Environ. Rev., 25(2), 150-162. doi:https://doi.org/10.1139/er-2016-0024
Qi, J., Li, S., Li, Q., Xing, Z., Bourque, C. P., & Meng, F.-R. (2016). A new soil-temperature module for SWAT application in regions with seasonal snow cover. J. Hydrol., 538, 863-877. doi:https://doi.org/10.1016/j.jhydrol.2016.05.003
Qian, B., Gregorich, E. G., Gameda, S., Hopkins, D. W., & Wang, X. L. (2011). Observed soil temperature trends associated with climate change in Canada. J. Geophys. Res.: Atmos., 116(D2). doi:https://doi.org/10.1029/2010JD015012
Ren, D., Engel, B., Mercado, J. A., Guo, T., Liu, Y., & Huang, G. (2022). Modeling and assessing water and nutrient balances in a tile-drained agricultural watershed in the U.S. Corn Belt. Water Res., 210, 117976. doi:https://doi.org/10.1016/j.watres.2021.117976
Ruggiero, R., Ross, D., & Faulkner, J. W. (2022). Tile drainage flow partitioning and phosphorus export in Vermont USA. Agriculture, 12(2), 167. doi:https://doi.org/10.3390/agriculture12020167
Sadhukhan, D., Qi, Z., Zhang, T., Tan, C. S., Ma, L., & Andales, A. A. (2019). Development and evaluation of a phosphorus (P) module in RZWQM2 for phosphorus management in agricultural fields. Environ. Model. Softw., 113, 48-58. doi:https://doi.org/10.1016/j.envsoft.2018.12.007
Shokrana, M. S., Ghane, E., & Qi, Z. (2023). Calibration and validation of RZWQM2-P model to simulate phosphorus loss in a clay loam soil in Michigan. J. ASABE, 66(1), 1-12. doi:https://doi.org/10.13031/ja.15283
Singh, M., Prasher, S. O., Tan, C. S., & Tejawat, C. M. (1994). Evaluation of drainmod for southern Ontario conditions. Can. Water Resour. J., 19(4), 313-326. doi:https://doi.org/10.4296/cwrj1904313
Skaggs, R. W. (1978). A water management model for shallow water table soils. Report No. 134. Raleigh, NC: University of North Carolina, Water Resources Research Institute. Retrieved from https://repository.lib.ncsu.edu/handle/1840.4/1618
Sobol, I. M. (1967). On the distribution of points in a cube and the approximate evaluation of integrals. USSR Comput. Math. Math. Phys., 7(4), 86-112. doi:https://doi.org/10.1016/0041-5553(67)90144-9
Spomer, R. G., & Hjelmfelt Jr., A. T. (1983). Snowmelt runoff and erosion on Iowa loess soils. Trans. ASAE, 26(4), 1109-1111. doi:https://doi.org/10.13031/2013.34086
Stämpfli, N., & Madramootoo, C. A. (2006). Dissolved phosphorus losses in tile drainage under subirrigation. Water Qual. Res. J., 41(1), 63-71. doi:https://doi.org/10.2166/wqrj.2006.007
Su, J. J., van Bochove, E., Thériault, G., Novotna, B., Khaldoune, J., Denault, J. T., . . . Chow, L. (2011). Effects of snowmelt on phosphorus and sediment losses from agricultural watersheds in Eastern Canada. Agric. Water Manag., 98(5), 867-876. doi:https://doi.org/10.1016/j.agwat.2010.12.013
Tait, R., Madramootoo, C. A., & Enright, P. (1995). An instrumented, field-scale research facility for drainage and water quality studies. Comput. Electron. Agric., 12(2), 131-145. doi:https://doi.org/10.1016/0168-1699(94)00043-P
Tan, C. S., & Zhang, T. Q. (2011). Surface runoff and sub-surface drainage phosphorus losses under regular free drainage and controlled drainage with sub-irrigation systems in southern Ontario. Can. J. Soil Sci., 91(3), 349-359. doi:https://doi.org/10.4141/cjss09086
Tan, C. S., Drury, C. F., Gaynor, J. D., & Welacky, T. W. (1993). Integrated soil, crop and water management system to abate herbicide and nitrate contamination of the great lakes. Water Sci. Technol., 28(3-5), 497-507. doi:https://doi.org/10.2166/wst.1993.0453
Valero, C. S., Madramootoo, C. A., & Stämpfli, N. (2007). Water table management impacts on phosphorus loads in tile drainage. Agric. Water Manag., 89(1), 71-80. doi:https://doi.org/10.1016/j.agwat.2006.12.007
van der Kamp, G., Hayashi, M., & Gallén, D. (2003). Comparing the hydrology of grassed and cultivated catchments in the semi-arid Canadian prairies. Hydrol. Process., 17(3), 559-575. doi:https://doi.org/10.1002/hyp.1157
Van Esbroeck, C. J., Macrae, M. L., Brunke, R. I., & McKague, K. (2016). Annual and seasonal phosphorus export in surface runoff and tile drainage from agricultural fields with cold temperate climates. J. Great Lakes Res., 42(6), 1271-1280. doi:https://doi.org/10.1016/j.jglr.2015.12.014
Wang, X. L., Xu, H., Qian, B., Feng, Y., & Mekis, E. (2017). Adjusted daily rainfall and snowfall data for Canada. Atmos. Ocean, 55(3), 155-168. doi:https://doi.org/10.1080/07055900.2017.1342163
Wang, Z., Qi, Z., Xue, L., & Bukovsky, M. (2016). RZWQM2 simulated management practices to mitigate climate change impacts on nitrogen losses and corn production. Environ. Model. Softw., 84, 99-111. doi:https://doi.org/10.1016/j.envsoft.2016.06.016
Wang, Z., Zhang, T. Q., Tan, C. S., Wang, X., Taylor, R. A., Qi, Z. M., & Yang, J. W. (2019). Modeling the impacts of manure on phosphorus loss in surface runoff and subsurface drainage. J. Environ. Qual., 48(1), 39-46. doi:https://doi.org/10.2134/jeq2018.06.0240
Wever, N., Würzer, S., Fierz, C., & Lehning, M. (2016). Simulating ice layer formation under the presence of preferential flow in layered snowpacks. Cryosphere, 10(6), 2731-2744. doi:https://doi.org/10.5194/tc-10-2731-2016
Whiteley, H. (2004). Influence of method of measurement of daily snowfall on climate normals in Ontario, Canada. Proc. 61st Eastern Snow Conf. Retrieved from https://static1.squarespace.com/static/58b98f7bd1758e4cc271d365/t/5e617e05f8a4f25611953240/1583447557995/05+Whiteley.pdf
Yang, J., Zhou, W., Liu, J., & Hu, X. (2014). Dynamics of greenhouse gas formation in relation to freeze/thaw soil depth in a flooded peat marsh of Northeast China. Soil Biol. Biochem., 75, 202-210. doi:https://doi.org/10.1016/j.soilbio.2014.04.006
Zhang, T. Q., Tan, C. S., Zheng, Z. M., Welacky, T. W., & Reynolds, W. D. (2015). Impacts of soil conditioners and water table management on phosphorus loss in tile drainage from a clay loam soil. J. Environ. Qual., 44(2), 572-584. doi:https://doi.org/10.2134/jeq2014.04.0154
Zhang, T. Q., Tan, C. S., Zheng, Z. M., Welacky, T., & Wang, Y. T. (2017). Drainage water management combined with cover crop enhances reduction of soil phosphorus loss. Sci. Total Environ., 586, 362-371. doi:https://doi.org/10.1016/j.scitotenv.2017.02.025
Zhang, X., Harvey, K. D., Hogg, W. D., & Yuzyk, T. R. (2001). Trends in Canadian streamflow. Water Resour. Res., 37(4), 987-998. doi:https://doi.org/10.1029/2000WR900357
Zhang, X., Vincent, L. A., Hogg, W. D., & Niitsoo, A. (2000). Temperature and precipitation trends in Canada during the 20th century. Atmos. Ocean, 38(3), 395-429. doi:https://doi.org/10.1080/07055900.2000.9649654