Estimation of temporal and spatial variations in groundwater recharge in unconfined sand aquifers using Scots pine inventories

Climate change and land use are rapidly changing the amount and temporal distribution of recharge in northern aquifers. This paper presents a novel method for distributing Monte Carlo simulations of 1-D sandy sediment profile spatially to estimate transient recharge in an unconfined esker aquifer. The modelling approach uses data-based estimates for the most important parameters controlling the total amount (canopy cover) and timing (thickness of the unsaturated zone) of groundwater recharge. Scots pine canopy was parameterized to leaf area index (LAI) using forestry inventory data. Uncertainty in the parameters controlling sediment hydraulic properties and evapotranspiration (ET) was carried over from the Monte Carlo runs to the final recharge estimates. Different mechanisms for lake, soil, and snow evaporation and transpiration were used in the model set-up. Finally, the model output was validated with independent recharge estimates using the water table fluctuation (WTF) method and baseflow estimation. The results indicated that LAI is important in controlling total recharge amount. Soil evaporation (SE) compensated for transpiration for areas with low LAI values, which may be significant in optimal management of forestry and recharge. Different forest management scenarios tested with the model showed differences in annual recharge of up to 100 mm. The uncertainty in recharge estimates arising from the simulation parameters was lower than the interannual variation caused by climate conditions. It proved important to take unsaturated thickness and vegetation cover into account when estimating spatially and temporally distributed recharge in sandy unconfined aquifers.


Introduction
Eskers are permeable, unconfined sand and gravel aquifers (Banerjee, 1975).In addition to water supply, they support groundwater-dependent ecosystems and provide recreational services (Kløve et al., 2011).Esker hydrology is important as eskers and other glaciofluvial aquifer types cover large areas of the north and are among the dominant aquifer types in the boreal zone.Management of these complex aquifers has gained recent attention (Bolduc et al., 2005;Karjalainen et al., 2013;Koundouri et al., 2012;Kurki et al., 2013).The European Groundwater Directive requires such systems to be characterized in order to determine their quality status, so knowledge of how to estimate groundwater recharge in esker aquifers is becoming increasingly important (EC, 2006).Esker aquifers are commonly covered with managed pine forests, where the forest canopy is likely to influence recharge amounts.The soil surface profile of eskers is complex and highly variable, consisting of kettle holes and sand dunes, resulting in variable thickness of the unsaturated zone (Aartolahti, 1973), a factor which also needs to be accounted for in recharge estimation.
Computational methods to estimate groundwater recharge vary from simple water balance models, where water stores and fluxes are represented conceptually and related with adjustable parameters (Jyrkama et al., 2002), to physically based models using the Richards' equation (Assefa and Woodbury, 2013;Okkonen and Kløve, 2011) to solve water fluxes through an unsaturated zone.Computational methods solving the Richards' equation are often limited to smallscale areal simulations (Scanlon et al., 2002a) and shallow unsaturated zones, and they commonly lack the soil freeze, thaw, and snow storage sub-routines relevant at higher northerly latitudes (Okkonen, 2011).However, computational approaches can be employed to produce the values on spatial and temporal variability in recharge often needed in groundwater modelling (Dripps and Bradbury, 2010).The methods commonly rely on a geographic information system (GIS) platform for spatial representation and calculation approaches based on water balance to create the temporal dimension of recharge (Croteau et al., 2010;Dripps and Bradbury, 2007;Jyrkama et al., 2002;Sophocleous, 2000;Westenbroeck et al., 2010).Neglecting variations in thickness of the unsaturated zone is common practice in many water balance models used in recharge estimations.However, the residence time in the unsaturated zone may play an important role, especially in the timing of recharge in deep unsaturated zones (Hunt et al., 2008), as has been acknowledged in recent work (Assefa and Woodbury, 2013;Jyrkama and Sykes, 2007;Scibek and Allen, 2006;Smerdon et al., 2008).
In numerical recharge models, actual evapotranspiration (ET) is a difficult variable to estimate accurately from climate, soil, and land use data.The vegetation is commonly parameterized from land use or land cover maps (Assefa and Woodbury, 2013;Jyrkama et al., 2002;Jyrkama and Sykes, 2007;Keese et al., 2005), where the vegetation characteristics and leaf area index (LAI) are estimated based solely on vegetation type.In addition to tree canopy transpiration, soil evaporation (SE), i.e. evaporation from the pores of soil matrix, can constitute a large proportion of total ET.Soil evaporation from the forest floor is generally reported to range from 3 to 40 % of total ET (Kelliher et al., 1993), although values as high as 92 % have been recorded (Kelliher et al., 1998).For conifer forest canopies, SE can largely compensate for low transpiration in areas with lower LAI (Ohta et  A simulated sediment profile is shown to give an example of how 1-D simulations are represented in the model domain and UZT represents the unsaturated zone thickness parameter.al., 2001;Vesala et al., 2005).Data on canopy-scale evaporation rates at latitudes above 60 • N are rare (Kelliher et al., 1993).A few studies have estimated ET from pine tree stands at patch scale (Kelliher et al., 1998;Lindroth, 1985), but none has extended this analysis to spatially distributed groundwater recharge.Forest management practices have the potential to affect the transpiration characteristics of coniferous forests, which typically leads to increased groundwater recharge (Bent, 2001;Lagergren et al., 2008;Rothacher, 1970).
The overall aim of the study was to provide novel information on groundwater recharge rates and factors contributing to the amount, timing, and uncertainty of groundwater recharge in unconfined sandy eskers aquifers.This study expands the application of physically based 1-D unsaturated water-flow modelling for groundwater recharge while taking into account detailed information on vegetation (pine, lichen), unsaturated layer thickness, cold climate, and simulation parameter uncertainty.Furthermore, this study considers the effect that forestry land use has on vegetation parameters and how this is reflected in groundwater recharge.

Study site
Groundwater recharge was estimated for the case of the Rokua esker aquifer in northern Finland (Fig. 1).Rokua is an unconfined aquifer consisting of unconsolidated sandy sediments underlain by crystalline bedrock (Fig. 2).Aquifer was formed during previous deglaciation when rivers under the melting ice sheet deposited sandy sediments in the river bed (Aartolahti, 1973).The Rokua esker has a rolling surface topography in the aquifer recharge area rising about 60 m above the flat peatland areas surrounding the esker.In the groundwater discharge areas, the aquifer is locally confined by peat soil with low hydraulic conductivity (Rossi et al., 2012).The climate at the Rokua aquifer is characterized by precipitation exceeding ET on an annual basis and statistics of the annual climate for the study period 1961-2010 in terms of precipitation, air temperature and FAO reference ET according to Allen et al. (1998) is presented in Table 1.Another important feature of the climate is annually recurring winter periods when most precipitation is accumulated as snow.

Leaf area index from forestry inventories
Forestry inventory data from the Finnish Forest Administration (Metsähallitus, MH) and Finnish Forest Centre (Metsäkeskus, MK) were used to estimate LAI for the Rokua esker groundwater recharge area.The available data consisted of 2786 individual plots covering an area of 52.4 km 2 (62.4 % of the model domain).The forestry inventories, performed mainly during 2000-2011, showed that Scots pine (Pinus sylvestris) is the dominant tree in the model area (94.2 % of plots).The forest inventory data include a number of data attributes and the following data fields, included in both the MH and MK data sets, were used in the analysis: -plot area (p A ; ha); -main canopy type; -average tree stand height (h; m); -average stand diameter at breast height (d bh ; cm); -number of stems (n stm ; 1 ha −1 ); -stand base area (b A ; m 2 ha −1 ); -stand total volume (V; m 3 ).
Inventory plots were excluded from the analysis if (1) main canopy type was not pine forest, (2) data were missing for d bh and h or n stm , or (3) the MH and MK data sets overlapped, in which case MH was retained.However, several plots in the MH data set were lacking n stm data, which would have created a large gap in data coverage.Therefore, the n stm variable was estimated with a log-transformed regression equation using data on d bh , p A , and V as independent variables.This regression equation was built from 280 plots (R 2 = 0.88) and used to estimate n stm for 288 plots.LAI was estimated as described by Koivusalo et al. (2008).Needle mass for an average tree in stand/plot was estimated from h and d bh using empirical equations presented by Repola et al. (2007).LAI for a stand was calculated as where N t is needle mass per average tree in stand (kg), n stm is number of stems per hectare (1 ha −1 ), and S LA is specific leaf area = 4.43 m 2 kg −1 = 4.43 × 10 −4 ha kg −1 (Xiao et al., 2006).Detailed information on LAI was used to obtain an estimate of how different forest management options, already actively in operation in the area, could potentially affect groundwater recharge.Three scenarios were simulated testing the potential impact of forestry operations on groundwater recharge: 1.The first baseline scenario simulated the current situation by using a LAI pattern at the site (Fig. 3) estimated with Eq. (1).
2. The second scenario simulated the impact of intensive forestry operations as clear-cutting of the tree stand.Clear-cutting is an intensive land use form in which almost the entire tree stand is removed, and it has been carried out in some parts of the study area.Low LAI values of 0-0.2 for the whole study site were used in simulating this scenario.
3. The third scenario simulated the impact of no forestry operations, i.e. absence of forestry cuttings.The hypothetical mature stand covering the study site was assumed to have high LAI values of 3.2-3.5 found at the study site and reported in the literature (Koivusalo et al., 2008;Rautiainen et al., 2012;Vincke and Thiry, 2008;Wang et al., 2004).

Lichen water retention
An organic lichen layer covers much of the sandy soil at the Rokua study site (Kumpula et al., 2000), so this lichen layer was introduced in SE calculations.Although lichens do not transpire water, their structural properties allow water storage in the lichen matrix and capillary water uptake from the soil (Blum, 1973;Larson, 1979).In this study lichen layer was explicitly included in the simulations to create an additional storage for water before the mineral sandy sediments.Water interception storage by the lichen layer was estimated from lichen samples.In total, six samples (species Cladonia stellaris and C. rangiferina) were taken in May 2011 from two locations 500 m apart, close to borehole MEA506 (see Fig. 1).These samples were collected by pressing plastic cylinders (diameter 10.6 cm) through the lichen layer and extracting intact cores, after which mineral soil was carefully removed from the base of the sample.Thus, the final sample consisted of a lichen layer on top and a layer of organic litter and decomposed lichen at the bottom, and was sealed in a plastic bag for transportation.To obtain estimates of water retention capacity, the samples were first wetted until saturation with a sprinkler, left overnight at +4 • C to allow gravitational drainage and weighed to determine field capacity.The samples were then allowed to dry at room temperature and weighed daily until stable final weight (dry weight) was reached.The water retention capacity (w r ) of the sample was calculated as where m f c is the field capacity weight (M), m dry is the final dry weight (M) at room temperature, ρ w (M L −3 ) is the density of water, and r (L) is the radius of the sampling cylinder.
The mean water retention capacity of the lichen samples was found to be 9.85 mm (standard deviation (SD) 2.71 mm) and approximations for these values were used in model parameterization (Table 2).Measured lichen water retention capacity was introduced to the simulations using parameters for soil porosity and layer thickness.Lichen porosity values varied between 7.5 and 12.5 % in simulation Monte Carlo runs (see Sect. 2.2.1) while keeping thickness of the lichen layer at 100 mm.In this manner the maximum amount of water retained by the lichen layer after gravitational drainage was adjusted to vary between 7.5 and 12.5 mm, as seen in the measurement data.To acknowledge the lack of information about Brooks and Corey (B&C) parameter estimates for lichen, the parameters were included in the simulation Monte Carlo runs with ranges which in our opinion produced a reasonable shape of the pressure-saturation curve allowing easy drainage of the lichen.

Sandy sediment hydraulic properties
Sediment texture was determined by sieving (ISO 3310-1 standard sieve, US sieve numbers 5,10,18,35,60,120 ) 26 samples taken from five boreholes at various depths (Fig. 1).Fourteen of the samples were also analysed for pressure saturation curves.Samples were characterized as fine or medium sand, while sediment texture in the other boreholes (Fig. 1) had previously been characterized as medium, fine or silty sand throughout the model domain by the Finnish Environmental Administration as expert in situ analysis during borehole drilling.Therefore, the sediment samples from the five boreholes were considered to be representative of the sediment type in the area.Pressure saturation data from the samples were then used to define parameter ranges for the B&C equation used in the simulations (Table 2).Furthermore, texture values were employed to calculate the range of saturated vertical hydraulic conductivity for the samples, using empirical equations by Hazen, Kozeny-Carman, Breyer, Slitcher, and Terzaghi (Odong, 2007).The hydraulic conductivity for a given sample ranged at approximately 1 order of magnitude between the equations.When using the five equations for the 26 samples in total, the calculated values were within 1.99 × 10 −5 -1.47 × 10 −3 (m s −1 ) for all but one sample.The obtained range was considered to reasonably represent the hydraulic conductivity variability in the study area and simulations.

Unsaturated zone thickness
The thickness of the unsaturated zone at each model cell was estimated by subtracting the interpolated water Kettle hole lakes in the area are imbedded in the aquifer and thus reflect the level of the regional water table (Alaaho et al., 2013).The lake stage was extracted as the DEM elevation for a given lake, while for large lakes several interpolation points were scattered around the lake shore to better steer the interpolation locally.
Wetland elevation was used as a proxy for the water table elevation in locations where more certain observations (piezometers, lake levels) were lacking.If a wetland was present in the topographical depression, the water table was considered to lie at the depression bottom, in order to sustain the conditions needed for wetland formation.Wetlands were detected from the base map and the value for the water table proxy was assigned from the DEM.
Finally, the land surface elevation was considered to give a reasonable estimate of the water table position in the transition zone between aquifer recharge area (model domain) and groundwater discharge areas covered by peatlands (Fig. 2).The Rokua aquifer is phreatic in the recharge area and Rossi et al. (2012) demonstrated that the peatlands partially confine the aquifer and can create artesian conditions in the discharge area.Even though some local overestimation of the water table may have resulted from the approximation method at the transition zone, it was found to be important to have some points to guide the interpolation at the model domain boundary in order to acknowledge the characteristics of the sloping water table towards the discharge area (Fig. 2).The proxy used for the water table was extracted from the DEM to points approximately 250 m apart at the boundary of the model domain.

Climate data
Driving climate data for the simulations were taken from Finnish Meteorological Institute databases for the modelling period 1 January 1961-31 October 2010.Daily mean temperature ( • C) and sum of precipitation (mm) were recorded at Pelso climate station, 6 km south of the study area (Fig. 1).The most representative long-term global radiation data (kJ m −2 d −1 ) for the area were available as interpolated values in a 10 km × 10 km grid covering the entire area of Finland.The interpolation data point was found to be at approximately the same location as the borehole MEA2110 (Fig. 1).Long-term data on wind speed (m s −1 ) and relative humidity (%) were taken from Oulunsalo and Kajaani airports, located 60 and 40 km from the study site, respectively.The data from the airports were instantaneous observations at 3 h intervals, from which daily mean values were calculated.All the climate variables were recorded at the reference height of 2 m except for wind speed, which was measured at 10 m height.The wind speed data were therefore recalculated to correspond to 2 m measurement height according to Allen et al. (1998) by multiplying daily average wind speed by 0.748.The suitability of long-term climate data for the study site conditions was verified with observations made at a climate station established at the study site in an overlapping time period (December 2009-October 2010), and the agreement between the measurements was found to be satisfactory.
Data on long-term lake surface-water temperature were needed to calculate lake evaporation (see Sect. 2.2.3), but were not available directly at the study site.However, surface-water temperature was recorded at Lake Oulujärvi by the Finnish environmental administration (2013) 22 km from the study site in the direction of the Kajaani climate station (Fig. 1).The Oulujärvi water temperature was found to be closely correlated (linear correlation coefficient 0.97) with daily lake water temperature recorded at Rokua during summer 2012.Daily lake surface temperature data for Lake Oulujärvi starting from 21 July 1970 were used in lake evaporation modelling.However, the data series had missing values for early spring and some gaps during 5 years in the observation period.These missing values were estimated with a sine function, corresponding to the average annual lake temperature cycle, and a daily time series was established for subsequent calculations.
Snowmelt was calculated with a degree-day approach model in Jansson and Karlberg (2004).Snow routines were calibrated separately using bi-weekly snow water equivalent (SWE) data from Vaala snow line measurements (Finnish environment administration, 2011) for the period 1960-2010 (Fig. 1).This separately calibrated snow model was used for all subsequent simulations.

Water-flow simulation in 1-D unsaturated sediment profile
Water flow through an unsaturated one-dimensional (1-D) sandy sediment profile (Fig. 2) was estimated with the Richards' equation using the CoupModel (Jansson and Karlberg, 2004).The CoupModel was selected as the simulation code because of its ability to represent the full soil-plantatmosphere continuum adequately and to include snow processes in the simulations (Okkonen and Kløve, 2011).The simulated sediment profile was vertically discretized into 61 layers with increasing layer thickness deeper in the profile.Layer thickness was 0.1 m for the first 16 layers (until 1.6 m), where the topmost 0.1 m was represented as a lichen layer.Layer thickness was progressively increased by defining 0.2 m thickness for the next 7 layers (between 1.6 and 3 m), 0.5 m for the next 14 layers (between 3 and 10 m), 1 m for the next 7 layers (between 10 and 17 m) and 2 m for last 17 layers ranging from 17 m to the bottom of the profile (51 m).The time variable boundary condition for water flow at the top of the column was defined by driving climate variables and affected by sub-routines accounting for snow processes with daily time step.The short time step was chosen to fully capture the main recharge input from snowmelt.All water at the top of the domain was assumed to be subjected to infiltration.Deep percolation as gravitational drainage was allowed from sediment column base using the unit-gradient boundary condition (see e.g.Scanlon et al., 2002b).Simulations for the unsaturated 1-D sediment profile were made for the period 1970-2010, and before each run 10 years of data (1960)(1961)(1962)(1963)(1964)(1965)(1966)(1967)(1968)(1969)(1970) were used to spin up the model.2. For each individual simulation homogeneity in the vertical direction in terms of sediment hydraulic properties was assumed.The parameters for which values were randomly varied were chosen beforehand by trial and error model runs exploring the sensitivity of parameters with respect to cumulative recharge or ET.The parameter ranges were specified from field data when possible; otherwise, we resorted to literature estimates or in some cases used ±50 % of the CoupModel default providing a typical parameter for the used equation.
The sensitivity of the parameters varied in the simulations was tested with Kendall correlation analysis, by testing the correlation between each model parameter and cumulative sums of different ET components and infiltration for the 400 model runs.Individual simulation with unique parameter values did not produce a groundwater recharge value due to the assembling strategy for recharge; therefore, the ET components and infiltration were selected as variables for comparison.In addition, correlations were examined as scatter plots to ensure that possible sensitivity not captured by the monotonic correlation coefficient was not overlooked.

Method to distribute 1-D simulations spatially
Groundwater recharge was estimated for a model domain of 82.3 km 2 (Fig. 1).To distribute the simulations in 1-D sediment column spatially, the simulation domain was subdivided into different recharge zones, similarly to, e.g.Jyrkämä et al. (2002).Zonation in the model was based on two variables: LAI and unsaturated zone thickness (UZT).The calculation of spatially distributed values for LAI and UZT is presented in detail in Sect.Variation in the LAI and UZT parameters was used to allocate the 1-D sediment profile simulations spatially to the study site.LAI class in model cell specified a subset of the 400 1-D simulations that were applicable for a given cell.UZT class for each cell (Fig. 2) specified the depth in the simulated 51 m sediment profile where the water flux output was extracted.Using this approach each unique recharge zone (a combination of UZT and LAI class) had on average 27 water-flow time series (number of total model runs -400 divided by number of LAI cell classes -15) produced by different random combinations of parameters (Table 2).Equation (3) was used to propagate the variability in the 27 time series into the final areal recharge.
where R i,j is the final sample of areal recharge (mm day −1 ), i is the index for simulation time step (e.g. 1 : 14 975), j is the index for sample for a given time step (1 : 150), l is the index for unique recharge zone, n(l) is the number of cells in a given recharge zone, Rs is the recharge sample (mm d −1 ) for a given recharge zone at time step i, k is the number of time series for a given recharge zone, A c is the surface area of a model raster cell (e.g.20 m × 20 m = 400 m 2 ), and A tot is the surface area of the total recharge area.
The resulting R matrix has 150 time series for areal recharge produced by simulations with different parameter realizations.The variability between the time series provides an indication of how much the simulated recharge varies due to different model parameter values.The method allows for computationally efficient recharge simulations, because the different recharge zones do not all have to be simulated separately.
The simulation approach assumes the following: (1) over the long-term, the water table remains at a constant level, i.e. the unsaturated thickness for each model cells stays the same.Monitoring data from 11 boreholes and seven lakes with more than 5 years of observation history shows a level variability of 1-1.5 m, with depressions and recoveries of the water table.This variability is within the accuracy of water table estimation by interpolation.(2) The capillary fringe in the sandy sediment is thin enough not to affect the water flow before arriving at the imaginary water table at the centre of each UZT class.(3) Only vertical flow takes place in the unsaturated sediment matrix, a typical assumption in recharge estimation techniques (Dripps and Bradbury, 2010;Jyrkama et al., 2002;Scanlon et al., 2002a).(4) Surface runoff is negligible primarily due to the permeable sediment type (as noted by Keese et al., 2005), and also due to the lichen cover inhibiting runoff by increasing surface roughness (Rodríguez-Caballero et al., 2012).The maximum observed daily rainfall for the area has been 57.4 mm.Further assuming that rain for the day fell only during 1 h, it would be equal to 1.59 × 10 −5 m s −1 input rate of water, which is close to the lower range of saturated hydraulic conductivity at the study site (1.99 × 10 −5 m s −1 ).Therefore, rainstorms at the site very rarely exceed the theoretical infiltration capacity.Being a field verification, surface runoff has not been observed during field visits and the area lacks intermittent or ephemeral stream networks.

Estimation of evapotranspiration
Four different evaporation processes were considered in this study (Fig. 5): SE (evaporation from the topmost soil layer, i.e. the lichen matrix), snow evaporation (evaporation from snow surface), transpiration (evaporation through the vascular system of tree canopy), and lake evaporation (evaporation from a free water surface).The first three components were estimated, along with water-flow simulations, using Coup-Model.However, as 3.6 % (2.9 km 2 ) of the surface area of the study site consists of lakes (Fig. 1), lake evaporation from P. Ala-aho et al.: Temporal and spatial recharge estimation free water surfaces was calculated independently from the CoupModel simulations.
Transpiration from the Scots pine canopy (L v E tp ) was calculated using the Penman-Monteith (P-M) combination (Eq.A1).Whenever possible, all the parameters relating to the P-M equation were estimated based on data, namely LAI of the canopy.Surface resistance and saturation vapour pressure difference are the main factors controlling conifer forest ET, while the aerodynamic resistance is of less importance (Lindroth, 1985;Ohta et al., 2001).In the calculation of aerodynamic resistance with the P-M equation, roughness length is related to LAI and canopy height, according to Shaw and Pereira (1982).Other parameters governing the aerodynamic resistance, except for LAI, were treated as constant.The surface resistance of the pine canopy was estimated with the Lohammar equation (see e.g.Lindroth, 1985), accounting for effects of solar radiation and air moisture deficit in tree canopy gas exchange.Because LAI values have a strong influence in the surface resistance Lohammar equation, the other parameters governing the surface resistance were excluded from the Monte Carlo runs.Distribution of root biomass with respect to depth from the soil surface was presented with an exponential function, because most Scots pine roots are concentrated in the shallow soil zone.A typical root depth value of 1 m was used for the entire canopy (Kalliokoski, 2011;Kelliher et al., 1998;Vincke and Thiry, 2008).Soil and snow evaporation were calculated using an empirical approach (Eq.A4) based on the P-M equation, as described in detail in Jansson and Karlberg (2004).Soil evaporation is calculated for the snow-free fraction of the soil surface, and the snow evaporation is solved separately as a part of snowpack water balance.
In areas where the water table is close to the ground surface, the water table can provide an additional source of water for ET (Smerdon et al., 2008).To take into account the decreased recharge for areas with near-surface-water tables, the recharge for cells with an unsaturated zone of < 1 m (8.3 % of the study site, 6.8 km 2 ) was estimated with a water balance approach.We assumed that for areas with a shallow water table, soil water content was not a limiting factor for transpiration.Therefore, an additional water source for transpiration was considered by making the transpiration rate equal to simulated potential transpiration (T ) during times when the actual transpiration was simulated (T > 0.05 mm).The increasing effect of the water table located at 1 m depth on SE was tested with simulations and found to be 5-10 % higher with than without a water table.Therefore, a 7 % addition was made to the simulated actual SE for cells with a shallow water table.Daily recharge (R 1 m , L T −1 ) for cells with unsaturated thickness below 1 m was estimated as where I is infiltration water arriving to lake/soil surface, including both meltwater from the snowpack and precipitation (L T −1 ), T adj (mm d −1 ) is adjusted transpiration, and ES adj (mm d −1 ) is adjusted SE.Kettle hole lakes in esker aquifers often lack surface-water inlets and outlets and are therefore an integral part of the groundwater system (Ala-aho et al., 2013;Winter et al., 1998), so we considered these lakes as contributors to total groundwater recharge.In other words, rainfall per lake surface area is treated equally as an addition to the aquifer water storage as groundwater recharge.As a difference, lake water table is subjected to evaporation unlike the groundwater table.Lake evaporation (E lake ) was estimated with the mass transfer approach (see e.g.Dingman, 2008) according to Eq. (A7).The mass transfer method was selected because of its simplicity, daily output resolution, low data requirement, and physically based approach.However, various calculation methods could easily be used in the modelling framework, depending on the data availability (see e.g.Rosenberry et al., 2007).If lake percentage in the area of interest is high, more sophisticated methods may be required to better represent the system.For the Rokua site the bias introduced by a simplistic approach was considered minor.

Model validation
Model performance was tested by comparing the simulated recharge values with two independent recharge estimates on a local and regional scale: the water table fluctuation (WTF) method and baseflow estimation, respectively.The WTF method is routinely used to estimate groundwater recharge because of its simplicity and ease of use.It assumes that any rise in water level in an unconfined aquifer is caused by recharge arriving at the water table.For a detailed description of the method and its limitations, see e.g.Healy and Cook (2002).The recharge amount (R, L T −1 ) is calculated based on the water level prior to and after the recharge event and the specific yield of the sandy sediments: where S y is the specific yield, h is the water table height (L), and t is the time of water table rise (T).
The WTF method requires groundwater level data with adequate resolution for both time and water level, to identify periods of rising and falling water table.The water table was monitored using pressure-based data loggers (Solinst Levelogger Gold) recording at hourly interval from six water table wells with average unsaturated zone thicknesses of 1.2, 1.6, 5.0, 8.0, 9.3, and 14.7 m (Fig. 1).Wells where the water table was < 2 m from the ground surface responded to major precipitation events.In the deeper wells, only the recharge from snowmelt was seen as water table rise.Estimates of the sandy sediment-specific yield are required for the calculations (Eq.5), but no sediment samples were available from the wells used in the water table monitoring.Drilling records for these wells reported fine and medium sand, which was consistent with the particle size distribution for other wells in the area.Therefore, an estimated value of 0.20-0.25 for the specific yield of all wells was used, according to typical values for fine and medium sand (Johnson, 1967).
The recharge estimated with the WTF method was compared with the simulated recharge during the recorded water level rise in the well.For each well, the cumulative sum of simulated water flow was extracted from sediment profile depth corresponding to well water table depth.As an example, the simulated recharge in well ROK1 (unsaturated thickness on average 14.7 m) was extracted from the UZT class 12, corresponding to recharge for unsaturated thickness of 14-16 m.All 400 model runs were used, providing 400 estimates for recharge for each time period of recorded water level rise.
A regional estimate of groundwater recharge was estimated as baseflow of streams originating at the groundwater discharge area.Because the Rokua esker aquifer acts as a regional water divide, stream flow was monitored around the esker, in a total of 18 locations (Fig. 1).The flows were measured a total of eight times between 6 July 2009 and 3 August 2010 (see Rossi et al., 2014).The lowest total outflow during 9-10 February 2010 was recorded after 3 months of snow cover period, when water contribution to streams from surface runoff was minimal.The minimum outflow was considered as baseflow from the aquifer reflecting long-term groundwater recharge in the area.

Model validation with the WTF and baseflow methods
Model validation showed that the modelling approach could reasonably reproduce (1) the main groundwater recharge events when compared to the WTF method (Fig. 6) and (2) the regional level of recharge compared to stream baseflow.
The WTF method agreed well with the simulated values, with overlapping estimates between the methods for all but two recharge events.Also the median value of simulations was close to the WTF method, with some bias to higher estimates from the simulations.The discrepancy can be due to very different assumptions behind the methods and uncertainty in local parameterization: in the WTF method for the specific yield (S y ) and for simulations mainly the hydraulic conductivity which dictates the simulated timing of recharge.Uncertainty in the S y estimate is acknowledged by showing S y as a range rather than a single value (Fig. 6), but still S y is not truly known for the location of observation boreholes.Simplifying assumptions and subjective interpretation of both timing and height of water table rise create additional inaccuracies in the WTF estimate.Independent regional estimate of recharge, stream baseflow, was 70 500 m 3 s −1 , or 312.7 mm a −1 when related to the recharge area.The order of magnitude agreed with the long-term simulated average of 362.8 mm a −1 .Typical error in individual stream flow measurements is within 3-6 % of the measured value (Sauer and Meyer, 1992), which brings minor uncertainty in the baseflow value.The smaller value for stream baseflow compared to simulated long-term average recharge can be explained with conceptual understanding of site hydrogeology (Ala-aho et al., 2013(Ala-aho et al., , 2015;;Rossi et al., 2012Rossi et al., , 2014)).Part of the recharged groundwater does not discharge to the small streams whose baseflow was measured, but flows underneath the stream catchments and seeps out to regional surface bodies (Oulujärvi and Oulujoki) further away from the recharge area (Fig. 2).Fully integrated surface-subsurface hydrological modelling study of the same site presented in Ala-aho et al. ( 2015) simulated an outflow of 79 mm a −1 to regional surface-water bodies.

Temporal variations in groundwater recharge
When recharge simulation time series were summarized to annual values (1 October-30 September), recharge rates covaried with annual infiltration with linear correlation coefficient of 0.89 (Fig. 7) as expected based on previous work in humid climate and sandy aquifers (Keese et al., 2005;Lemmelä, 1990).Both annual recharge and infiltration displayed an increasing trend.The plot also showed the level of uncertainty in annual recharge values introduced by differences in model parameterization (black area).The difference between minimum and maximum value for simulated annual recharge was on average 23.0 mm.Thus, the variability in recharge estimates was 6.3 % of mean annual recharge 362.8 mm.
According to the simulations, the effective rainfall, i.e. the percentage of corrected rainfall resulting in groundwater recharge annually, was on average 59.3 %.This is in agreement with previous studies on unconfined esker aquifers at northerly latitudes, in which the proportion of annual precipitation percolating to recharge is reported to be 50-70 % (71 % by Zaitsoff (1984), 54 % by Lemmelä and Tattari (1986) and 56 % by Lemmelä, 1990).The percentage of effective rainfall varied considerably, by almost 30 %, between different hydrological years, from 44.8 % in some years up to 73.1 % in others.

Influence of LAI on spatial variation of groundwater recharge
The spatial distribution of groundwater recharge was mostly due to variations in LAI, as well as influenced by distance to water table, and distribution of lakes (Fig. 8).Higher evaporation rates from lakes led to lower recharge in lakes (see red spots in Fig. 8).Similarly, high LAI led to high ET and resulted in low recharge in plots with high LAI.Other areas of low recharge, although not as obvious at the larger spatial scales shown in Fig. 8, were cells with a shallow water table (Sect.2.2.3).The effect of high ET at locations with a shallow water table can best be seen in south-east parts of the aquifer.Kendall correlation analysis of simulation parameters and annual average model outputs identified LAI as the most important parameter controlling ET and infiltration (Table 3).Parameters related to sediment hydraulics and evaporation  showed some sensitivity to simulation results, while the parameters for lichen vegetation were only slightly sensitive or insensitive to simulation output variables.The LAI parameter governed the level of evaporation for different ET components (Fig. 9).Evaporation from soil (and snow) compensated for mean annual ET for LAI values up to around 1.0, after which total ET increased as a function of LAI.
The scenarios for low (0-0.2) and high (3.2-3.5)LAI changed the groundwater recharge rates compared to the current LAI distribution (in Fig. 7).In the high LAI scenario the annual recharge was on average 101.7 mm lower than in the low LAI scenario.These results suggest that management of the Scots pine canopy has a significant control on the total recharge rates in unconfined esker aquifers.Average land surface ET components remained relatively constant between years, but the simulated ET displayed a wide difference between simulations (Fig. 10).Estimated annual ET (mean 237.6 mm) was somewhat lower than previous regional estimates of total ET (300 mm; Mustonen, 1986).Lake evaporation rates were generally higher than ET from the land surface (420.0 mm).The variation in simulated lake evaporation was considerably lower than that in ET, as a different approach was used to account for uncertainty in the simulations.Transpiration showed greater variation between simulations than SE and total land surface ET.On average, transpiration also comprised a slightly larger share of total evaporation than SE.Simulated snow evaporation was a small, yet not insignificant, component in the total ET.

Discussions
The method used here to estimate LAI from forestry inventories introduces a new approach for incorporating large spatial coverage of detailed conifer canopy data into groundwater recharge estimations.LAI values reported for conifer forests in Nordic conditions similar to the study site are in the range 1-3, depending on canopy density and other attributes (Koivusalo et al., 2008;Rautiainen et al., 2012;Vincke and Thiry, 2008;Wang et al., 2004).The LAI values obtained for the study site (mean 1.25) were at the lower end of this range.Furthermore, the data showed a bimodal distribution, with many model cells with low LAI (< 0.4) lowering the mean LAI.The low LAI values were expected because of active logging and clear-cutting activities in the study area.Although the equations to estimate LAI are empirical in nature and based on simplified assumptions, the method can outline spatial differences in canopy structure.Wider use of this method in Finland is practically possible, as active forestry operations in Finland have yielded an extensive database on canopy coverage, which could be used in groundwater management.However, the LAI estimation method could be further validated with field measurements or lidar techniques (Chasmer et al., 2012;Riaño et al., 2004).
Plant cover, represented as LAI, proved to be the most important model parameter controlling total ET, and thereby the amount of groundwater recharge (Table 3, Fig. 9).The LAI parameter was included in the equations controlling both transpiration and SE, and therefore the sensitivity of the pa-P.Ala-aho et al.: Temporal and spatial recharge estimation rameter is not surprising.While SE partly compensated for the lower transpiration with low LAI values, the total annual ET values progressively increased as a function of LAI (Fig. 9).Interestingly, the simulations suggested that ET remains at a constant level in the LAI range 0-1, potentially due to the sparse canopy changing the aerodynamic resistance and partitioning of radiation limiting SE, while still not contributing much to transpiration in total ET.This implies that the maximum groundwater recharge for boreal Scots pine remains rather constant up to a threshold LAI value of around 1.This knowledge can be used when co-managing forest and groundwater resources in order to optimize both.
Importance of LAI has been reported in earlier studies estimating groundwater recharge (Dripps, 2012;Keese et al., 2005;Sophocleous, 2000), but here the vegetation was represented with more spatially detailed patterns and a field databased approach for LAI.According to previous studies, average ET from boreal conifer forests is around 2 mm d −1 during the growing season (Kelliher et al., 1998), which is similar to our average value of 1.6 mm d −1 for the period 1 May-31 October.Some earlier studies have claimed that the influence of LAI on total ET rates from boreal conifer canopies is minor (Kelliher et al., 1993;Ohta et al., 2001;Vesala et al., 2005), but our simulation results indicate that higher LAI values lead to higher total ET values.The simulations showed that variable intensity of forestry, from low canopy coverage (LAI = 0-0.2) to dense coverage (LAI = 3.2-3.5)resulted in a difference of over 100 mm in annual recharge (Fig. 7).It can be argued that the scenarios are unrealistic, because high LAI values, covering the whole study site, may not be achieved even with a complete absence of forestry operations.Nevertheless, the result demonstrates a substantial impact of forestry operations on esker aquifer groundwater resources.
The lichen layer covering the soil surface was explicitly accounted for in the simulation set-up, which to our knowledge is a novel modification.Kelliher et al. (1998) concluded that precipitation intercepted by lichen was an important source of understorey evaporation, especially directly after rain events.In addition, Bello and Arama (1989) reported that lichen could intercept light rain showers completely and that only intense rain events caused drainage from lichen canopy to mineral soil.While the lichen layer might have an increasing effect on SE through interception storage, Fitzjarrald and Moore (1992) suggested that a lichen cover may in fact have an insulating influence on heat and vapour exchange between soil and atmosphere, therefore impeding evaporation from the mineral soil.In the present study, the lichen layer appeared to have a minor influence on total evaporation, SE, and infiltration, as these variables showed only little sensitivity to lichen B&C parameters (Table 3).However, the approach to represent lichen with the B&C model needs to be better examined, as water retention capacity of the lichen layer was introduced to the simulations using the concept of total porosity, which is not strictly coherent with the B&C model.Nevertheless, the used approach successfully produced an additional dynamic interception storage of water in the correct range (generally 3-7 mm depending on random parameterization, data not shown).The performed laboratory measurement of lichen water retention should be supplemented with detailed analysis of lichen pressuresaturation curve and hydraulic conductivity to clarify the role of lichen in SE, and thereby groundwater recharge.
Stochastic variation of selected model parameters illustrated the uncertainties relating to numerical recharge estimation using the Richards' equation in one dimension.The capability and robustness of the Richards' equation to reproduce soil water content and water fluxes have been demonstrated extensively in various studies (Assefa and Woodbury, 2013;Scanlon et al., 2002b;Stähli et al., 1999;Wierenga et al., 1991).Therefore, we considered that model calibration and validation with point observations of variables such as soil volumetric water content or soil temperature would not provide novel insights into water flow in unsaturated porous media.Instead, we incorporated the parameter uncertainty ranges, usually used in model calibration, to the final recharge simulation output.An important outcome was that the uncertainty in the model output caused by different model parameterizations was small in comparison with the interannual variation in recharge.The error caused by uncertainty in the model assumptions or driving climate data was not addressed in this study.
The sensitivity analysis focused on total cumulative values of fluxes and did not address the temporal variations in the variables.Sediment hydraulic parameters mainly influenced the timing of recharge through residence time in the unsaturated zone, not so much the total amount.Therefore, the sediment hydraulic parameters showed only minor sensitivity, perhaps misleadingly.It should be noted that vertical heterogeneity in the unsaturated sediment profile hydraulic parameters can reduce the total recharge rates (Keese et al., 2005).However, vertical heterogeneities were ignored in this study in order to simplify the model, and because the drilling logs showed only little variation in the area.Work by Wierenga et al. (1991) supports the simplification by showing that excluding moderate vertical heterogeneities does not significantly affect the performance of water-flow simulations with the Richards' equation.
Simulations acknowledged shallow water table contribution to ET in an indirect, conceptual approach.Including a water table fixed at different depths in the sediment profile would have been possible in the CoupModel set-up.Influence of the water table fixed at 2 m depth was tested and found to increase ET 3.5 % for LAI values of 3, but for LAI values of 0.5 and 1.5 the increase in ET was only trivial.We expect only a minor increase in ET with deeper water table configuration (with the given sediment texture), and therefore argue that excluding the water table results in only minimal overestimation of total recharge at the study site.It should be noted that upward water fluxes were not excluded from the water-flow time series and negative fluxes were considered as negative recharge at any depth.The simplification is made that water available for upward fluxes comes only from the soil moisture storage, not from the water table.
According to the simulations, the percentage of precipitation forming groundwater recharge varied considerably between years, as also reported in previous studies on transient recharge (Assefa and Woodbury, 2013;Dripps and Bradbury, 2010).Even though annual recharge was correlated with annual precipitation, and therefore years with high precipitation resulted in higher absolute recharge (Fig. 7), the percentage of effective rainfall did not increase as a function of annual sum of precipitation.This is somewhat surprising, because the rather constant evaporation potential between years (Fig. 10) and high sediment hydraulic conductivity could be expected to result in a higher percentage of rainfall reaching the water table in rainy years.Some studies (Dripps and Bradbury, 2010;Okkonen and Kløve, 2010) have suggested that when the main annual water input arrives as snowmelt during the low evaporation season, it is likely to result in a higher percentage recharge than in a year with little snow storage and precipitation distributed evenly throughout summer and autumn, which may contribute to the variability in the effective rainfall coefficient.However, when the maximum annual SWE value was used as a proxy for annual snow storage, there was no evidence of snow amount explaining the interannual variability in the recharge coefficient.Other factors contributing to recharge coefficient variability may be related to soil moisture conditions prior to snowfall, or the intensity of summer precipitation events (Smerdon et al., 2008;Stähli et al., 1999).
The above-mentioned reasons make the concept of effective rainfall, which is currently routinely used to estimate groundwater recharge for groundwater management in, e.g.Finland (Britschgi et al., 2009), susceptible to over-or underestimation of actual annual recharge.This applies especially for aquifers with a thick unsaturated zone, where rainy years produce higher average recharge with some delay and for a longer duration (Zhou, 2009).

Conclusions
A physically based approach to simulate groundwater recharge for sandy unconfined aquifers in cold climates was developed.The method accounts for the influence of vegetation, unsaturated zone thickness, presence of lakes, and uncertainty in simulation parameters in the recharge estimate.It is capable of producing spatially and temporally distributed groundwater recharge values with uncertainty margins, which are generally lacking in recharge estimates, despite understanding of uncertainty related to recharge estimates being potentially crucial for groundwater resource management.However, the parameter uncertainty defined for the study area was of minor significance compared with inter-annual variations in the recharge rates introduced by climate variations.
The simulations showed that Scots pine canopy, parameterized as leaf area index (LAI), was important in controlling the total amount of groundwater recharge.Forestry inventory databases were used to estimate and spatially allocate the LAI and the results showed that such inventories could be better utilized in groundwater resource management.Forest cuttings were demonstrated to increase groundwater recharge significantly.A sensitivity analysis on the parameters used showed that SE could compensate for low LAI-related transpiration up to a LAI value of approximately 1, which may be important in finding the optimal level for forest management in groundwater resource areas.The concept of effective rainfall gave inconsistent estimates of recharge in annual timescales, showing the importance of using physically based recharge estimation methods for sustainable groundwater recharge management.
The Supplement related to this article is available online at doi:10.5194/hess-19-1961-2015-supplement.

Figure 1 .
Figure 1.Recharge area of the Rokua esker aquifer.Boreholes in the area were used for model validation and sediment type characterization.Baseflow was measured from streams originating outside the groundwater recharge area.Profile of cross section A-B is presented in Fig. 2.

Figure 2 .
Figure 2. Cross section A-B (Fig. 1) to demonstrate the geometry of the unsaturated zone and the aquifer (vertical axes exaggerated).A simulated sediment profile is shown to give an example of how 1-D simulations are represented in the model domain and UZT represents the unsaturated zone thickness parameter.

Figure 3 .
Figure 3. Spatial distribution of leaf area index (LAI) and a 20 m × 20 m cell-based histogram of LAI values.In areas where forestry inventory data were lacking, a weighted average value of 1.25 was used in simulations.

Figure 4 .
Figure 4.Estimated thickness of the unsaturated zone in the model area and interpolation points for estimation of water table elevation.
2.1.1 and 2.1.4.Both vari-Hydrol.Earth Syst.Sci., 19, 1961-1976, 2015 www.hydrol-earth-syst-sci.net/19/1961/2015/ ables were presented as a grid map with 20 m × 20 m cell size with a floating point number assigned to each cell, resulting in a total of 205 708 cells for the model domain.The small model cell size was selected to ensure full exploitation of the forest inventory plots in LAI determination.The spatially distributed data were then divided into 15 classes for LAI and 30 classes for UZT.The classes are primarily equal intervals, which was convenient in the subsequent data processing, but in addition the frequency distributions of LAI and UZT cell values were used to assign narrower classes for parameter ranges with many values (see histograms in Figs. 3 and 4).Class interval for LAI was 0.2 units up to a value of 2 (class 1: LAI = 0-0.2) and 0.3 to the maximum LAI value of 3.5.The class interval for UZT was 1 to 10 m depth and 2 m to the final depth of 51 m.Finally, the classified LAI and UZT data were combined to a raster map with 20 m × 20 m cell size, producing 449 different zones with unique combinations of LAI and UZT values.Spatial coupling was done with the ArcGIS software (ESRI, 2011).

Figure 5 .
Figure5.Flow chart of different evaporation processes considered in the study.Total evapotranspiration (ET) is comprised of soil evaporation from the topmost soil layer (i.e. the lichen matrix) snow evaporation from snow surface, transpiration through the vascular system of tree canopy and lake evaporation from free water surface.

Figure 6 .
Figure 6.Assemblage of simulated recharge for individual recharge events, shown as box plots where circles represent the median, bold lines 25-75th percentiles of the simulations, thin lines the remaining upper and lower 25th percentiles and plus symbols are outliers.The location of the box plots on the x axis is the WTF estimate for a given recharge event using a specific yield value of 0.225.The dashed lines indicate the uncertainty in the WTF estimates caused by the selection of specific yield.The two estimates would agree perfectly (given the uncertainty in S y ) if all simulations shown as box plots fell between the dashed lines.

Figure 7 .
Figure 7. Annual recharge time series from simulations where the black area covers the minimum and maximum values for different recharge samples.The annual recharge pattern closely followed trends in infiltration.Effects of different land use management practices over time on annual recharge rates are shown as high and low leaf area index (LAI) scenarios.

Figure 8 .
Figure 8. Spatial distribution of mean annual recharge, which was influenced mainly by the Scots pine canopy (LAI), the presence of lakes and, to some extent, areas with a shallow water table.

Figure 9 .
Figure 9. Example of scatter plots with the mean annual ET components are plotted as a function of the variable leaf area index (LAI), showing clear dependence of all ET components on LAI.

Figure 10 .
Figure 10.Values of different evapotranspiration (ET) components (mean and standard deviation) simulated for the study period.

Table 1 .
Characteristics of the study site annual climate.

Table 2 .
Jansson and Karlberg (2004) related equations and parameter ranges included in the model runs.For full description of parameters and equations, seeJansson and Karlberg (2004).

aho et al.: Temporal and spatial recharge estimation 1965 230
, and Hydrol.Earth Syst.Sci., 19, 1961-1976, 2015 www.hydrol-earth-syst-sci.net/19/1961/2015/ P. Ala- Water table borehole observations give the most accurate and reliable estimate of the water table position because they provide direct measurements on the water table.The water table elevation in a given piezometer was estimated here as the average value of the entire measurement history of each piezometer.
table level from the digital elevation model (DEM) topography calculated based on lidar data (National Land Survey of Finland, 2012).The water table elevation was estimated with the ordinary Kriging interpolation method from four types of observations (Fig. 4): water table boreholes (n = 19), stages of kettle hole lakes (n = 82), elevation of wetlands located in landscape depressions (n = 36), and land surface elevation at the model domain boundary (n = 229).

Table 3 .
Kendall correlation coefficient for simulation parameters and average annual sum of simulation output variables.ET is evapotranspiration, E is evaporation, for other symbols see Table2.