Effects of climate change on depression‐focused groundwater recharge in the Canadian Prairies

Small topographic depressions are ubiquitous in the glaciated terrain of the Northern Prairies characterized by a cold semiarid climate. Groundwater recharge in this region is focused in topographic depressions, which receive lateral inputs of snowmelt runoff in addition to vertical inputs of precipitation. The response of depression‐focused recharge to climate change is largely unknown due to the difficulty of representing the complex interaction between depressions and surrounding uplands in hydrological models. This study evaluates climate change effects on recharge using a process‐based hydrological model and the pseudo‐global warming (PGW) approach representing a climate of 2092–2100, which has a higher mean annual temperature (+4.9 °C) and precipitation (+135 mm or +27%) compared with the present climate. The recharge model is conditioned using field data from an instrumented grassland site in the Canadian Prairies. Under the present climate, snowmelt runoff occurred in March–April over frozen soil, which accounted for the majority of water transfer from the upland to the depression. Under the PGW scenario, the amount of snowmelt runoff in March–April decreased due to lower snow accumulation and shorter periods of frozen soil. Instead, runoff occurred in midwinter and also in summer months due to increased intensity and duration of summer storms. Despite the increased precipitation, mean annual recharge rates decreased from 10.2 to 3.2 mm yr–1, indicating the importance of considering the complex effects of warmer and wetter climate on hydrological processes in cold semiarid regions.

management (Pierce et al., 2013;Zhou, 2009), but its response to climate change has a large degree of uncertainty (Green et al., 2011;Holman, 2006). In cold regions of Eurasia and North America, where snow and seasonally frozen soil have a major influence on hydrological processes, climate change will cause earlier snowmelt and delayed snow accumulation (Jungqvist et al., 2014;Rasouli et al., 2014), shorter frost periods (Henry, 2008), and changes in snowmelt rates (Musselman et al., 2017). All of these could affect snowmelt runoff and infiltration (Sutinen et al., 2008), and subsequently groundwater recharge. Previous studies on climate change impacts on groundwater recharge were mostly conducted in temperate regions (Crosbie et al., 2013;Eckhardt & Ulbrich, 2003;Jyrkama & Sykes, 2007), with a few exceptions in cold regions-for example, in Finland (Okkonen et al., 2010) and Canada (Zhang et al., 2020). Hence, there is a need to advance our understanding of climate change effects on groundwater recharge in cold regions.
The spatial distribution of recharge is largely dependent on the balance between infiltration and evapotranspiration. In humid regions, where precipitation exceeds potential evaporation, recharge can occur over wide areas within the landscape, and it is considered "direct" (de Vries & Simmers, 2002). In contrast, in arid and semiarid regions, where potential evaporation exceeds precipitation, recharge is restricted to localized areas. In these areas, lateral transfer of surface water increases infiltration in localized areas such as ephemeral ponds (Greenwood & Buttle, 2018;Scanlon & Goldsmith, 1997) and intermittent streams (Izbicki et al, 2000). This mode of recharge is called "localized" or "focused" (Scanlon et al., 2006). Focused recharge may occur in humid regions, but it is most dominant in arid and semiarid regions.
The Northern Prairies or Prairie Pothole Region of North America ( Figure 1a) is a major physiographic region characterized by a cold semiarid climate and glaciated landscape containing numerous topographic depressions (see Section 2). Annual potential evaporation greatly exceeds annual precipitation in this region, and much of the precipitation inputs are consumed by evapotranspiration during the growing season. Therefore, groundwater recharge is small and focused in small topographic depressions (10 2 -10 3 m 2 ) that receive lateral inputs of snowmelt runoff generated over frozen ground (Bam et al., 2020;Berthold et al., 2004;Hayashi et al., 1998).
The depressions in the Northern Prairies have a density of a few to few tens per square kilometer (e.g., Figure 1b), and the size and density of depressions have a strong influence on the degree to which the water input is focused (Noorduijn et al., 2018). Therefore, for reliable assessment of climate change impacts on groundwater recharge in the Northern Prairies, it is necessary to understand cold-region hydrological processes at a scale of individual depressions and, at the same time, develop a numerically efficient model that can be used to simulate these processes for a large number (10 3 -10 4 ) of depressions distributed over a large region. This is because models using diffuse recharge concept and large grid cells cannot adequately represent depression-focused recharge (Zhang et al., 2020), which is the essential feature of the Northern Prairies region.
A wide range of models have been used to estimate groundwater recharge, from simple empirical models (Chen et al., 2002) to physically based models with various complexity levels (Loukili et al., 2008). Soil water balance models have simple algorithms for fast computation, and an intermediate level of complexity to allow for representation of essen-

Core Ideas
• Cold semiarid climate and glaciated landscape focus infiltration in depressions. • Depression-focused recharge is driven by snowmelt runoff and summer rain. • Groundwater recharge will decrease as the climate becomes warmer and wetter.
tial hydrological processes such as snow accumulation and melt, overland flow, infiltration, evapotranspiration, and soil moisture redistribution and drainage. They are widely used to estimate spatially variable, transient groundwater recharge fluxes (Jyrkama & Sykes, 2007;Scibek & Allen, 2006). In this study, a process-based soil water balance model is used to simulate depression-focused groundwater recharge involving snowmelt runoff over frozen soil (see Section 3). A numerical weather model is used to downscale dynamically the outputs of general circulation models and provide a climate change scenario to the groundwater recharge model. The objective of this study is to evaluate potential effects of climate change on hydrological processes in depressions and their catchments, and to examine the interaction and feedback between the processes and their overall influences on depression-focused recharge. To achieve these goals, the longterm (2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015) data collected at an instrumented study site in the Canadian Prairies are used to test a newly developed recharge model, downscaled climate model outputs are corrected for biases using the current climate data, and the biascorrected climate change scenario for 2092-2100 is used to simulate recharge processes.

STUDY SITE AND FIELD METHODS
The Canadian Prairies are part of the Northern Prairies, which cover more than 750,000 km 2 of North America (Figure 1a). They are characterized by numerous topographic depressions of varying sizes, which play an important role in hydrological process (Hayashi et al., 2016). The data used in this study were collected at a grassland site on the Spyhill Farm, near Calgary, AB, Canada (51˚10′31" N, 114˚13′44″ W). Located in the western edge of the Canadian prairies, the site has a cold semiarid climate with frequent occurrence of midwinter melt events (Pavlovskii, Hayashi, & Itenfisu, 2019 (Zaitlin et al., 2007). The site was previously used as a summer pasture for cattle grazing but has not been grazed since 2006. This study mainly uses the data collected from depression GP, and supplemental data from depressions N12, G16, and G17 ( Figure 1b). The water table under the upland surrounding GP is approximately 9 m below the ground surface. The water table under GP is approximately 3.5 m below the ground surface during dry periods and rises to the surface when the depression is flooded .
An automated weather station was installed in an upland adjacent to depression GP ( Figure 1b) in June 2006 to monitor meteorological and soil variables on a 30-min interval, including air temperature and relative humidity (Vaisala, HMP45C), wind speed (RM Young, 05103), long-and shortwave radiation (Kip & Zonnen, CNR1), precipitation (Geonor, T200B), turbulent fluxes (Campbell Scientific eddy-covariance system consisting of CSAT3 and KH20), soil water content at depths of 0.1, 0.3, 0.6, and 1 m (Campbell Scientific, CS616), soil temperature at depths of 0, 0.2, 0.4, 0.6, 0.8, 1.0, and 1.5 m (thermocouples), and ground heat flux at a depth of 0.05 m (Campbell Scientific, HFT3). Snow surveys were conducted on approximately biweekly to monthly intervals on a 100-m survey line to monitor snow water equivalent, and the amount of snowmelt runoff was estimated from the volume of water collected in depressions. This method cannot account for the initial loss of runoff water to air-filled pore spaces in the topsoil of the depression during the early stage of snowmelt runoff and hence tends to underestimate runoff (Mohammed et al., 2013). However, it is more reliable than point estimates of runoff using runoff traps , as it integrates local-scale variability over the entire uplands. The water level in depressions was monitored manually or by using pressure transducers (in situ, Mini-Troll; Solinst, Levelogger). Details on sensor heights, calibration, and data processing procedures are found in Mohammed et al. (2013).

Groundwater recharge model
Groundwater recharge in the Canadian Prairies is strongly focused in numerous depressions distributed over the landscape ( Figure 2a). Among possible modeling approaches are integrated hydrological models solving the three-dimensional (3D) Richards equation at the most complex end of the spectrum and one-dimensional (1D) soil water balance models at the least complex end (Figure 2b). The former has an advantage of representing all hydrological processes in the most rigorous manner. Examples include HydroGeoSphere (Shilling et al., 2019), GEOTop (Endrizzi et al., 2014), and Advanced Terrestrial Simulator (Jan et al., 2020). The latter has an advantage of fast computational time while capturing essential processes. This study uses coupled 1D models consisting of upland and depression ( Figure 2c) because of the intended model application involving 10 3 s of depressions to obtain regional-scale recharge estimates (see Section 5.2). The model has been developed from the Versatile Soil Moisture Budget (VSMB) model, which is widely used to simulate soil water balance in the Canadian Prairies for agricultural applications (Government of Alberta, 2020). The original VSMB was developed by Baier and Robertson (1966) to estimate crop-available water from daily air temperature and precipitation data using the water balance of multiple soil layers. The model has undergone a series of improvements by Baier et al. (1972Baier et al. ( , 1979, Dyer and Mack (1984), Akinremi et al. (1996), Hayashi et al. (2010), and Mohammed et al. (2013) to represent relevant hydrological processes in a more physically based manner, including an improved representation of snow and frozen soil processes and macropore infiltration in frozen soil.
The VSMB in this study calculates daily soil water balance for seven soil layers (five for depression) ( Open arrows indicate gravitational drainage, and thin solid arrows indicate moisture diffusion. The drainage at the bottom of the soil profiles gives recharge from upland (R u ) and depression (R d ). (d) Flow chart of simulation steps using Utah Energy Balance (UEB) and VSMB models. 1D, one-dimensional; 3D, three-dimensional 0.8, 1.2, and 2.0 m. The model is forced by meteorological data consisting of hourly air temperature, relative humidity, precipitation, wind speed, and shortwave incoming radiation. During the winter season (1 November-30 April), the data are input to the Utah Energy Balance (UEB) model (Tarboton & Luce, 1996) to simulate snowpack evolution and calculate daily snowmelt input to the two VSMB models representing upland and depression ( Figure 2d). In nonwinter seasons, hourly data are used to calculate daily values of meteorological variables and directly input to the VSMB models. Rainfall is applied the soil surface layer after subtracting canopy interception. The plant canopy storage (i.e., maximum interception loss) is estimated to be 1.1 mm in this study based on a detailed laboratory and field study of mixed grass prairie in Saskatchewan by Couturier and Ripley (1973). Within the VSMB model, runoff is estimated using two different methods depending on frozen-unfrozen status of the top soil layer. When the top layer is unfrozen, runoff is estimated using the curve number method (USDA-NRCS, 2004) as described in Mohammed et al. (2013). When the top soil layer is frozen and unsaturated, snowmelt water cannot infiltrate into the next layer until the top layer reaches saturation (Mohammed et al., 2013). The liquid water infiltrates into the next layer if the added water exceeds the saturation limit; the infiltration rate in this process is limited by a userspecified constant (f lxm , m s −1 ) representing the effects of soil macropores (Mohammed et al., 2013). This sequence of processes is consistent with results of recent field observations (Mohammed et al., 2019) and laboratory experiments (Pittman et al., 2020). The remaining excess of water turns to generate runoff.
After the infiltration into the top soil layer, water is subsequently distributed to lower layers by first allowing the soil to drain to field capacity and then allowing water to move between layers by gradient-driven moisture diffusion. The evapotranspiration rate from each soil layer is individually calculated based on soil moisture, meteorological data, and plant growth stage (see Section 3.2). The coupled soil water model is called VSMB Depression-Upland System (VSMB-DUS) and uses a number of input parameters representing soil and plant properties as described in Section 3.3. The following describes those model algorithms that are particularly relevant to VSMB-DUS, whereas details of other process algorithms are found in Akinremi et al. (1996), Hayashi et al. (2010), Mohammed et al. (2013), and Noorduijn et al. (2018).
The VSMB-DUS transfers runoff from uplands to a depression as commonly observed in the Canadian Prairies (Noorduijn et al., 2018) (Figure 2c). Runoff from the upland accumulates in a pond when the volume of water inputs to the depression exceeds the infiltration capacity of the soil column (Noorduijn et al., 2018). The relationship between the volume (V p ), area (A p ), and depth (h p ) of the pond is given by (Hayashi & van der Kamp, 2000): (1) where s is a scaling constant (m 2 ), p is a dimensionless parameter representing the overall shape (i.e., degree of concavity) of the depression, and h 0 is a unit depth (=1 m). Table 1 lists the values of s and p determined by detailed elevation surveys of the four depressions, as well as their catchment areas (A c ) determined by detailed elevation survey for GP and from high-resolution (2 m) digital elevation maps for other depressions. Depressions in the Northern Prairies have a concentric distribution of plant species reflecting the soil moisture condition, which is wettest in the center (Stewart & Kantrud, 1972). The VSMB-DUS considers the effective area of depression (A d ) as the area occupied by those plant species that are adapted to the wet environment. Field observation at the four depressions indicates that the outer edge of these species roughly corresponds to elevation of 0.3-0.6 m above the lowest point in the depression ( Figure 1b). Therefore, A d is calibrated within the bounds defined by Equation 1 for h p = 0.3-0.6 m. The pond area can reach its maximum (A pmax ) when pond water rises to the spill point, and as a result, the excess water generates overflow (O) leaving the depression, although the pond depth in the four depressions never exceeded 0.5 m during 2003-2020. Losses from the ponded water can occur through direct evaporation from the water surface at the potential evaporation rate. Other losses can occur by infiltration to the underlying soil layer, including lateral flow to the unflooded area within the depression. The underlying soil layers may reach saturation under sustained infiltration. As a result, the saturated hydraulic conductivity (K s , m s −1 ) sets the upper limit of vertical flux between model layers. The amount of groundwater recharge is given by the drainage flux from the deepest soil layer. The bottom drainage flux is restricted by a model parameter (f bmax , m s −1 ) representing the combined effects of the hydraulic gradient and the hydraulic conductivity of clay-rich till underlying the soil column.
The drainage fluxes from the bottom of the two soil columns (R u and R d , Figure 2c) are weighted according to the relative areas of upland and depression to estimate catchmentaverage groundwater recharge ( Figure 2d). The lateral runoff from the upland to the depression provides a link between the upland and depression, and the area ratio of depression to upland, determines the degree to which atmospheric water inputs are focused in the depression ( Figure 2). As a result, the amount of depression-focused recharge is sensitive to the area ratio and the retention capacity of depressions, in addition to soil hydraulic parameters (Noorduijn et al., 2018).

Modification of potential evapotranspiration algorithm
Prior to this study, the VSMB model used the Priestly and Taylor (1972) equation to calculate potential evapotranspiration (PET) (Akinremi et al., 1996). Preliminary simulations showed that calculated PET was relatively insensitive to climate warming, as its temperature dependence was only reflected in changes in the slope of the saturation vapor pressure-temperature curve. Therefore, this study uses the Penman-Monteith equation in the form presented by Allen et al. (1998): is wind speed, and e s − e a (kPa) is vapor pressure deficit, all measured at 2m height. Note that the numbers appearing in Equation 3 are the result of unit conversion using specific units, as well as the following assumptions. It assumes that the aerodynamic resistance is equal to 208/u 2 s m −1 and the surface resistance is equal to 70 s m −1 representing a uniform and neutral boundary layer with fixed roughness lengths and a hypothetical crop height of 0.12 m (Allen et al., 1998, pp. 19-26).
The VSMB estimates actual evapotranspiration (E) by multiplying E p by empirical factors dependent on soil moisture and plant conditions: where r i is the root extraction coefficient for layer i, which is plant specific and linked to plant growth stages (Baier et al., 1979). The dimensionless drying curve function (f DC ) reduces evapotranspiration as the ratio of plant available water (S i ) to available water capacity (S Ci ) decreases and has the following form: where x = S i /S ci ; and C m , C n , C h , and C r are dimensionless fitting parameters. Note that S i and S Ci are defined as where θ i , θ WPi , and θ FCi are volumetric water content, wilting point, and field capacity, respectively, and Δz i is the thickness of soil layer i.

Model calibration
The soil and crop parameters in VSMB are calibrated within a plausible range encompassing the values reported at the Spyhill site by Mohammed et al. (2013), Noorduijn et al. (2018), and Mohammed et al. (2019). Calibrated parameters and their lower and upper bounds are listed in Table 2a (upland) and Table 2b (depression). In addition to the soil and crop parameters, several parameters related to snow processes are calibrated (z 0 , ρ s , D e , λ s , and r v ; see Table 2a for definition) and used for both upland and depression. These parameters have the strongest influence on snow accumulation and melt simulation based on the model sensitivity analysis by Mohammed et al. (2013). The model was calibrated using an automated parameter estimation software PEST (Doherty et al., 2010). Using the initial parameter set obtained by a preliminary trial-anderror calibration, the model was first calibrated with PEST by adjusting all model parameters (Table 2a and 2b) to minimize the objective functions, defined as a sum of weighted squared residuals between observed and simulated variables. They included evapotranspiration, snow water equivalent, cumulative snowmelt runoff for each year, and volumetric liquid water content at four soil depths for the calibration of the upland model; and pond water level and annual cumulative runoff for the calibration of depression model. The balance weighting strategy was used in assigning the observation weights to achieve similar contribution from different observation groups. Model sensitivity analysis (see below) was conducted after the first PEST calibration, and a set of parameters that had the most significant influence on model output variables were selected for the second round of PEST calibration to determine the optimal set of parameters.
The sensitivity of model output variables to the variance of model input parameters was evaluated using the global sensitivity analysis method of Morris (1991) implemented in PEST++ suite (Welter et al. et al., 2012). The method uses one-at-a-time approach, in which input parameters are varied in fixed increments to evaluate the elementary effects (EE i ) of each parameter on the variance of model output. For example, the elementary effect of input parameter X i on model output Y is defined as where k is the number of input parameters. The mean (μ) and standard deviation (σ) of EE i and the mean (μ*) of |EE i | are calculated for randomly sampled values of (X 1 , . . . , X k ), and used as the measures of model sensitivity to input parameters (Saltelli et al., 2008, p. 110). In general, a higher value of μ* indicates a higher model sensitivity to the parameter examined (Saltelli et al., 2008).

Dynamical downscaling of climate change scenario
The outputs of general circulation models are too coarse for hydrological applications. Weather Research and Forecasting (WRF) is a numerical weather prediction system used for atmospheric forecasting and dynamical downscaling of lowresolution climate products (Skamarock et al., 2008). Highresolution atmospheric models such as WRF are needed to represent fine-scale processes such as convective summer storms common in the Canadian Prairies . Reliable downscaling of precipitation with complex characteristics and higher spatiotemporal variability can reduce the uncertainty in groundwater recharge estimation.
Despite recent improvements in the representation of atmospheric processes in climate models, there are still large biases (i.e., mismatch) between atmospheric model outputs and observations, which need to be corrected (Teutschbein & Seibert, 2012). Most bias-correction approaches are applied to individual atmospheric variables without taking into account the dependence between different variables. This impairs the spatial and temporal consistency of the atmospheric fields by violating the mass and energy conservation principles (Ehret et al., 2012). This study uses a bivariate quantile mapping method (Cannon, 2018), in which a two-dimensional probability density function is used to correct the model biases for air temperature and precipitation in the control period first and then the future projection. This bias correction method keeps the consistency and dependency between air temperature and precipitation that were previously neglected by conventional methods using univariate bias correction approaches (Maraun, 2013).
The bias-corrected WRF model outputs with 4-km horizontal grids and 37 vertical levels were used under the present climate (HY2007-2015) and under a future warming projection (HY2092-2100) to study the future changes in an uplanddepression hydrological system. Note that a hydrological year (HY) in this study starts on 1 November and ends on 31 October, coinciding with the start of soil freezing in most years.
The 4-km resolution WRF model was set up for a 7.2 × 10 6 -km 2 area over western Canada and was forced with a 0.703˚× 0. simulate the present climate . For the future climate, the WRF model was forced with a future projection of global climate generated from the same ERA-Interim data perturbed by climatological change factors for surface variables (e.g., surface temperature, soil temperature, etc.) and three-dimensional field variables (e.g., temperature, specific humidity, wind, etc.). The change factors were calculated monthly by comparing the base period with the projection of global climate under the representative concentration pathway of 8.5 W m −2 (RCP8.5) for 2085-2100 . This future projection used an ensemble of more than 60 general circulation models (GCMs) in the Coupled Model Intercomparison Project Phase 5 (CMIP5) for HY2092-2100. This approach, referred to as the pseudo global warming (PGW) method, is commonly used to dynamically downscale climate model projections (Kawase et al., 2009). In this paper, PGW runs using HY2007-2015 as the base period are equivalent to what are expected to occur during HY2092-2100.

Introduction of uncertainty in climate change scenario
The PGW approach allows for the representation of important small-scale processes such as convection. However, it requires extremely large computing resources and hence is not suitable for generating many climate scenarios for uncertainty analysis. To further capture the uncertainty in future groundwater recharge, air temperature and precipitation records at the Spyhill site during HY2007-2015 were perturbed using monthly change factors obtained from the GCMs participated in CMIP5 (Taylor et al., 2012). Out of the 39 GCMs that had monthly data readily available, 15 GCMs were selected to cover different models developed independently in different countries and to minimize the similar results from the same GCMs with different configurations (Table 3). Following the standard statistical downscaling methods, monthly projected changes in air temperature (i.e., temperature change factors) were calculated by subtracting HY2007-2015 values from future values under the RCP8.5 scenario over HY2092-2100. The monthly projected changes in precipitation (i.e., precipitation change factor) were calculated as the ratio of future to present precipitation values. Next, temperature change factors were added to observed temperature at the Spyhill site, and observed precipitation was multiplied by precipitation change factors. The perturbed climate change scenarios were then used to run VSMB-DUS and estimate a range of changes in groundwater recharge under future climate scenarios.

Evaluation of the Penman-Monteith PET implementation
Parameters in the drying curve function (Equation 5) had been determined by Hayashi et al. (2010) using the Priestly and Taylor (1972) equation to calculate PET. After the change of PET algorithm to the Penman-Monteith equation (Equation 3), the drying curve parameters were determined by model calibration. Preliminary sensitivity analysis indicated that only C h and C m needed calibration when C n and C r were fixed at 1.0. Therefore, only C h and C m were included in model calibration (Table 2a). Figure 3 shows observed and simulated vapor flux (E) using the calibrated model parameters for HY2007-2015.
Note that E is calculated using Equations 3-6 only when the top soil layer is unfrozen and free of snow cover. Otherwise, snow surface evaporation or sublimation is simulated by the Utah Energy Balance model. Compared with observed data, simulated E captured the seasonality of vapor flux reasonably well (Figure 3a) and had a mean bias error (MBE) of   (Figure 3b, Table 4). These are relatively small compared with the mean growing-season evapotranspiration rate of 2.3 mm d −1 .

Snow, soil, and runoff processes on the upland
The performance of the VSMB upland model was evaluated using the field data from the Spyhill upland west of GP depression (Figure 1) for HY2007-2015. Daily mean air tem-perature in the region generally remains well below 0˚C from early November to mid-March. However, occasional Chinook (i.e., foehn) events cause midwinter melt of shallow snowpacks (Pavlovskii, Hayashi, & Itenfisu, 2019), leading to a patchy snow cover, which varies with vegetation heights and topography. As an example, the Chinook event on 26-27 Jan. 2011 was particularly pronounced with a 2-d average air temperature of 6.7˚C (Figure 4a), which caused complete depletion of snowpack before it started to accumulate again on 29 January. Figure 4b shows a comparison of observed and simulated snow water equivalent (SWE), and Table 4 lists Vadose Zone Journal F I G U R E 4 (a) Observed daily precipitation (precip.) and daily mean air temperature (temp.). Vertical lines indicate 1 January of each calendar year. (b) Observed (obs.) and simulated (sim.) snow water equivalent (SWE). (c) Observed and simulated soil temperatures at 10-cm depth, with simulated temperature represented by the average of 0-to-10-cm and 10-to-20-cm model layers. (d) Soil temperatures at 30-cm depth. (e) Soil temperatures at 60-cm depth. (f) Observed and simulated soil liquid water contents (LWC) at 30-cm depth. (g) Soil liquid water contents at 60-cm depth statistical indicators of goodness of fit: RMSE, MBE, and Pearson correlation coefficient (CC). During the simulation period (HY2007-2015), SWE was measured 77 times at the study site. Compared with observations, VSMB simulated SWE with RMSE and MBE of 7.96 and −2.60 mm, respectively (Table 4). In the context of this study, it is desirable to estimate the peak SWE as it controls the magnitude of snowmelt inputs to the ground. The observed peak SWE ranged from 6 to 34 mm (Figure 4b) with an average of 23.3 mm. Compared with this value, the magnitudes of RMSE and MBE were relatively small. Therefore, the match between simulated and observed SWE is deemed reasonable, considering the high spatiotemporal variability of snowpack in the prairies and the fact that the VSMB model does not include algorithms for wind-driven snow redistribution over the complex landscape (Fang & Pomeroy, 2008).
Simulated and observed upland soil temperatures (Figure 4c, d, e) had similar responses to snowmelt infiltration and soil freezing and thawing. Simulated temperatures of the top soil layers (0-0.2 m) coincided with observations as they reached the freezing point in mid-to late November in all simulated winters (Figure 4c). Simulated and observed temperatures of the top layer reached the thawing point in late March or in early April in all simulated years. The simulated soil temperature had RMSEs ranging from 1.17 to 1.67˚C and MBE ranging from 0.54 to 0.74˚C for the first 1.5 m of the soil where the observation data were available (Table 4) These values are reasonably small in comparison to ranges of observed temperature variability (Figure 4c-e). Figures 4f and 4g show simulated and observed soil liquid water contents at 0.3-and 0.6-m depths. The VSMB overestimated liquid water contents in top soil layers (0-0.4 m) and underestimated them for the deeper layers (0.4-1.0 m). However, the model captured the temporal variability over the simulated period. The simulated liquid water content in the T A B L E 4 Versatile Soil Moisture Budget Depression upland model performance statistics for evapotranspiration (ET), snow water equivalent (SWE), annual cumulative runoff, daily mean soil temperature, and volumetric liquid water content (VLWC) over the simulation period of hydrological years HY2007-2015 0-to-1.2-m zones had RMSE ranging from 0.03 to 0.04 and MBE ranging from −0.01 to 0.00 (Table 4). The magnitude of error was greater in shallower layers, but the model captured soil water dynamics reasonably at all depths, with RMSE and MBE much smaller than observed ranges of variability (Figure 4f-g). Snowmelt runoff over frozen soil occurs in March-April of years with a sufficiently large amount of snow accumulation. Hayashi and Farrow (2014) found that snowmelt runoff was generally small when winter (1 November to soil thaw) precipitation was <150 mm at the Spyhill grassland and nearby sites during HY2003-2013. Similarly, in this study, snowmelt runoff was observed when winter precipitation was greater than 150 mm, which happened in HY2007 (176 mm Snowmelt runoff in the Canadian Prairies is strongly influenced by soil macropores (Pittman et al., 2020;van der Kamp et al., 2003), which is represented by the model parameter f lxm . As expected, simulated snowmelt runoff was most sensitive to f lxm , where the Morris's sensitivity measure μ* indicated that roughly 50% of runoff variability was attributed to f lxm (see Supplemental Figure S2c). The best match between observed and simulated snowmelt runoff was achieved by f lxm = 10.7 mm d −1 with a MBE of −0.74 mm (Table 4), which was relatively small in comparison with the mean F I G U R E 5 Observed (obs.) and simulated (sim.) runoff caused by snowmelt and summer rainfall. Simulated runoff includes contributions from both upland and depression. Observation data were not available (*) for rainfall runoff in 2011, and for both snowmelt and rainfall runoff in 2012 snowmelt runoff of 5.4 mm. However, simulated and observed runoff had substantial differences in individual years (Figure 5), resulting in RMSE of 3.48 mm ( Table 4).
Regardless of f lxm values, no groundwater recharge occurred on the upland, consistent with the results of previous field studies in the Canadian Prairies (Berthold et al., 2004;.

Depression water balance and groundwater recharge
The runoff generated by the upland model was transferred to the depression model ( Figure 2c) to simulate the water balance (i.e., pond water level) and depression-focused recharge. Simulated pond water levels in the four depressions were compared with observations during HY2007-2015. Two examples are shown in Figure 6 to demonstrate the effects of snowmelt runoff (HY2007) and summer runoff (HY2013), whereas simulation results for other years are presented in Supplemental Figures S4-S8. Note that no observable runoff occurred in HY2010 and HY2015. The model was first calibrated for GP depression using all parameters listed in Table 2b, and subsequently calibrated for other depressions using the GP parameter values as the initial estimates.
In HY2007 ponds formed in all depressions on 7 March (Figure 6b-e). The pond persisted in GP until 10 May with two brief interruptions, and the ponds still existed in N12, G16, and G17 when the last measurement was made on 30 April. The VSMB model captured the timing of the snowmelt runoff peaks reasonably well but undersimulated water levels during the snowmelt season (Figure 6b  ing, but pond water levels were not measured during summer 2007. In HY2013, snowmelt ponding occurred for a brief period in GP (25-27 March) and G17 (24-30 March) (Figure 6g, j), whereas no observation data were available for N12 and G16. The depression model overestimated pond depths at GP and G17 in March-April (Figure 6g, j) as a result of the oversimulated upland runoff ( Figure 5). An exceptionally large amount of precipitation fell over a large area of southern Alberta during 23 May-21 June. For example, 262 mm was recorded at the Spyhill site, which is the largest 30-d precipitation recorded during 2003-2020. This caused summer runoff, which flooded N12, G16, and G17 starting on 23 May (Figure 6h-j). For GP depression, observations were only available from 23 July (Figure 6g). In the Canadian Prairies upland soils generally have enough capacity to absorb summer infiltration (e.g., Figure 4f). Therefore, much of the summer runoff in 2013 and other years ( Figure 5) was likely generated in the depression and its fringe areas, where soils are generally moist (Gerla, 1992;Hayashi et al., 1998). This was reflected in the model result, where the majority of runoff (18.7 mm) in May-August 2013 was generated within the depression ( Figure 5). Simulated summer runoff in other years was also mostly generated within the depression ( Figure 5). Note that runoff values shown in Figure 5 are calculated by dividing runoff volumes by the catchment area (A c in Table 1). The MBE (0.57 mm) and RSME (0.94 mm) of simulated annual runoff were reasonably small (Table 4) in comparison with average observed annual runoff of 13.6 mm ( Figure 5).
Area-averaged groundwater recharge was calculated for the catchments of the four depressions (Figure 7). Recharge amounts varied between the depressions depending on the upland-to-depression area ratio, but they had similar patterns of interannual variability. As the amount of depressionfocused recharge in the region increases with snowmelt runoff and summer precipitation (Hayashi & Farrow, 2014), little or no recharge was simulated in years of low runoff and relatively low May-August precipitation, such as 2010 and 2015. The average annual recharge over the 9-yr simulation period was 10.2 mm yr −1 (Figure 7), which was comparable with 10 mm yr −1 based on the detailed water balance during HY2007-2013 at the study site (Hayashi & Farrow, 2014). Considering the simplicity of the VSMB-DUS and the model parameter uncertainty, the match between estimated and observed recharge rates is deemed reasonable.

Bias correction of WRF for the present climate
The study area is located near the foothills of the Canadian Rockies with rapidly changing weather conditions and frequent deep convection events. Even the convection-permitting WRF model cannot resolve the convection issue completely and shows biases when its outputs are compared against observations. The WRF model oversimulated precipitation (i.e., wet bias), especially in June and July, and consistently undersimulated air temperature (i.e., cold bias) compared with observed data at the Spyhill site during HY2007-2015 (Figure 8). Modeled recharge is particularly sensitive to wet bias, and the use of biased precipitation data can lead to substantial oversimulation of recharge. Therefore, the biases in air temperature and precipitation were simultaneously corrected using the bivariate quantile mapping method (Cannon, 2018), and the same bias corrections were applied to the WRF projection for the PGW scenario representing the future climate at the Spyhill site during HY2092-2100. The PGW conditions were warmer and wetter than the historical observation in HY200-2015 (Figure 9).

Soil and hydrological changes under the PGW scenario
Mean annual precipitation increased from 491 mm under the present climate to 626 mm under the PGW scenario (Figure 9b). The largest change was in April-August, with 5mo total precipitation increasing from 349 to 438 mm. Mean annual air temperature increased from 4.0 to 8.9˚C. The VSMB model differentiates precipitation to snow and rain based on air temperature (T a ): all snow for T a < −1˚C and all rain for T a > 3˚C (Tarboton & Luce, 1996). The number of days with mean air temperature below 0˚C decreased from 142 d in the present climate to 69 d in the PGW scenario (Figure 9a). As a result, the number of snowfall events decreased, and precipitation that fell as snow in early winter (November) and spring (March) under the present climate fell as rain under the PGW scenario. The increased precipitation and fluctuating air temperatures near 0˚C from mid-November to early F I G U R E 1 0 Simulated daily mean soil temperature and total soil water content at the Spyhill site from November to October, averaged over the 9-yr period.  (Figure 9) had a large influence on simulated hydrological processes. Figure 10 shows simulated daily soil temperatures and total water contents during HY2007-2015 (present) and HY2092-2100 (PGW) using the average of nine hydrological years for each day to display a general pattern. Note that Figure 10 averages out the interannual variability of soil temperature, resulting in an underrepresentation of the actual days with temperature below 0˚C (i.e., frozen soil) for individual years. The soil was generally frozen down to a depth of 1.2-2.0 m under the present climate (Figure 10a). The number of days when the soil was frozen at some depth (i.e., frozen-soil days) varied between 121 d in HY2010 and 205 d in HY2014 with an average of 173 d under the present climate. Under the PGW scenario, the soil froze down to 0.2-0.4 m and had freeze-thaw cycles in midwinter, which did not occur under the present climate (Figure 10a). The number of frozen-soil days varied between 29 and 146 d with an average of 74 d. An increase in mean annual air temperature of 4.9˚C under the PGW was reflected in temperatures of the entire soil column, with particularly strong warming in March-May (Figure 10c).
Snowmelt infiltration in midwinter under the present climate was limited to the top 0.2 m (Figure 10d) due to the presence of frozen soil. In contrast, midwinter infiltration of snowmelt and rain reached down to 0.6 m in the PGW scenario ( Figure 10e). When the present climate and the PGW was compared (Figure 10f), total water contents in the shallow zone (<0.4 m) decreased substantially in April due to the smaller amount of snowmelt infiltration. The increase in April-August precipitation (Figure 9b) did not increase soil water content (Figure 10f) due to an increase in plant water uptake in response to the increased PET. The PET during the growing-season (April-September) was 488 mm under the present climate and 583 mm under the PGW. Figure 11 shows the simulated snowpack for the present climate (HY2007-2015) and the PGW scenario (HY2092-2100), in which the warmer climate substantially contributed to decreasing snowpack. Relatively rapid melt of snowpack in December and January under the PGW scenario resulted in a sizable midwinter snowmelt runoff in some years-for example, in HY2099 ( Figure 11)-that was not observed under the present climate. Except for these midwinter runoff events,

Uncertainty in recharge under climate change scenarios
The 15 perturbations of future climate scenarios had a range of air temperature ( Figure 9a) and precipitation (Figure 9b). Of the 15 perturbations, seven of them had lower or marginally higher (<50 mm) mean annual precipitation than the present climate (Table 3). These low precipitation scenarios resulted in no recharge, presumably because of increased evapotranspiration resulting from higher PET. The remaining eight perturbation cases had a wide range of recharge (Figure 12), and the mean of all perturbations for each year was considerably smaller than under the present climate (Figure 7). Based on these observations, it is likely that the reduction in recharge simulated by the PGW scenario represents an overall trend of an ensemble of climate change scenarios represented by the models included in CMIP5.

Climate change effects on depression-focused recharge
There is a broad consensus about increases in air temperature in the Canadian Prairies . This study uses the HY2092-2100 climate scenario based on RCP8.5 representing the upper end of temperature prediction to demonstrate the maximum potential effects of warming. Not surprisingly, the model forced by this scenario predicts a major reduction in snow accumulation ( Figure 11) and frozen-soil period (Figure 10), which is consistent with the findings of previous studies (Mellander et al., 2007). However, the effects of the warmer and wetter climate on runoff and recharge are difficult to predict due to the complex feedback of multiple processes-for example, an increase in precipitation may be compensated by an increase in evapotranspiration (Rasouli et al., 2019).
Under the present climate, groundwater recharge is focused in depressions, whereas no recharge occurs in uplands because infiltration is completely used up by evapotranspiration (Figure 10a). The PGW simulations show that recharge remains focused in depressions despite a substantial increase in precipitation under the warming scenario (Figure 9b). As a result, the response of recharge fluxes to the warmer and wetter climate is more complex than the cases of diffuse or blanket recharge reported in previous studies that predict increased recharge with increased precipitation (Crosbie et al., 2013;Okkonen et al., 2010;Zhang et al., 2020).
Given the continued dominance of depression-focused recharge, snowmelt runoff will continue to be an important factor under the future climate. The PGW simulations indicate an increased frequency of midwinter melt events (Figure 11), which partially shifts snowmelt runoff from March-April to December-February. The PGW simulations also indicate the occurrence of March-April runoff due to rain on unfrozen soil and a substantial increase in May-August runoff. This is partly due to an overall increase in precipitation amount (Figure 9b), but also due to the increased intensity and duration of summer storm events under the warmer climate (Shook & Pomeroy, 2012). For example, the large simulated runoff in August 2099 is caused by a large storm event with 170 mm of rain falling in 5 d, which is extremely rare under the present climate. These results imply that the timing of depressionfocused recharge will shift later in the year under the warmer climate.

Model limitation and uncertainty
This study is a first attempt to evaluate the potential effects of climate change on depression-focused groundwater recharge in cold semiarid environments, where snow and frozen soil play a major role in hydrology. It is founded on two decades of detailed field studies of depression-focused recharge in the prairies of Saskatchewan (Hayashi et al., 1998Berthold et al., 2004) and Alberta , and incremental model improvements incorporating field-based understanding of snow and frozen soil processes (Mohammed et al., 2013;Noorduijn et al., 2018). However, the simple one-dimensional algorithms of VSMB-DUS are unable to capture all aspects of complex feedback and their responses to climate change. For example, VSMB-DUS uses simple soil water balance to simulate soil moisture distribution in the soil column, and the upland-depression interaction is represented by one-way transfer of runoff ( Figure 2). This is far less sophisticated than fully three-dimensional integrated hydrological models (Kollet et al., 2017), as it lacks the ability to fully couple soil water processes with groundwater processes considering both horizontal and vertical fluxes. The lack of coupling could result in overestimation of net recharge, as it cannot simulate the removal of groundwater by phreatophytes and other deep-rooted plants surrounding depressions (van der Kamp & Hayashi, 2009). However, among many integrated hydrological models, only a few can simulate soil freeze-thaw and its effect on infiltration in a fully coupled, physically consistent manner (Jan et al., 2020). In addition, 3D solutions of the Richards equation with freeze-thaw are computationally demanding due to the additional nonlinear feedback between water and energy transfer. Therefore, VSMB-DUS offers a reasonable compromise for simulating thousands of upland-depression systems for watershed-scale recharge estimation (Klassen et al., 2018). Noorduijn et al. (2018) compared simulations of depressionfocused recharge by VSMB-DUS and Hydrus2D (Šimůnek et al., 2006), which solves the 2D Richards equation, and showed that simulated recharge values were nearly identical.
To demonstrate the usefulness of VSMB-DUS, it was used to estimate the spatial variability of groundwater recharge under the present climate in the agricultural region of Alberta, Canada (∼190,000 km 2 ), using the statistical upscaling methodology developed by Pavlovskii et al. (2020) to represent a wide variety of depression-catchment area ratios (A d /A c ) and retention capacity. The model sensitivity analysis indicated that simulated groundwater recharge is most sensitive to A d , (Supplemental Figure S3). The statistical upscaling approach incorporates the uncertainty in A d by running VSMB-DUS for many combinations of A d and A c and calculating an average recharge for various surficial sediment types (Pavlovskii et al., 2020).
Over 7,000 sets of VSMB-DUS simulations were conducted following the methods described in the supplemental material for the unmanaged grassland setting using the soil and crop parameters listed in Tables 2a and 2b. Simulation results indicate a regional gradient of recharge from roughly <5 mm yr −1 in southern parts of the region having drier and warmer climates to >20 mm yr −1 in northwestern parts of the region (Figure 13). This type of information is useful for regional-scale groundwater management.
In addition to the model simplicity and parameter uncertainty, a major limitation of this study is its focus on a single land use, namely unmanaged grassland. Previous studies in the Northern Prairies have shown that land uses have a strong influence on the partitioning of snowmelt and rain into runoff and infiltration (Renton et al., 2015;van der Kamp et al., 2003). Therefore, the findings of this study may not be directly transferrable to other land-use types such as cereal crops or managed pasture. Nevertheless, the basic understanding obtained from the study is expected to be applicable to other land uses-for example, earlier snowmelt, shorter frozen-soil period, and more intense summer storms will F I G U R E 1 3 Spatial distribution of mean annual recharge for the unmanaged grassland setting simulate by the Versatile Soil Moisture Budget Depression-Upland System (VSMB-DUS). Recharge in each polygon is simulated using the data from a weather station (WS) representing the polygon. Variability within polygons reflect the difference in surficial sediment types (see Supplemental Figure S1) have similar hydrological effects among different land-use types.
Another omission in this study is changes in plant species and phenology under the warming climate. The VSMB model simulates plant phenology using a simple approach based on air temperature, soil temperature, and latitude (i.e., potential duration of sunshine hours). The model simulations indicate that the average start date of the growing season shifts from 26 April under the present climate to 18 March under the PGW scenario, whereas the end of the growing season remains unchanged on 18 October. A longer growing season is expected under the warmer climate. The VSMB growthstage algorithm (Baier & Robertson, 1966) has been successfully used in the Canadian Prairies for many decades, but its applicability has not been tested for a change in mean annual temperature as large as 4.9˚C and will require further examination.

CONCLUSIONS
Previous field studies have shown that groundwater recharge in the Northern Prairies is focused under numerous topographic depressions that receive lateral inputs of snowmelt runoff from uplands. This study used a process-based soil water balance model, VSMB-DUS, to simulate hydrological processes controlling depression-focused groundwater recharge. The VSMB-DUS model captured the seasonality of hydrological processes reasonably well in comparison with field observation. Simulated snowmelt runoff was sensitive to the model parameter representing macropore flow in frozen soil, implying the importance of characterizing the hydrological functions of soil macropores and its dependence on land use. The effect of climate change on groundwater recharge was evaluated using the PGW approach representing a climate of 2092-2100 under the RCP8.5 scenario. The dynamically downscaled and bias-corrected climate model outputs indicated that the PGW scenario had substantially higher mean annual air temperature (+4.9˚C) and precipitation (+135 mm or +27%) compared with the present climate. Such a large increase in precipitation would increase diffuse (or blanket) recharge in humid regions. However, the model simulations in this study demonstrated that the complex interplay between increasing air temperature and precipitation would reduce the amount of depression-focused recharge. This is primarily because the warmer winter generates less snowmelt runoff, and the warmer summer increases potential evaporation, which counters the effects of higher precipitation. The reduction in snowmelt runoff is associated with smaller accumulation of snowpack and shorter periods of frozen soil under the warmer climate.
This study demonstrates the importance of the nonlinear nature of groundwater recharge in the Northern Prairies, where the occurrence of recharge requires the lateral transfer of runoff from uplands, which allows the infiltration in depressions to exceed evapotranspiration demands. The climate change effects on snow and frozen soil processes are complex, and accurate simulation of snowmelt runoff over frozen soil is particularly difficult. Further research is required to improve the representation of these processes in recharge models and incorporate the effects of changes in land use, including crop types and phenology under the future climate.

A C K N O W L E D G M E N T S
We thank numerous field technicians and students who assisted with field data collection at the Spyhill site, Alberta Infrastructure for the land access, Larry Bentley and Edwin Cey for helpful discussion, Igor Pavlovskii and Polina Abdrakhimova for GIS assistance, Aaron Mohammed for soil moisture data, Yanping Li for WRF output data, Alex Cannon for help with bias correction, and Daniel Itenfisu for meteorological data used for WRF bias correction and regional-scale recharge simulation. Funding support for the study was provided by Alberta Innovates, Alberta Environment and Parks, Alberta Energy Regulator/Alberta Geological Survey, Global Water Futures (Prairie Water), and the Natural Sciences and Engineering Research Council. Constructive comments by two anonymous reviewers and the associate editor improved the quality of the article.

C O N F L I C T O F I N T E R E S T
There are no conflicts of interest.