Soil water recharge for grassed and forested land covers on the Oak Ridges Moraine, southern Ontario, Canada

Land cover affects recharge (R) through its influence on water delivery to and evapotranspiration (ET) from the soil (Healy 2010) and its potential to modify soil properties controlling water storage and transfer (Greenwood & Buttle 2014a). For example, grasslands and forests vary in their interception (I), partitioning and use of precipitation and soil water. Reforestation of grasslands increases I (Vose et al. 2011) and decreases throughfall (TF) and snow accumulation (Gelfan et al. 2004). Evaporation of intercepted rainfall (Le Maitre et al. 1999) and sublimation of intercepted snow (Gelfan et al. 2004) in forests is insignificant in grasslands. This, coupled with greater transpiration (Farley et al. 2005), leads to enhanced ET (Zhang et al. 2001; Ladekarl et al. 2005) and reduced R in forests (Ladekarl et al. 2005; Nosetto et al. 2005; Wattenbach et al. 2007).

Forest composition and age also influence canopy water partitioning, ET and R. Annual I for conifers exceeds that for deciduous species (Houle et al. 1999; Barbier et al. 2009), reducing water delivery to the soil. This difference is most pronounced in winter when most conifers retain their needles (Pomeroy et al. 1998; Houle et al. 1999; Komatsu et al. 2008), such that snow accumulation and water availability for R is greater in deciduous forests (Komatsu et al. 2008; LaMalfa & Ryle 2008). Branch architecture influences the spatial distribution of TF and stemflow (SF), with implications for R. Branches sloping toward the tree bole (typical of many deciduous species) promote TF focusing close to boles (Carlyle-Moses & Price 2006) and large SF fluxes (Voigt 1960). Conversely, trees with branches sloping away from the bole (typical of many conifers) have more uniform sub-canopy TF (Loustau et al. 1992) or TF focusing toward the crown edge (Johnson 1990) along with less SF (Barbier et al. 2009). Focusing of TF and SF near deciduous tree boles may enhance R via root-induced macropore flow (Johnson & Lehmann 2006). Stand age for a given forest species can affect I and ET with implications for water availability for R. Conifer I may increase with stand age (Bryant et al. 2005; Buttle & Farnsworth 2012), and Rosenqvist et al. (2010) noted increased I and declining R for ageing Norway spruce stands with no change in either I or R with age for an oak stand.

While hydrological implications of land use changes such as grassland conversion to forest plantations are broadly understood, specific consequences for R are generally unknown despite their important ecological and water resource implications (Scanlon et al. 2005; Zhang & Schilling 2006). Such consequences are particularly significant on the Oak Ridges Moraine (ORM) in southern Ontario. This key hydrogeologic feature supplies potable water to many of the region's residents as well as key aquatic ecosystems services (Howard et al. 1995; Buttle et al. 2015), and accurate R estimates are critical to land use planning in the context of the ORM's groundwater resources (Gerber & Howard 2002). Nevertheless, the influence of land cover on R on the ORM has received relatively little attention (but see Meriano & Eyles 2003). This paper examines point-scale R for major land covers (grasslands, mature mixed hardwood forests and red pine plantations of differing ages) in a managed forest (the Ganaraska Forest, (GF)) on the crest of the ORM by testing the following hypotheses:

1. Annual R will increase in the order old red pine plantation (OP) → young red pine plantation (YP) → mixed hardwoods (MH) → grassland (OPEN), consistent with greater annual I for conifers relative to hardwoods and similar soil ET for all land covers.

2. Dormant season R will increase as in (1) given greater I in pine plantations relative to deciduous stands, whereas growing season R will increase from MH → OP → YP → OPEN, assuming greater growing season I for deciduous species relative to red pines and for OP relative to YP in the GF (Buttle & Farnsworth 2012) and similar ET from the soil for all land covers.

Recharge below 1 m depth was estimated at 10 sites from June 23, 2009 to July 12, 2010: two grasslands (OPEN1 – abandoned pasture, OPEN2 – hay field), three mixed hardwood (MH) stands, three young RP (YP) and two old RP (OP) plantations (Figure 1; Table 1). All were on outcropping glaciofluvial sands and gravels on either Pontypool sand or Pontypool gravely sand soils, and were on level ground to minimize potential for lateral water flow. Areas with minimal undergrowth were selected to limit the effects of undergrowth interception and increased permeability associated with undergrowth root networks (Greenwood & Buttle 2014a).

Species, basal area, height, density, leaf area index and canopy density was measured for forest stands (Buttle & Farnsworth 2012; Table 1). Vertical profiles of bulk density, porosity and soil texture were obtained from horizontal cores (6 cm long, 4.8 cm i.d.) at 0.1-m intervals to 1 m from a 1.5 m deep soil pit in the centre of the site, with the exception of OPEN2, YP1 and OP3 where profiles were from a vertical core to 1 m depth at the centre of the stand and sectioned every 0.1 m. Analytical methods are described in Buttle et al. (2014). A single ring infiltrometer and Guelph Permeameter (Reynolds & Elrick 1985) were used to determine saturated hydraulic conductivity (K_H) at the surface and at depth (every 0.1 m to a depth of 0.8 m), respectively.

Above-canopy precipitation (P_g) was measured at OPEN1 and OPEN2 and at two other clearings with a standard rain gauge (Figure 1) while TF was the average catch of 10 randomly sited rain gauges in each forest site (Buttle & Farnsworth 2012). Ten MH trees at MH1 and four RP trees at both YP1 and OP3 (Table 2) were used to estimate SF as described by Buttle et al. (2014). All SF infiltrated the soil within ∼0.2 m of the bole (Buttle et al. 2014); thus, SF volume was divided by the area of a 0.2 m wide ring around each bole to obtain SF depth (mm). Measurements of TF and SF were not made during winter due to the inability to regularly access some sites, and SF was assumed to be negligible.

Above-canopy ET in 2009 was from the MOD16 global ET product (Mu et al. 2011), which simulates spatially distributed 8-day ET at 1 km² resolution. The MOD16 grid was overlain on LANDSAT 7 imagery to identify pixels that best represented each land cover, and pixels composed predominantly of a single land cover were selected to minimize convolution from mixed pixels. One pixel was used for MH while two pixels were averaged for OPEN ET. Nine pixels covering a 3 × 3 km area and dominated by RP were examined. Young or old RP or white pine stands could not be differentiated, so pixels likely contained a mix of the three types. Pixels with total ET in 2009 > 1 standard deviation above the mean were excluded as outliers, and 8-day ET values for remaining pixels were averaged to estimate ET for YP and OP stands. Estimates were divided by eight to obtain daily ET. Canopy I was subtracted from MOD16 ET to estimate direct soil evaporation plus transpiration. MOD16 data were unavailable for 2010, and regression relationships between 2009 MOD16 ET for the land covers and 2009 potential ET from the Hamon model (Dingman 2002) were used to predict 2010 ET. Temperature data for ET estimation were obtained from a BML-TS-7 Thermilinear Air Temperature Sensor at a clearing used to measure P_g (Figure 1).

Cumulative R was estimated for a 384-day period subdivided into growing season, fall and winter–spring. Growing season was June 23–October 1, 2009 and May 5–July 12, 2010. The fall period was October 2–December 4, 2009 (decreasing MH leaf area index and ET for all land cover), while winter–spring was December 5, 2009–May 4, 2010 (MH leaf area index approaching its maximum). Recharge was estimated for snow-free and snow-covered conditions.

Basal area was estimated in two or three 314 m² circular plots in each stand, and the 0.2 and 0.2–0.5 m focusing zones areas around each tree were calculated. Weighting factors a and b ranged from 0.012 to 0.022 and 0.036 to 0.073 in MH stands, respectively, while a ranged from 0.012 to 0.063 in RP stands (b = 0 given no focused TF).

Total P_g for July 1, 2009–June 30, 2010 (1,199–1,335 mm, Table 3) was well above the 840 mm annual average at Peterborough Airport (Environment Canada 2011) ∼30 km to the northeast. Total snowfall for the 2009–2010 winter was 82 cm at Tapley, approximately half the normal depth. Mean TF as a fraction of P_g was similar between forest sites (0.8–0.87); however, MH TF:P_g ratios increased when TF focusing near the bole was considered. MH TF:P_g ratios increased between the growing season and winter–spring; conversely, TF:P_g for YP and OP generally decreased between these periods, reflecting greater snow interception efficiency. MH stands generated more SF than either YP or OP, whether SF was expressed relative to stand area or the 0.2 m wide focusing zone adjacent to the bole. Both TF focused within 0.5 m of the bole and SF focused within 0.2 m of the bole exceeded P_g for MH stands.

Maximum SWE for all but one site (OPEN2) occurred at the beginning of March 2010. Snowmelt ended by March 18, 2010 with the exception of MH2 where it ended a week later. OPEN and MH sites generally had more SWE than RP sites (peak SWE of 95–125 mm at the former vs. 48–90 mm at the latter), with least SWE at OP sites. Differences in peak SWE for RP relative to OPEN sites suggest I_s between 25 and 40% of P_g.

Period I ranged from 159 mm (YP3) to 260 mm (OP3), representing 13% and 19% of P_g, respectively (Table 3). Mean MH, YP and OP I for the study period was 220 ± 31 mm, 177 ± 30 mm and 227 ± 47 mm, respectively.

There was a general increase in R:P_g from growing season to winter–spring (Table 5), and 55–100% of P_g recharged during the latter. Growing season R and R:P_g at OPEN and MH sites generally exceeded RP values. MH R:P_g ratios were greater relative to OPEN sites, reflecting coarser soil textures at the former (Figure 5). There were no clear distinctions between fall or winter–spring R and R:P_g at OPEN-MH relative to RP sites. Sites with marked increases in S_t in fall (OPEN1, OPEN2, MH3, YP2 and YP3, Figure 6) converted less P_g to R relative to sites with little storage change. Conversely, sites with the smallest fall R:P_g values had larger winter–spring R;P_g values, consistent with rapid declines in S_t following snowmelt (Figure 6).

Despite greater total R and R:P_g from OPEN and MH relative to RP sites (Table 4, Figure 8), Z-scores indicated few significant inter-site differences in R (Table 6). OPEN1 had greater R than two RP sites, while MH3 R exceeded that at OPEN and RP sites. There were no significant differences in R or R:P_g between sites in a given land cover, or between YP and OP. Inter-site R differences were largely driven by growing season conditions, with greater R at OPEN1 relative to YP1 and OP1, and at MH3 relative to all RP sites. Greater fall R at sites with little S_t change relative to those with marked S_t increases (Table 5) was often balanced by less R at the former during winter–spring.

Hypothesis 1 was partly supported (Table 4) by greater total R from OPEN and MH (594–794 mm) relative to RP stands (480–565 mm), consistent with previous studies. Forests have greater water use (and less R) relative to grasslands (Houghton-Carr et al. 2013), and Dripps & Bradbury (2010) found more R from grasslands compared to coniferous forests but greater R from coniferous relative to hardwood forests due to the latter's greater ET. However, they did not account for focused TF and neglected SF when estimating hardwood forest R. Conversion of pine to oak forest increased R by 4.8% largely due to decreased ET from the latter (Wattenbach et al. 2007), while oak forest R exceeded that from spruce forests approximately 20–25 years following reforestation (Rosenqvist et al. 2010). However, the large associated errors meant that only OPEN1 and MH3 had total R significantly greater than some (YP2 and OP1 in the case of OPEN1) or all (MH3) of the RP stands (Table 6). Greater total R at MH3 partly resulted from larger P_g and TF inputs (Table 3); nevertheless, inter-site differences in P_g and ET exerted a secondary control on R variations between land covers. Thus, total R was ∼800 mm at MH3 and <600 mm at OP3 despite identical P_g (Figure 8(b) and 8(d)), while total period ET was only ∼70 mm greater at RP relative to MH stands (Figure 3).

Annual R for land covers was assumed to be primarily influenced by I, with greatest I for mature conifers and least for non-forested areas (Houle et al. 1999; Barbier et al. 2009). However, R did not follow the hypothesized increase in water inputs along this land cover gradient, with no significant differences in I:P_g ratios between OP, YP and MH stands (Table 3). Thus, while water inputs to the soil were greatest at OPEN sites (1,205–1,239 mm), there was considerable overlap between inputs (TF + SF) at MH (960–1,106 mm), YP (1,027–1,080 mm) and OP (1,005–1,075 mm) stands.

This hypothesis also presumed similar ET from the soil across land covers (Sun et al. 2008). Study period above-canopy ET was least for OPEN (564 mm) and greatest for RP stands (679 mm), while maximum daily total ET ranged from 3.2 mm d⁻¹ (OPEN) to 4.4 mm d⁻¹ (MH). Forest ET often exceeds that for grasslands (Eugster & Cattin 2007; Houghton-Carr et al. 2013), and greater ET for RP relative to MH stands is consistent with Mackay et al. (2002). However, much of the greater RP and MH ET relative to OPEN sites resulted from evaporation of intercepted water. Soil ET when adjusted for I was greatest for OPEN sites (564 mm), intermediate for RP stands (467–550 mm), and least for MH stands (362–427 mm). Total R suggests larger net precipitation inputs more than compensated for greater ET from the soil at OPEN sites.

MH and OPEN R values were similar, partly as a result of less soil ET at MH sites. Greater MH R relative to RP stands also reflects focused SF and TF near MH boles. Many forest recharge studies assume spatially uniform TF (Dripps & Bradbury 2010) and SF (Neary & Gizyn 1994) inputs to the soil, or neglect SF altogether (Dripps & Bradbury 2010). While this is reasonable in RP where SF contributions to R were negligible (Table 4), MH branch architecture focuses TF and SF towards the bole (Carlyle-Moses & Price 2006; Sato et al. 2011; Bialkowski & Buttle 2015). MH TF_F was >2× the average catch of the 10 TF gauges, and MH generated more SF relative to RP stands. MH SF relative to the 0.2 m wide ring around the bole (Table 3) produced SF fluxes ∼ 1.3 × P_g. These were less than the focused SF:P_g ratio of 2.7 for beech trees (Voigt 1960), which are more efficient SF producers than maple trees (Carlyle-Moses & Price 2006). Nevertheless, focused near-bole TF and SF represent an important R source (Johnson & Lehmann 2006). Although areas for R_0.5 and R_0.2 estimations were only ∼ 4–7% and 1–2% of total stand area, respectively, their inclusion increased MH R:P_g by 11–18%. This highlights the importance of including focused inputs in R studies in MH forests.

Inability to monitor SF in all stands during winter–spring led to the assumption of negligible SF for this period. While this may be reasonable for RP (Table 3), Herwitz & Levia (1997) and Levia (2004) showed significant winter SF generation in deciduous forests. Inclusion of focused SF during this period would have increased total (Table 4) and winter–spring (Table 5) MH R; however, the increase is of unknown magnitude and needs further study.

Omission of such SF inputs may partly explain why Hypothesis 2 was not supported. There was no consistent difference in R between land covers in fall and winter–spring, despite considerable I_s in RP sites consistent with other findings (Pomeroy et al. 1998; Houle et al. 1999) and reductions in rainfall I between leaf-on and leaf-off periods for deciduous species (Gerrits et al. 2010). Larger peak SWE in OPEN and MH sites relative to RP stands was also expected to enhance R (Komatsu et al. 2008; LaMalfa & Ryle 2008); however, consistent differences between land covers were not seen. Differential ability of sites to store and release soil water also contributed to the failure to support Hypothesis 2. Despite similar soil types, sites showed considerable variability in soil texture and hydrological properties (Figures 3, 4 and 5). All access tubes were >1 m from the bole, and observed changes in S_t should reflect a site's soil textural and hydrologic properties rather than focused TF and SF inputs (Buttle et al. 2014). Thus, some sites (OPEN1, OPEN2, MH3, YP2, YP3) showed marked increases in fall S_t regardless of land cover, and had relatively less R compared to sites with limited storage response and greater water transmission through the upper 1 m of soil (Table 5). Conversely, the former sites were relatively wetter prior to snowmelt, suggesting greater unsaturated K_H compared to drier sites and more effective transmission of water inputs through the profile. Drainage of soil water supplemented R at these sites, leading to large R:P_g ratios (Table 5). Results highlight the importance of considering soil type in land cover comparisons of R, and the potential for inaccurate R estimates using low resolution soils maps (Faust et al. 2006; Crosbie et al. 2009).

Growing season MH R was ≥ OPEN sites, with no consistent differences between YP and OP stands. Growing season TF and TF:P_g were generally smallest (and I was largest) in MH relative to RP stands (Table 3), suggesting that R should be least at MH stands. Minimum R:P_g might also be expected to occur at all sites during the growing season due to greater ET during this period (Figure 7; Faust et al. 2006; Dripps & Bradbury 2010). Instead, the greatest fraction of total MH R occurred in the growing season relative to fall and winter–spring. This contrasted with OPEN and RP sites (with the exception of OP1) where winter–spring R contributed the greatest portion of total R. This may reflect the wet summer in 2009 and 2010, combined with enhanced MH R through near-bole SF and TF focusing (Table 4). However, the caveat that unmonitored winter–spring SF in MH stands may have increased R during this period must be considered.

Many hardwood-conifer forests in the GF were converted to agriculture in the 19th century, and the pine plantations intended to suppress soil erosion from agricultural fields and reduce downstream flooding are currently managed to promote re-conversion to mixed hardwood stands (Buttle 2011). Annual R during this land cover transition may decline from high values for grasslands (mean annual R = 575 mm) to minima under red pine plantations (mean annual R = 375 mm), and subsequently increase for mixed hardwoods (mean annual R = 616 mm). The shift in R during this land cover change (∼200 mm) represents ∼ 20% of annual precipitation and 50% of the ∼400 mm of average annual recharge for exposed sand and gravel deposits on the ORM (Gerber & Howard 2000, 2002), with important implications for groundwater resources. This may also be accompanied by a change in the seasonality of R, with winter–spring dominating annual R for grasslands and red pine plantations, but with growing season R of increasing importance for mixed hardwoods partly as a result of enhanced R due to focused near-bole inputs of TF and SF.

These R values were for level sites, and superimposed on the suggested temporal trajectory of R during land cover change are topographic controls on the spatial pattern and magnitude of R. Depression-focused R in agricultural areas along the ORM's crest can greatly exceed R on level sites (Greenwood & Buttle 2014b), attributed to frozen winter soils and overland flow of spring snowmelt and rainfall to topographic depressions (Hayashi et al. 2003). Its apparent absence in RP plantations and mixed hardwood stands likely results from smaller SWCs prior to freeze-up combined with warmer soil. Thus, the magnitude and spatial pattern of R on the permeable soils mantling the ORM's crest may depend on the combined effects of topography and land cover, and merit further study.

Accuracy of recharge estimates from 1-d water balances such as those used here depends on the exactitude with which other water balance components have been determined (Healy 2010). Nevertheless, recharge values for major land covers in the GF are in the order of the 400 mm of mean annual recharge estimated for outcropping sand and gravel deposits along the ORM's crest. Differences between maximum and minimum recharge for these land covers represent ∼30% of mean annual precipitation and up to 66% of mean annual recharge on the ORM's outcropping sands and gravels, and should be considered when evaluating the ORM's water resources. These differences were not simply a function of inter-site contrasts in precipitation or the ability of different forest types to intercept and subsequently evaporate precipitation. Large errors associated with the recharge values complicated efforts to discern differences between land covers in the GF, as did inter-site variations in the ability of nominally similar soil types to store and transmit water. Nevertheless, annual recharge for grasslands and mixed hardwood stands exceeded that from young and old red pine plantations. This suggests that historical land cover changes along this portion of the ORM from mixed hardwood-conifers → agricultural fields → red pine plantations → regenerating mixed hardwoods were accompanied by shifts in annual recharge. Focused throughfall and stemflow inputs near the boles of mixed hardwoods appear to promote differences between land covers in terms of the fraction of total recharge supplied in growing vs. dormant periods, suggesting that changes in the seasonality of recharge may also accompany this land cover transition.

Recharge values reported here are among the first direct estimates for major land covers on the ORM, and highlight the need to consider the role of land cover when assessing spatial variations in recharge along the ORM. They should also assist efforts to model groundwater flow and understand inter-basin differences in streamflow for this key hydrogeologic feature in southern Ontario.

This work was supported by the Natural Sciences and Engineering Research Council of Canada. We thank the Ganaraska Region Conservation Authority for continued support; Peter Lafleur, Andrew Farnsworth, Clément Le Saux and Alexandra Ryland for their assistance; Ken Hill and Marie and John Toon for access to their property; Qiaozhen Mu (University of Montana) for the MOD16 data; and two anonymous reviewers for comments on an earlier version of this paper.