Modeling forest evapotranspiration and water balance at stand and catchment scales: a spatial approach

Vegetation is known to have strong influence on evapotranspiration (ET ), a major component of terrestrial water balance. Yet hydrological models often describe ET by methods unable to sufficiently include the variability of vegetation characteristics in their predictions. To take advantage of increasing availability of high-resolution open GIS-data on land use, vegetation and soil characteristics in the boreal zone, a modular, spatially distributed model for predicting ET and other hydrological 5 processes from a grid cell to a catchment level is presented and validated. An improved approach to upscale stomatal conductance to canopy scale using information on plant type (conifer / deciduous) and stand leaf-area index (LAI) is proposed by coupling a common leaf-scale stomatal conductance model with a simple canopy radiation transfer scheme. Further, a generic parametrization for vegetation-related hydrological processes for Nordic boreal forests is derived based on literature and data from a boreal FluxNet site. With the generic parametrization, the model was shown to well reproduce daily ET measured by 10 eddy-covariance technique at ten conifer-dominated Nordic forests whose LAI ranged from 0.2 to 6.8 m2m−2. Topography, soil and vegetation properties at 21 small boreal headwater catchments in Finland were derived from open GIS-data at 16 x 16 m grid size to upscale water balance from a stand to catchment level. The predictions of annualET and specific discharge were successful in all catchments, located from 60 to 68 ◦N, and daily discharge also reasonably well predicted by calibrating only one parameter against discharge data measurements. The role of vegetation heterogeneity on soil moisture and partitioning of 15 ET was demonstrated. The proposed framework can support e.g. forest trafficability forecasting and predicting the impacts of climate change and forest management on stand and catchment water balance. With appropriate parametrization it can be generalized outside the boreal coniferous forests. Copyright statement. Author(s) 2019. CC-BY 4.0 License


Introduction
The boreal region, encompassing ca. 12 % of world's land area, is characterized by mosaic of peatlands, lakes and forests of different ages and structures.The landscape heterogeneity has major influence on hydrological cycle, carbon balance and land-atmosphere interactions in the region (McDonnell et al., 2007;Govind et al., 2011;Stoy et al., 2013;Chapin et al., 2000;McGuire et al., 2002;Karlsen et al., 2016).Understanding spatial and temporal patterns in hydrological fluxes and state variables is becoming increasingly important in the context of intensifying use of boreal forests under the pressures from climate change (Bonan, 2008;Gauthier et al., 2015;Price et al., 2013;Spittlehouse, 2005;Laudon et al., 2016).Thus, model approaches that can effectively utilize available environmental data, open high-resolution GIS-data and remote-sensing products in hydrological predictions are necessary for climate-smart and environmentally sustainable use of boreal ecosystems (Mendoza et al., 2002).
Diverse modeling approaches are used to predict point scale, catchment and regional hydrological balance, which reflect the broad spectrum of practical needs and research questions addressed, as well as the historical development of hydrological models.The approaches range from lumped rainfall-runoff schemes (Beven and Kirkby, 1979;Bergström, 1992) to semi and fully distributed physically-based models (Vivoni et al., 2011;Panday and Huyakorn, 2004;Ivanov et al., 2004;Ma et al., 2016;Clark et al., 2015a, b;Best et al., 2011;Niu et al., 2011;Ala-aho et al., 2017).Lumped models are often based on conceptual representation of hydrological processes and calibrated against a few integrative measures, such as the stream discharge.They are computationally inexpensive but cannot well address the spatial variability of hydrological fluxes and state variables within a catchment.Distributed mechanistic models, on the other hand, use first principles to predict water flow and state variables through the landscape and can incorporate topography, soil texture and vegetation heterogeneity in their predictions (Samaniego et al., 2010).However, high computational costs and challenges to estimate spatially variable parameters hampers their use and performance (Freeze and Harlan, 1969;Montanari and Koutsoyiannis, 2012;Grayson et al., 1992;Reed et al., 2004).It has been questioned whether fully distributed models are suitable for applications requiring operative hydrological forecasts over large areas (Khakbaz et al., 2012) and semi-distributed models combining physical and conceptual elements are often suggested as practical solutions (Khakbaz et al., 2012;Gao et al., 2013;Savenije, 2010).
The common scientific questions in hydrological modelling are, as proposed by Clark et al. (2015a), related to describing and parametrising water and energy fluxes and representing landscape variability and hydrological connectivity at the spatial and temporal scales of the model discretization.The availability of good-quality high-resolution open data on land use, topography, vegetation, and soil characteristics has increased significantly in the recent decade.In Finland, for instance high-accuracy digital elevation models (DEM's) are openly available at 2m and 10m resolution (NSLF, 2017), reasonably good soil maps cover the country at scale of 1:20 000 or 1:200 000 (GSF, 2015), and the multi-source National Forest Inventory (mNFI, Mäkisara et al. (2016); Kangas et al. (2018)) provides information on various forest and site type attributes at 16m resolution throughout the country.To take full advantage of open GIS data and fine-resolution (i.e.tens to hundreds of meters) remote-sensing products of hydrological fluxes and state variables (Ryu et al., 2011;Herman et al., 2018), computationally efficient models capable of accounting for the landscape variability are necessary.Further, these models should be sufficiently generic in their parametrization and use standard meteorological data to allow their use on large, often data-sparce areas.As the appropriate level of detail is strongly driven by the research question or practical application at hand (Clark et al., 2011;Savenije, 2010), effective development of hydrological models requires moving from a specific model towards modular frameworks (Clark et al., 2015a;Wagener et al., 2001;Clark et al., 2011).
The increasing availability of high-resolution data on vegetation and its functioning paves way to improve the description of spatial and temporal variability of evapotranspiration (ET ), a major component of terrestrial water balance.Within a specific biome and climatic region, vegetation characteristics such as species composition and leaf-area index (LAI) have major influence on variability of ET (Williams et al., 2012;Zeng et al., 2018;Launiainen et al., 2016).In modern Land Surface Models, ET components are described either using a big-leaf framework or by describing the microclimatic gradients and exchange rates explicitly throughout a multi-layer canopy-soil system and upscaling these directly to ecosystem scale (Katul et al., 2012;Bonan et al., 2014).In both cases the upscaling of stomatal conductance g s and transpiration rate from leaf to canopy scale is based on physical arguments, and constrained by plant carbon economy (Cowan and Farquhar, 1977;Katul et al., 2012;Medlyn et al., 2012) and hydraulic architecture (Sperry, 2000;Tyree and Zimmermann, 2002).The non-linear dependency of ET components on vegetation characteristics and microclimate, however, remain mostly unresolved or are highly parametrised in most of the hydrological models where the bulk ET is commonly computed by using Penman-Monteith equation or as a crop or vegetation type dependent fraction of potential evaporation (Zhao et al., 2013;Fisher et al., 2005;Allen et al., 1998).
Thus, improving ET description by more physiologically-phased approach can be proposed as one potential area to reduce uncertainties in the predictions of hydrological budget and resulting streamflow and soil moisture patterns.Motivated both by scientific needs and potential practical applications, this study addresses two independent but inter-related objectives: First, we develop a generic model for daily ET in boreal forest and peatland ecosystems, and explore how daily and annual ET can be predicted based on plant functional traits, canopy LAI and open data on landscape structure and meteorological forcing.We distinguish between evaporative fluxes and transpiration, and predict canopy conductance G c and canopy transpiration rate by coupling the unified stomatal model (Medlyn et al., 2012) with simplified canopy radiative transfer theory (Saugier and Katerji, 1991;Kelliher et al., 1995;Leuning et al., 2008).We perform parameter sensitivity analysis and validate the model predictions against eddy-covariance (EC) measurements of stand-scale ET at ten boreal forest and peatland sites in Finland and Sweden (Launiainen et al., 2016).
Second, we extend the analysis to catchment scale and propose modular, semi-distributed Spatial Forest Hydrology (SpaFHy) model for predicting spatial and temporal patterns of hydrologic fluxes and state variables across the boreal catchments.The SpaFHy aims to provide reasonably simple, practically applicable and extensible framework that can effectively use open GIS data and basic meteorological data.We apply SpaFHy to 21 headwater catchments located throughout Finland to validate its predictions against daily stream discharge and annual ET derived from catchment water balance.The spatial variability of ET , snow water equivalent (SW E) and soil moisture, and temporal variability of stream discharge are demonstrated and potential applications finally discussed.Although developed for boreal ecosystems, the proposed methods can be extended to other biomes with appropriate parametrizations.

Model description
SpaFHy framework consists of three sub-models (Fig. 1).Hydrological processes in vegetation and two-layer topsoil are explicitly modelled for each grid cell, while Topmodel concept (Beven and Kirkby, 1979) is used to link grid cell and catchment water budgets, and to describe the baseflow and returnflow generation mechanisms.The SpaFHy submodels and equations are presented in the next sections, and complemented in Supplementary material.Throughout, we use notation where angle parenthesis y indicate spatial and y temporal averages of quantity y, and units mm correspond to kg H 2 O m −2 of surface area.

Canopy sub-model: above-ground water budget and fluxes at a grid cell
The hydrological processes in vegetation canopy, forest floor, snowpack, and in organic moss/humus layer and underlying root zone are solved for each grid cell using information on stand structure and soil type (Fig. 1; Canopy and Bucket -submodels).

P-M equation and ET
Total evapotranspiration is defined as sum of physiologically controlled transpiration (T r ) and physically regulated evaporation from wet canopy (E) and forest floor (E f ).To account for different controls of these processes, a three-source model is applied to describe ET at a grid-cell scale (Fig. 1).The Penman-Monteith (hereafter referred as P-M) equation gives each component of ET (mm d −1 ) as (Monteith and Unsworth, 2008) where L v is latent heat of vaporization (J kg −1 K −1 ), ρ w and ρ a are densities of liquid water and air (kg m −3 ), c p is the heat capacity of dry air at constant pressure (J kg −1 K −1 ), D is the vapor pressure deficit at air temperature (Pa) and A e is the available energy (W m −2 ), and ∆t daily timestep (86400 s).Depending on specific ET component E i , the surface conductance G i (m s −1 ) and aerodynamic conductance G a have different forms.For canopy layer, which contributes to T r and E, the G a represents efficiency of within-canopy turbulent transport and transport through laminar boundary layers on leaf surfaces, and is computed as a function of wind speed U , canopy height h c and LAI (Magnani et al., 1998;Leuning et al., 2008) (Suppl.S2).

Transpiration and canopy conductance
To calculate T r and resulting water uptake from the root zone, an estimate of the canopy conductance G c is needed.Analysing large corpus of leaf gas-exchange data through stomatal optimization arguments, Medlyn et al. (2012) proposed that leaf-scale stomatal conductance (g s , mol m −2 s −1 ) is related to leaf net photosynthetic rate (A, µmol m −2 s −1 ) as where C a is the atmospheric CO 2 mixing ratio (ppm), D (kPa) is vapor pressure deficit, g o residual (or cuticular) conductance and g 1 a species-specific parameter that depends on plant water use strategy.Noting that g o g s (Medlyn et al., 2012) and representing photosynthetic light response by saturating hyperbola (Saugier and Katerji, 1991), eq. ( 2) can be approximated as where A max (µ mol m −2 s −1 ) is the light-saturated photosynthesis rate, b (W m −2 ) the half-saturation value of photosynthetically active radiation (P AR), and molar density of air C air (mol m −3 ) converts units of g s to m s −1 .
Assuming PAR decays exponentially within the canopy, P AR(L) = P AR o exp(−k p L) (where L is the cumulative leaf area from canopy top, k p the attenuation coefficient and P AR o the incoming P AR above canopy) and neglecting vertical variations in D, eq. ( 3) can be integrated analytically over L yielding canopy conductance (m s −1 ) as where the first term of eq. 4 is the canopy-scale light and D response, the f (θ rew ), and f Y (-) are dimensionless scaling factors introduced to represent the effect of soil moisture availability (eq.6) and phenology (eq.8).Equation ( 4) shows that at a given LAI, G c is constrained leaf by leaf water use traits (via g 1 ), and photosynthetic capacity and light-response (via A max and b).Such traits are readily measurable by leaf gas-exchange and widely available in literature and in plant trait databases such as TRY (Kattge et al., 2011).Derivation of parameters is presented and predictions of eq. 4 compared against a common gas-exchange model in Suppl.S3.
Water use strategies and to lesser extent photosynthetic capacity of common coniferous and deciduous species in boreal forests differ (see e.g.Lin et al. (2015)).Thus, LAI-weighted effective values of g 1 and A max are calculated for a grid cell as where p is the parameter, the subscript c and d refer to conifer and deciduous trees, respectively, and the contribution of deciduous trees on total LAI.The seasonal cycle of LAI d is described using scheme based on accumulated degree-days (Launiainen et al., 2015) calibrated using leaf phenology observations in Southern and Northern Finland.
The effect of soil moisture availability on G c is based on sap-flow study on Scots pine and Norway Spruce in central Sweden (Lagergren and Lindroth, 2002) where θ rew (m 3 m −3 ) is relative plant available water θ rew , and x r (m 3 m −3 ) threshold below which reduction in G c occurs.
θ rew relates volumetric water content θ (-) in the root zone to soil-type depended field capacity (θ f c ) and wilting point (θ wp ) as The phenology factor f Y (-) describes seasonal acclimation of photosynthetic capacity as a function of delayed temperature sum (Kolari et al., 2007a) where t is time, Y (t) ( • C) describes the stage of development (Kolari et al., 2007b), T 0,y a threshold temperature and Y max the value at which full recovery of photosynthetic capacity occurs.The Y is accumulated from the beginning of the year and its rate of change is where τ is time constant of the recovery (Table 1).

Evaporation from forest floor
Forest floor evaporation E f is extracted from the organic moss/humus layer at the forest floor (Fig. 1).We compute E f (mm ) as where E f,0 is evaporation rate when moisture supply in the organic layer does not limit E f , and is calculated by eq. 1 where ), G a depends on surface roughness length and modeled U 0.5 m above the forest floor and G i now represents the conductance of saturated ground surface (G f ) and is calibrated against EC data from a boreal fen as described later.The factor f accounts for the decay of E f in drying organic matter as where θ org (m 2 m −2 ) is organic layer water content, and x r,o = 0.8 θ f c,org based on linear decrease of moss evaporation below threshold moisture content (Williams and Flanagan, 1996).

Interception, throughfall and evaporation from canopy storage
The canopy water storage is described as a single pool filled by interception I c of precipitation and snowfall (P ), and drained by evaporation/sublimation E and snow unloading U s (all in mm d −1 ).The change in canopy water storage W (mm) is The interception sub-model assumes that full storage is approached asymptotically (Aston, 1979;Hedstrom and Pomeroy, 1998) where c f (-) is canopy closure, W 0 the initial water storage, and canopy storage capacity W max = w max LAI (mm) linearly proportional to LAI.The empirical storage parameter w max (mm LAI −1 ) is known to be greater for rain and snow (Koivusalo and Kokkonen, 2002); if W exceeds W max of liquid water and daily mean temperature is above zero, the excess snow storage is removed as snow unloading and added into throughfall input to snow model.In snowfree conditions, all throughfall is routed to forest floor surface and provides input to Bucket sub-model.
Evaporation / sublimation from canopy storage is calculated by P-M equation (eq.1), where the G a is defined as for T r while the canopy surface conductance G i set infinite for evaporation from wet canopy, and computed for snow sublimation as in Essery et al. (2003) and Pomeroy et al. (1998) (Suppl.S4).

Snow accumulation and melt
Snowpack on the ground is modelled in terms of snow water equivalent (SW E, mm), a lumped storage receiving throughfall and unloading from the canopy and releasing water by snowmelt.The melt rate M (mm d −1 ) is based on temperature-index approach where T o = 0.0 • C is threshold temperature and T a daily mean air temperature.The melting coefficient K decreases with increasing canopy closure as Kuusisto ( 1984) where K m,o is the melting coefficient at open area.The snowpack can retain only a certain fraction of liquid water (Table 1 and Suppl.S4), and excess is routed to soil module.

Bucket model: topsoil water balance
The topsoil water balance at each grid cell is described as a two-layer bucket model (Bucket, Fig. 1).An organic layer of depth z org (mm), representing living mosses and poorly decomposed humus, overlies the root zone and acts as an interception storage for throughfall and snowmelt.Its volumetric water content θ org (m 3 m −3 ) is bounded between the field capacity θ f c,org and residual water content θ r,org and vary according to where I org is the interception rate, restricted either by throughfall or available storage space, and Q r,ex returnflow from the hillslope through the rootzone described below.All E f is extracted from the organic layer.
The water content θ in the root zone of depth z z (mm) changes according to where infiltration I f (mm d −1 ) and returnflow from the catchment sub-surface storage (sect.2.3) Q r provide input, and T r , and drainage D r outflows from the root zone.The maximum water storage is determined by root zone depth z s and porosity θ s , and I f restricted either by potential infiltration or available storage space.In case of infiltration or returnflow excess, the organic layer storage is first updated, and remaining routed to stream outlet without delay as surface runoff (Q s ).
Drainage D r (mm d −1 ) from root zone occurs whenever θ is above field capacity θ f c as (Campbell, 1985) where the saturated hydraulic conductivity K sat (mm d −1 ) and its decay parameter β depend on soil type.

Topmodel: integration from point to catchment level
To achieve computational efficiency and applicability at large-scale, lateral flow in the saturated zone is not explicitly solved but grid-cell and catchment water balances conceptually linked by Topmodel approach (Beven and Kirkby, 1979).In Topmodel, the catchment sub-surface storage is described as a single bucket (Fig. 1).The change in the average saturation deficit S (mm), i.e. the average amount of water per unit area required to bring the catchment sub-surface storage (below the root zone) to saturation, is where D r (mm d −1 ) is catchment average root zone drainage, Q b (mm d −1 ) the catchment baseflow and Q r (mm d −1 ) average return flow from the sub-surface storage.Assuming soil transmissivity is spatially uniform and decays exponentially with depth, the Q b becomes (Beven, 1997) The saturation deficit S (mm) of a grid cell is uniquely related to S by which implies that grid cells with high T W I have higher probability to become saturated, and the catchment saturated area fraction is related both to T W I distribution and to the amount of water in the catchment sub-surface storage.Furthermore, eq.
21 shows high value of T W I can result either from large contributing area or flat local topography.
At grid cells where saturation excess (S < 0) occurs, returnflow Q r = −S/∆t from the sub-surface storage is routed through the rootzone and organic layer and their water storages are sequentially updated at next ∆t.This creates an approximate feedback from local S, controlled by catchment water storage and topography, to topsoil water budget (sect.2.2) and delays drying of root zone and organic layer at lowland grid cells receiving Q r from the hillslope.
The specific discharge Q f (mm d −1 ) at catchment outlet is finally computed as where Q s is the catchment average surface runoff (sect.2.2).

Model inputs
SpaFHy requires daily mean air temperature T a ( • C), global radiation R g (Wm −2 ), relative humidity RH (%), wind speed (m s −1 ) and daily accumulated precipitation P (mm d −1 ) as forcing.The forcing can be either spatially uniform or vary for each grid cell in the spatial simulations.The available energy is computed from R g accounting for the effect of LAI on R n (Fig. 2a in Launiainen et al. ( 2016)), and PAR o = 0.5 × R g .
The model requires following variables to be provided at user-defined grid: 1. Canopy and Bucket -submodels -Conifer and deciduous tree 1-sided leaf-area index (LAI c and LAI d , respectively) -Organic layer depth, root zone depth and hydraulic properties (Table S1) 2. Topmodel -submodel topographic wetness index T W I masks of catchment area and permanent water bodies All the above variables are derived from open GIS-data available throughout Finland.The SpaFHy structure is modular and the three sub-models are linked via water fluxes, and feedbacks based on state variables such as θ rew (Fig. 1).Each sub-model can thus be used stand-alone when appropriate forcing data is provided.The model is written in pure Python 2.7/3.6 and uses element-wise operations of Numpy -arrays for all computations.The GIS-data for model initialization are given as raster arrays, while NetCDF -format is used for storing the model outputs that include daily grids of all state variables and fluxes.

Model parametrization and sensitivity analysis at stand scale
Parameters required by each sub-model are given in Table 1 with their generic values.We applied a sequential approach to determine the generic parameter set to describe above-ground hydrology of coniferous-dominated landscape.First, we derived likely ranges of Canopy sub-model parameters from the literature and predictions of a common leaf gas-exchange model (Suppl.S3).The rainfall interception capacity was calibrated against spatially averaged throughfall measurements (2001 -2010) made at the Hyytiälä research station in Juupajoki, Southern Finland (FIHy; Table .2, Fig. 2).An overview of the site is given by Hari and Kulmala (2005) while Ilvesniemi et al. (2010) describe the hydrological measurements.The parameter f in surface conductance for evaporation from wet soil surface (G f , eq. 10) was calibrated against eddy-covariance (EC) measured ET from a boreal fen site (FISii, Alekseychik et al. (2017)) located next to FIHy (Table 2).Monte-Carlo simulations (n=100), where parameter candidates were sampled from uniform distribution and objective function was set to minimize bias between modelled and measured values, were performed.
After parameter ranges were determined, a global sensitivity analysis was performed to identify the key parameters controlling annual ET and its components, and annual maximun SW E. For this analysis, the Canopy and Bucket modules were coupled and the resulting stand-scale model run with various parameter combinations using daily forcing data from FIHy (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010) as input.We used Morris method, a global extension of elementary effect test used to determine which model parameters are negligible, linear and additive, or non-linear or involved in interactions with other parameters (Morris, 1991;Campolongo et al., 2007).In Morris method, three sensitivity measures are calculated from the distribution of scaled elementary effects.The mean of distribution (µ) is the overall effect of a parameter on the output, and the standard deviation (σ) is effect of a parameter due to non-linearity or due to interactions with other parameters.Third measure is the mean of the distribution of the absolute values of the elementary effects (µ ) that provides ranking of parameters which is not biased by possible non-monotonic behavior of the model.The sensitivity measures are interpreted graphically together with rank parameters according to their overall influence; the intuitive interpretation is that the greater the absolute value of the measure the more important the parameter is for the studied model output.To ease graphical interpretation, standard error of the mean as where r is the number of trajectories, was estimated and used as suggested by Morris (1991).Analysis was conducted by using a Python package SALib (v.1.1.2;Herman and Usher, 2017).
The ranges of the 14 parameters considered in the sensitivity analysis are listed in Table 3.In the analysis, leaf area indices for conifers and deciduous were calculated from total 1-sided LAI and deciduous fraction.Each parameter was allowed to vary over eight levels, and 60 optimal trajectories were generated from 600 initial trajectories by the sampling scheme introduced by Ruano et al. (2012).In total, 900 samples were generated and the number of optimal trajectories was determined following Ruano et al. (2011).
After sensitivity analysis, most of the parameters could be fixed (those deemed less-influential), and only the 'generic' values for A max and g 1 in eq. ( 4) were confirmed by calibrating them against eddy-covariance (EC) -measured ET (years 2005ET (years -2007) ) at FIHy -site.The possible ranges of these parameters were constrained by physiological arguments.Monte-Carlo simulations (N=100), where parameters were sampled from uniform distribution and objective function was set to minimize bias between modelled and measured daily ET were performed.We considered only dry-canopy conditions, i.e. no rain during the current or previous day.

Model validation at stand scale
To validate how daily ET can be predicted across LAI, sitetype and latitudinal gradient using a single parameter set (Table 1), the stand-level model was run using daily meteorological data from nine additional EC-flux sites in Finland and Sweden (Table 2, Fig. 2).The sites range from dense mixed coniferous forests (SENor) to recently harvested stand (FICage4) and pristine fen peatland site (FISii), and the measurements, flux calculation and data post-processing have been described elsewhere (Launiainen et al., 2016;Minunno et al., 2016).For each site, LAI c , LAI d , h c and soil properties were set according to measured/inferred values, and predicted daily growing-season (May-Oct) ET in dry-canopy conditions (ET T r + E f ) was compared to measured.At FIHy, the soil moisture in the root zone was measured continuously, and SW E recorded bi-weekly during five winters and used to compare respective model predictions.

Model validation at catchment scale
To address how well SpaFHy can predict daily specific discharge and annual partitioning of P into ET and Q f at catchment scale, we applied the model to 21 small boreal headwater catchments located throughout Finland (Fig. 2, Table S2) using same generic parameter set in the stand-level validation (Table 1).All the catchments belong to the Finnish network for monitoring water quality impacts of forestry (Finér et al., 2017), and their characteristics can be found in Supplementary material (Table S2).The water levels at v-notch weirs were measured continuously at the catchment outlets by limnigraphs or pressure-sensors, and manual reference measurements ware taken ca.20 times per year adjacent to water quality sampling and used to calibrate the weir water level data whenever necessary.Weir equations and catchment area were used to convert water level to specific discharge Q f .In absence of in-situ weather data, daily 10 x 10 km grid data provided by Finnish Meteorological Institute were used as model forcing taking the values from a gridpoint nearest to the catchment outlets.Since wind speed was not available, it was set to constant value of 2 ms −1 resembling annual mean 2 m wind speed in Finland.The DEM pre-processing, defining of the catchment boundaries and the calculation of T W I based on the aggregated DEM were conducted in WhiteBox GIS programme (Lindsay, 2014).Pre-processing included consideration of the road and stream intersections derived from the Topographical Database (NSLF, 2017), which were burned into the DEM to account for culverts and ensure continuous stream network.Further, all water elements were burned into the DEM with 1 meter upper threshold and a decay factor accounting for possible miss-aligned stream data.The filling of artificial pits in DEM was conducted using 'Fast Breach Depressions' tool (Lindsay, 2016) and the flow direction and flow accumulation (a) rasters were calculated with the D-infinity method (Tarboton, 1997).The T W I was finally calculated by eq.21 and small lakes within the catchments, derived from the Topographic Database (NSLF, 2017), were reset as nodata and omitted from further computations.The needle and leaf mass rasters were converted into LAI c and maximum deciduous tree LAI LAI d,max using specific 1-sided leaf-areas for pine, spruce and birch (6.8, 4.7 and 12.0 m 2 kg −1 , respectively; Härkönen et al. (2015)).The canopy closure and h c were obtained directly from the mNFI data.
Topsoil classification was derived from soil maps and peatland boundaries.Soil information is provided for parts of Finland in 1:20 000 scale while the whole Finland is covered with a coarser 1:200 000 scale soil map (GSF, 2015).Peatland classification is available as detailed polygon elements from the Topographical Database.The soil information were transformed to the 16 m grid based on the majority principle, and then re-classified into four classes: coarse, medium and fine-textured mineral and organic peat soils whose hydrologic properties are given in Table S1.Fine-textured soils correspond to clayey and silt soils, whereas coarse-textured are fine sand and coarser.Majority of the mineral soils in the study catchments belong to the medium-textured class (Table S2).

Calibration of Topmodel against specific discharge
Catchment-specific calibration was performed to determine the effective soil depth m of Topmodel, a parameter that defines the shape of Q f recession and catchment average storage deficit S (eq.22).The parameter T o was fixed to 0.001 ms −1 since it was found not markedly affect the model performance, as also observed elsewhere (Beven, 1997).The m was calibrated against measured daily specific discharge using Monte-Carlo sampling from uniform distribution (N=100).We used modified Willmott's index of agreement (Krause et al., 2005) as an objective function to quantify the goodness of fit where d is in a range of 0 to 1, the higher the value, the better the match is; Q f,i and Q m,i are modelled and measured specific discharges at day i, and Q m represents temporal average of the measurements.This model goodness statistics provided visually determined better fits of streamflow recession than other commonly used statistical criteria, e.g.Nash-Suchliffe model efficiency that was overly sensitive to high-flow peaks and affected by potential biases in P .The initial state of the model was set through one year spin-up period.The value of m significantly affected the dynamics of specific discharge Q f (t) and S(t) but had negligible impact on catchment ET or Q m at annual scale.

Sensitivity analysis at stand scale
The sensitivity measures µ and σ for maximum SWE and annual ET and its components are shown in Table 4, and the ranking of parameters (via µ ) in Supplementary Fig. S1.
Total LAI was ranked the most influential parameter for all studied Canopy sub-model outputs.In addition to LAI, the parameters that affect leaf level water use (g 1 , z s , A max , and b) were among the most influential parameters for total ET and transpiration.The most influential parameters for ground evaporation E f were LAI and k p , which jointly define radiation availability at the ground layer.LAI also affects wind speed and thus aerodynamic conductance at the ground layer.In addition, surface conductance for wet forest floor G f and z s,org and θ f c,org that define water storage capacity of the organic layer were significant for E f .The most influential parameters for interception evaporation E were LAI, w max , f d , h c , and w max,snow that define interception capacity and subsequent evaporation/sublimation of rain and snow.The most influential parameters affecting annual maximum snow water equivalent were LAI, w max,snow , f d , w max , and h c .
LAI had also the largest σ meaning either interactions with other parameters or strong non-linearity.In case of ET and T r , coefficient of variation (σ / µ -ratio) was over 1.0 and for E, E f , and SW E it was smaller but over 0.5.The most influential parameters of all studied outputs had the coefficient of variation over 0.5.Non-monotonic behavior (i.e.µ / µ -ratio is significantly different from unity) of the model was only observed in case of ET for LAI.

Validation at stand-scale
The predicted daily dry-canopy ET and root zone moisture content are compared against 10 years of measurements at the pine-dominated FIHy -site in  (Pomeroy et al., 1998;Essery et al., 2003).
The ET predictions for the nine additional EC-sites are shown in Fig. 5.The growing season (doy 120 -273) dry-canopy ET is reasonably well predicted compared to independent observations across broad LAI -range (from 0.7 to 6.8 m 2 m −2 ) and over latitudinal and site-type gradient (Table 2, Fig. 2).At the youngest, recently clearcut site FICage4 the model underestimates ET, while slight overestimation is observed in particular at the northernmost, old-growth Scots pine site on coarse textured soil (FISod).In terms of explained variability, the model performance is the weakest at SESky2 (spruce), FIKal and FILet (drained peatland forests), potentially because ill-represented moisture limitations of transpiration and/or that of E f .The nonlinear behavior at SENor and less clearly at SESky2 and FILet is primarily caused by slower than observed spring recovery at these sites having high abundance of Norway spruce.As the Norway spruce has observed to recover more rapidly from winter dormancy than pine (Linkosalo et al., 2014;Minunno et al., 2016), this can be partly related to biased phenology-model that is based on Scots pine (Kolari et al., 2007a).
Also ET at the pristine fen peatland site FISii, where E f T r , was accurately predicted when moisture limitation of E f was neglected (f = 1 in eq.10).Such case can be expected due to strong capillary connection between peat moss (Sphagnum sp.) and shallow water table maintained by lateral inflows from the surrounding landscape and weak drainage (Rouse, 1998;Ferone and Devito, 2004).When the organic layer moisture content feedback to E f was activated, ET at FISii was frequently underestimated during summer dry spells (not shown); in the point-scale simulations this represent the case where the organic layer water storage is recharged only by P .
Overall, the model performance at stand scale was satisfactory and dry-canopy ET well predicted over range of forest sites and climatic gradient in Finland.This suggests that the proposed three-source ET formulation and its generic parametrization for T r , E f and snow interception should be scalable over landscape scale variability of LAI, site types and latitude-driven weather forcing.Since EC measurements are know to be problematic during rainfall events (van Dijk et al., 2015;Kang et al., 2018), the comparison of stand-level ET was restricted to dry-canopy conditions.

Catchment water balance and specific discharge
On annual scale, changes in catchment water storage are negligible compared to annual ET and Q f , and water balance approach provides an independent check for the upscaled ET predictions at catchment level.Fig. 6 shows the comparison of modelled and water-balance based annual evapotranspiration fraction ET /P for the 21 headwater catchments across Finland (Fig. 2, Table S2).Results show a close agreement between measured and modelled P partitioning across the catchment space, especially considering the uncertainties in both axis.The uncertainty range of modelled ET /P implies the impact of model parameter uncertainty.The uncertainty range in Fig. 6 was derived by varying the most influential parameters for total ET and its partitioning (LAI, g 1 , w max , w max,snow ) by ± 20 % and grouping the combinations into 'high' and 'low' ET scenarios, respectively.While the model is mass-conserving, uncertainty of ET /P derived from catchment water balance is linearly proportional to uncertainty of Q f derived from streamflow measurements and catchment area.Also systematic and random errors in the annual P cause respective uncertainties in ET /P .In Fig. 6 the horizontal errorbars correspond to modest 10 % uncertainty assumed for P and catchment area.
Overall, the model predictions are reasonably good across the catchment space.Stepwise linear regression was tested to explain the annual residuals by catchment characteristics in Table S2 but no significant relationships were found.Also interannual variability of ET /P was well captured for majority of the catchments (not shown).
Figure 7 shows specific discharge and modelled soil moisture at catchment C3 Porkkavaara in Eastern Finland (Table S2), over two years characterized by wet (2012, P 452 mm in June-Sept) and dry (2013, P 246 mm) growing seasons.In 2012, the high snow accumulation resulted into stronger streamflow peak, and frequent rainfall events kept the catchment average root zone moisture θ ≥ 0.3 m 2 m −2 throughout the year.Also Q f remained significantly higher throughout the summer compared to 2013, and responded rapidly to rainfall.During the drier 2013, transpiration depleted the root zone moisture well below field capacity and θ dropped frequently to ∼0.15 m 2 m −2 in June -August.The model was well able to predict the spring Q f peaks and recession curve, and also rainfall-induced peaks during the wet summer.During drier conditions, however, the small-magnitude peaks in summer Q f were not well captured by the model.This suggests too high Bucket storage capacity and thus underestimated fraction of saturated area that contributes to overland flow during and after precipitation events.This is, however, not a general behaviour of the model as better comparison between measured and modelled specific discharge was observed at several other catchments (Fig. S2).

Soil moisture
To illustrate how vegetation, soil and topography create within-catchment variability to local water fluxes and state variables, the relationship between θ and its spatial standard deviation σ θ at C3 Porkkavaara is shown in Fig. 8 for the two hydrologically contrasting years.Snapshots of spatial variability of θ and local saturation deficit of Topmodel (eq.22) are further shown for dry and wet conditions in Fig. 9.
During winter, root zone moisture content decreases and its spatial variability is dampened by slow drainage.The onset of snowmelt is followed by infiltration peak and saturated soils nearly throughout the catchment (Fig. 7 -9).This leads to rapid increase of σ θ , mainly because of spatial variation in soil porosity.In wet 2012, drainage rapidly decreased σ θ after snow melt, while the spatial variability was preserved until early July in the drier 2013.The latter result is because of spatially heterogeneous transpiration rate that created spatial variance of soil moisture and compensated for the dampening effect of drainage until ca.doy 180.After that σ θ started to decrease because transpiration at grid cells characterized by coarse and medium-textured soil and high LAI (Fig. 3) become soil-moisture limited (eq.7).In 2013 summer when θ was most of the time well below field capacity, the rainfall events tend to dampen spatial variability of soil moisture (Fig. 7).In wetter conditions (most of 2012, autumn 2013), however, the effect of infiltration is opposite and resembles that of spring snowmelt.
As a result, there is clear hysteresis of σ θ with respect to antecedent θ in the dry year while such patterns are less visible in moist conditions.This indicates soil and vegetation variability can in the model both create or destroy spatial variability of soil moisture, as has been proposed both by theoretical arguments (Albertson and Montaldo, 2003) and analysis of soil moisture observations (Teuling and Troch, 2005).To conclude, the role of landscape heterogeneity as driver of modelled soil moisture variability depends on antecedent soil moisture conditions and season; during drier spells variability is primarily driven by heterogeneous vegetation and plant water use, while soil type and topography become the primary controls in wet conditions and outside growing season (Seyfried and Wilcox, 1995;Teuling and Troch, 2005;Robinson et al., 2008).

ET and snow
The predicted spatial variability of evaporation fraction ET /P and its components are illustrated in Interception of rainfall and snow contributes from less than 5 to 30 % of long-term P , which is in line with measurements from boreal forests (Barbier et al., 2009;Toba and Ohta, 2005).The linear scaling of interception capacity with LAI and asymptotic approach of full storage (eq.13), and temporal distribution of precipitation lead to the near-linear increase of E/P with increasing LAI (Fig. 10).At grid cells with high fraction of deciduous trees, low wintertime LAI leads to weaker snow interception and smaller annual ET /P compared to coniferous-dominated stands.
The model predictions suggest transpiration contributes from < 10 to more than 35 % of annual P being the largest ET component in stands whose LAI > 1.5 m 2 m −2 (Fig. 10).The shape of T r to LAI -response is to most extent caused by saturation of G c because of light limitations in dense stands (eq.4).The more liberal water use strategy of deciduous species , Table 1) is reflected as higher transpiration rate at grid cells where deciduous trees form a significant part of total LAI.Moreover, the lower envelope of the points occur at grid cells corresponding to coarse-textured soils (Fig. 3), where drought limitations become most frequent.This is visible also in ET /P at LAI > 2 m 2 m −2 .
Evaporation from forest floor ET /P decreases asymptotically with LAI, showing complementary relationship to T r , as expected by decreased available energy in denser stands.The upper envelope curve corresponds to grid cells with high T W I and large tendency to be permanently saturated due returnflow from the hillslope.In these grid cells E f is mainly determined by available energy; however, rapid drying of the forest floor in sparse stands between rainfall events decreases E f /P and explain it less steep decrease with LAI at grid cells receiving less frequent or no returnflow (lower T W I).
The spatial pattern of maximum SW E (Fig. 11) indicate snow accumulation in the densest stands (LAI > 7 m 2 m −2 ) was ∼ 75 % of that on open areas; the exact fraction was found sensitive to winter weather conditions being lowest during mild winters in the southern catchments, and in years with smaller annual snowfall.The predicted impact of forest canopy on snow accumulation is in good agreement with observational studies from similar climatic conditions in Finland and Sweden (Koivusalo and Kokkonen, 2002;Lundberg and Koivusalo, 2003), although also higher snow interception losses have been reported (see Kozii et al. (2017) for summary).The near linear increase of snow interception and resulting decrease of SW E with LAI is supported by Hedstrom and Pomeroy (1998); Pomeroy et al. (2002).

Modeling ET at stand and catchment scales
At stand-scale, SpaFHy was shown to well reproduce daily ET measured by eddy-covariance technique at several forest and peatland sites in Finland and Sweden (Fig. 4 and 5).The good performance using generic parametrization derived mainly from literature sources and process-specific data suggests the model is capable of accounting for the key drivers of temporal and site-to-site variability of ET .The sensitivity analysis reveals that for given meteorological forcing, total LAI is the primary parameter affecting ET and its partitioning into component fluxes.In case of transpiration, the dominant role of LAI and parameters defining leaf water use efficiency (g 1 and A max ), and insensitivity to parameters related to aerodynamic conductance of the P-M equation (eq. 1) indicate variations in T r are mainly governed by that of canopy conductance (eq.4).The root zone depth, soil hydraulic properties and size of interception storage in the organic layer (z org and θ f c,org ) are important for probability of drought occurrence and consequent reduction of transpiration (Table 4, Fig. 10).As rooting depths vary across species, site types and ecosystems (Gao et al., 2014) and soil heterogeneity is not fully represented by existing soil maps, uncertainties of these properties are large in general.It was, however, recently shown the root zone storage capacity can be estimated from satellite-based evaporation (Wang-Erlandsson et al., 2016), which may in future provide data to constrain these model parameters.
The multiplicative formulation for canopy conductance (eq.4) was developed by coupling the commonly used unified stomatal model (Medlyn et al., 2012) and leaf-scale light response with simplified canopy radiative transfer scheme (see Suppl.

S3
).The approach accounts the non-linear scaling between G c and g s similarly as Saugier and Katerji (1991) and Kelliher et al. (1995).To derive bulk surface conductance for remote-sensing applications, Leuning et al. (2008) combined their G c scheme with a ground evaporation model based on equilibrium evaporation (Priestley and Taylor, 1972).They showed that after site-specific optimization, the dry-canopy ET was accurately predicted by P-M equation across different vegetation types.
That particular model, however, still requires an arbitrary and non-measurable maximum g s and few other parameters to be specified or calibrated.In our work g s and its response to D were derived from stomatal optimization arguments and are tightly constrained by plant water use traits and photosynthetic capacity.These traits start to be widely available in databases such as TRY (Kattge et al., 2011), and can also be readily measured using leaf gas-exchange techniques.Due these constraints to g s , we consider eq. ( 4) as a major advancement of the Leuning et al. (2008) approach.The good comparison between modelled and measured dry-canopy ET for sites having strongly different T r /ET and E f /ET -ratios (Fig. 4 and 5) are indeed supportive for the proposed G c formulation.However the comparison was done within a single vegetation type and further evaluation across ecosystem types are necessary to extend the approach outside boreal forests.
The sensitivity analysis (Table 4 and Fig. S1) proposes the P-M equation could be replaced with simpler approaches.Making the assumption that canopy is well-coupled to the atmosphere, reasonable for aerodynamically rough boreal forests, leads , where p a (kPa) is the ambient pressure.Also evaporation from the ground and canopy storage were found relatively insensitive to aerodynamic terms, which suggests they could be computed proportional to equilibrium evaporation ∆+γ , where α i (-) is a proportionality factor calibrated against measurements, and R n,i available energy at the ground level.Moving to such approaches would relax input data requirements by eliminating the canopy height and wind speed from model forcing.
Open GIS data on LAI, species composition, soil type and topography was used to apply SpaFHy at 16 x 16 m grid size to 21 headwater catchments in Finland.Results indicate the model well reproduces the variability of annual evaporation fraction across catchments (Fig. 6), as well as inter-annual variability at most of the studied catchments (not shown).It should be noted that the variability of annual ET /P across the catchment space is dominated by latitudinal climate gradient and further testing across different catchments on similar climatic conditions is needed.
Validation of spatial predictions of θ, ET or SW E (Fig. 9 -11) was not attempted in this work.This would require either extensive spatially distributed and continuous in situ measurements, or high-resolution (i.e.order of tens of meters) remote sensing data that can be already obtained by near-ground microwave radiometry or low-frequency radars using unmanned aerial vehicles as a platform (Robinson et al., 2008).Also ongoing advances in satellite-based soil moisture (Chen et al., 2014) and ET products (Hu et al., 2015) could be used to evaluate the modelled spatial patterns and temporal evolution of these hydrological components.
The results of site and catchment scale validation suggest that ET and water budget partitioning in boreal forest-dominated landscape can be reasonably well predicted by the model based on generic parametrization, which is advantageous for scalability and applicability of the model for areas and locations where data is scarce or lack for model calibration.Moreover, the model-data comparison at catchment scale propose ET components and water budgets can be upscaled from stand to catchment scale using relatively simple mechanistic approach that derives characteristics of the modelling domain from open GIS data.

Capabilities and limitations of the model framework
This study presented a semi-distributed model for boreal forest hydrology at stand and catchment scales (Fig. 1).The model consists of three independent components: a Canopy model for above-ground hydrology, a Bucket model for topsoil water balance and Topmodel for point to catchment integration.The modularity of SpaFHy provides clear advantages since all model components are independently parametrized which allows their stand-alone development and use, as well as inclusion to other distributed or lumped hydrological models.Moreover, parameters of each sub-model were obtained separately and calibrated based on good-quality data that clearly enhances the predictive power of SpaFHy by reasonably constraining the degree of freedom in model parametrization (Jakeman et al., 2006;Jackson-Blake et al., 2015).
In SpaFHy, the above-ground hydrology and root zone water balance (eq.16 & 17) are solved distributively (Fig. 1), which propagates the spatial variability of vegetation (LAI, c f , species composition) and soil type into the local hydrological fluxes, SW E and organic layer and root zone water contents.Applied stand-alone, such approach would assume grid cell water balances are independent from each other, and omits the lateral flows and the topographic position of a grid cell on a hillslope.The role of Topmodel (sect.2.3) can be considered as a non-linear streamflow generation routine, which delays average root zone drainage signal D r leading to realistic response of streamflow to P as controlled by TWI distribution.The other catchment properties are lumped into the parameter m, the effective subsurface water-conducting depth.It is this parameter that primarily controls both the shape of rainfall-runoff response and streamflow recession.The SpaFHy can thus be used as a simple catchment model to predict the signals of vegetation changes, forest management or varying climatic drivers on streamflow at daily or longer time scales.Indeed, the daily time series of streamflow (Fig. 7 and Fig. S2) were well reproduced for majority of the 21 studied catchments although m was the only parameter specifically calibrated for each catchment (Suppl.Table 2).
On the other hand, SpaFHy can assist in mapping how soil saturation may vary spatially and temporally as response to weather forcing (Fig. 9).The T W I-based scaling in Topmodel is used to predict magnitude and location of returnflow formation based on state of the catchment sub-surface storage.The spatial Q r field is then used to update Bucket sub-model water storages and θ at respective grid-cells.In this way, SpaFHy can be used to predict local soil saturation that depends on both local (via vegetation and soil characteristics) and approximative landscape (via topography) controls (Fig. 9).In essence, the effect of returnflow formation is to delay drying of gridcells that receive water from the surrounding landscape.Depending on T W I-distribution and value of m, this conceptualization implies that some gridcells never receive water from the surrounding landscape(those with low T W I) while some receive Q r in highflow conditions but not in baseflow conditions.At the other extreme, there are permanently inundated areas (high T W I) that contribute constantly to overland flow.
We emphasize that linking grid-cell water balances through Topmodel is conceptual rather than physically correct approach (Beven, 1997;Seibert et al., 1997;Kirkby, 1997), and driven by the goal to develop a simple and practically applicable representation of topographic controls of soil moisture.Future work should explore whether m can be related to catchment characteristics to derive a more generic parametrization for Topmodel, as well as analyse the impact of parameter uncertainty on streamflow and saturated area predictions.For applications requiring more rigorous treatment of sub-surface flows, the Topmodel can be replaced with 2D ground water flow schemes.Fig. 9 and 10 show that landscape position (accounted via T W I) can markedly affect grid cell soil moisture and ET .In this work, other topographic controls were omitted for simplicity.While likely to have small impact for annual catchment water balance, including topographic effects on radiation (Dubayah and Rich, 1995) is presumed to alter the spatial patterns of ET and θ.In addition, the shading by vegetation at the neighbouring grid cells should be considered to derive a more comprehensive understanding on hydrological variability on the landscape.Also adding sub-models to simulate spatial and temporal patterns of soil temperature and frost depth, vegetation productivity and carbon balance would be relatively straightforward future developments.
As shown in this work, the mNFI data (Mäkisara et al., 2016) can provide estimates of LAI, canopy height, site type and conifer/deciduous composition at 16 x 16m resolution throughout Finland.Härkönen et al. (2015) compared mNFI-based LAI estimates against ground-based estimates and MODIS AI, and found good agreement between the methods.Consequently, the mNFI data can provide an easy way to obtain vegetation characteristics for hydrological and biogeochemical models at spatial scale currently unresolved by e.g.MODIS and other satellite products.Similar high-resolution data on forest resources is openly available also from other Nordic countries (Kangas et al., 2018).

Potential applications
The proposed modular framework can provide support to a variety of questions benefiting from spatial and temporal hydrological predictions.These include, but are not limited to: (1) predicting soil moisture necessary for forecasting forest soil trafficability (Vega-Nieva et al., 2009;Jones and Arp, 2017), precision forestry, and confronting climate-induced risks (Muukkonen et al., 2015); (2) identifying how saturated areas, considered as biogeochemical and biodiversity hotspots particularly sensitive to negative environmental impacts of human activities, evolve in time (Laudon et al., 2016;Ågren et al., 2015); (3) addressing the impacts of forest structure, management and climate change on ET partitioning, streamflow dynamics and soil moisture (Zhang et al., 2017;Karlsen et al., 2016); (4) supporting water-quality modelling in headwater catchments (Guan et al., 2018); and (5) providing starting point for developing spatially distributed forest productivity and sustainability framework that combines open data streams, statistical approaches and mechanistic models.Moreover, we propose the Canopy sub-model, in particular the leaf-to-canopy upscaling of canopy conductance, to be tested more widely in other ecosystems.

Conclusions
A distributed hydrological model framework for predicting ET and other hydrological processes from a grid cell to a catchment level using open GIS-data and daily meteorological data was presented and validated for boreal coniferous-dominated forests and peatlands.SpaFHy consists of coupled, stand-alone modules for aboveground, topsoil and subsurface domains.
An improved approach to upscale stomatal conductance to canopy scale was proposed, and a generic parametrization of vegetation and snow-related hydrological processes for Nordic boreal forest and peatland ecosystems derived.With the generic parametrization, daily ET was well reproduced across conifer-dominated forest stands whose LAI ranged from 0.2 to 6.8 m 2 m −2 .Predictions of annual ET were successful for the considered 21 boreal headwater catchments in Finland located from where m (mm) is a scaling parameter reflecting the effective water-conducting soil depth, T o the soil transmissivity at saturation, and Q o (mm d −1 ) the baseflow rate when S is zero.The T W I represents the catchment average of local topographic wetness index T W I defined by the natural logarithm of the area draining through a grid cell a from upslope and tangent of the local surface slope α(Beven and Kirkby, 1979)

2. 7 . 1
Processing of GIS -data Example of GIS-data used to set up the model for catchment C3 Porkkavaara in Eastern Finland are shown in Fig. 3.The catchment boundaries and T W I were derived from DEM provided by National Survey of Finland (NSLF, 2017).The DEM original resolution was 2 m or 10 m depending on catchment location.The resolution was aggregated with the mean elevation value into 16 m resolution which corresponds to the resolution of the multi-source National Forest Inventory of Finland (mNFI) data.The mNFI data provides the essential data layers for the model, e.g.stand volume, basal area, mean height, age, site fertility class and estimates of root, stem, branch and needle/leaf biomasses for pine, spruce and aggregated deciduous trees at 16 m resolution throughout Finland Mäkisara et al. (2016) and thus 16 m was chosen as the baseline resolution.

Fig. 4 .
The results indicate the model reproduces well the observed seasonal patterns of ET and θ both during calibration (2005 -2007) and validation period.The regression plots indicate ET predictions have negligible bias and well represent the variability, while the soil moisture changes are not fully captured.The SW E was also well reproduced by snow model parametrized by literature values

Fig. 10 .
The model results, averaged over 2006 -2016 period, reveal strong sensitivity of component ET fluxes to stand leaf area index, and secondary impacts of soil type and topography.The model predicts ET /P increases non-linearly with LAI and varies from >0.25 at grid cells where LAI < 1 m 2 m −2 to ∼0.65 at locations where the standing tree volume and LAI (Fig. 3) are largest.The shape of LAI-response results from the non-linear scaling of component fluxes with LAI, which also explain the inflection point at LAI ∼3 m 2 m −2 .

Figure 1 .
Figure 1.Structure of SpaFHy.At each grid cell, above ground and topsoil hydrology is solved by Canopy and Bucket sub-models whereas lumped Topmodel is used to model saturated zone.The arrows correspond to interfacial fluxes: P precipitation; T f throughfall to virtual snowpack; Ip potential infiltration to organic layer; I f infiltration to root zone; Dr drainage to saturated zone; E evaporation/sublimation of canopy storage; E f evaporation from ground; Tr transpiration; Qr returnflow; Qs surface runoff; Q b baseflow.

Figure 2 .
Figure 2. Location of the forest and peatland eddy-covariance sites and the 21 boreal headwater catchments used in the study.

Figure 3 .
Figure 3. Spatial data at 16 m resolution used to set up the model for the catchment C3 Porkkavaara in Eastern Finland (see TableS2).LAI Figure 3. Spatial data at 16 m resolution used to set up the model for the catchment C3 Porkkavaara in Eastern Finland (see TableS2).LAI is total 1-sided leaf-area index; f d deciduous fraction; hc canopy height; Elev elevation; T W I topographic wetness index; soiltype refers to TableS2.Rasters overlay topographic basemap provided by National Survey of Finland and the scale of x and y axis is meters.

Figure 4 .
Figure 4. Modeled vs. measured dry-canopy ET at FIHy (top), root zone water content θ (middle) and snow water equivalent SW E (bottom).As soil freezing is not modelled, comparison of θ is restricted to conditions when measured soil temperature was ≥ 0 • C.

Figure 5 .
Figure 5. Scatter plots between modeled and observed daily stand-level ET during growing-season at the eddy-covariance flux sites in Finland and Sweden (Table 2).The title of each panel shows LAI (maximum deciduous LAI in parenthesis).The slope s and R 2 of linear regression forced through the origin and mean error M E are given and dashed line is the 1:1 line.Only dry canopy conditions, i.e. no rain during the day or previous day are included.At pristine fen peatland site FISii, E f was assumed non-limited by organic layer moisture.Color coding is according to transpiration to ET -ratio Tr/(Tr + E f ).

Figure 6 .
Figure 6.Modeled annual catchment evaporation fraction ET mod /P compared to that inferred from catchment water balance ET wb /P .The vertical and horizontal errorbars show effect of parameter uncertainty and that of catchment area and P , respectively(see text).The colors refer to latitude and symbol size to catchment mean LAI (from 0.2 to 4.6 m 2 m −2 ).Using median year for each catchment (N=21), the respective statistics are: slope s 0.99±0.30,R 2 0.67, RM SE 0.06, M E -0.003.

Figure 7 .
Figure 7. a) Measured (black) and modeled (red) specific discharge Q f , daily precipitation P (black bars) and mean snow water equivalent SW E (blue); b) mean volumetric soil moisture θ and its spatial standard deviation σ θ (blue) over two hydrologically contrasting years at C3 Porkkavaara, Eastern Finland.In b) the the grey range shows the inter-quartile range.The points correspond to dates in Fig. 9.The Willmot's index of agreement (eq.23) for specific discharge over 2012 -2013 period is 0.77.

Figure 8 .
Figure 8.The relationship between catchment daily mean soil moisture θ and its spatial standard deviation σ θ over the course a wet (2012) and dry (2013) year at C3 Porkkavaara.The color scale shows day of year (doy).

Figure 9 .
Figure 9. Snapshots of soil moisture patterns during wet and dry conditions (Fig. 7) at C3 Porkkavaara.Water content θ in the root zone (top) and local saturation deficit S of Topmodel (bottom).

Figure 10 .
Figure 10.Spatial variability of evaporation fraction ET /P and its components at C3 catchment in Eastern Finland from a long-term (2006-2016) run (left).The relationship of component fluxes interception evaporation E, transpiration Tr and forest floor evaporation E f on LAI (right) is modified with spatial variability of soil type, proportion of deciduous trees (LAI d /LAI), and topographical wetness index (T W I).

Figure 11 .
Figure 11.Spatial variability of maximum snow water equivalent (SW E) at C3 (left).The SW E relative to open area scales near-linearly with wintertime leaf-area index LAI (right).
sided leaf-area index (deciduous LAI in parenthesis); P is mean annual precipitation; Ta mean annual air temperature and soil type refers to Table2in the main document.* for runs at FIHy, site-specific field capacity θ f c = 0.3 and wilting point θwp = 0.1 corresponding to main root zone(Launiainen et al., 2015) were used.

Table 4 .
Sensitivity of Canopy sub-model predictions to parameter variability.Mean (µ) and standard deviation (σ) of the distribution of elementary effects for evapotranspiration (ET ), transpiration (Tr), evaporation from canopy interception (E), ground evaporation (E f ), and annual maximum snow water equivalent (SW E).Negative sign of µ indicate output variable decreases when parameter value increases.Units are in mm a −1 except for SW E (mm).

Table 1 .
Generic parameter set used in stand and catchment scale simulations.