Towards improved parameterization of a macroscale hydrologic model in a discontinuous permafrost boreal forest ecosystem

Modeling hydrological processes in the Alaskan sub-arctic is challenging because of the extreme spatial heterogeneity in soil properties and vegetation communities. Nevertheless, modeling and predicting hydrological processes is critical in this region due to its vulnerability to the effects of climate change. Coarse-spatial-resolution datasets used in land surface modeling pose a new challenge in simulating the spatially distributed and basin-integrated processes since these datasets do not adequately represent the smallscale hydrological, thermal, and ecological heterogeneity. The goal of this study is to improve the prediction capacity of mesoscale to large-scale hydrological models by introducing a small-scale parameterization scheme, which better represents the spatial heterogeneity of soil properties and vegetation cover in the Alaskan sub-arctic. The small-scale parameterization schemes are derived from observations and a sub-grid parameterization method in the two contrasting sub-basins of the Caribou Poker Creek Research Watershed (CPCRW) in Interior Alaska: one nearly permafrost-free (LowP) sub-basin and one permafrost-dominated (HighP) sub-basin. The sub-grid parameterization method used in the small-scale parameterization scheme is derived from the watershed topography. We found that observed soil thermal and hydraulic properties – including the distribution of permafrost and vegetation cover heterogeneity – are better represented in the sub-grid parameterization method than the coarse-resolution datasets. Parameters derived from the coarse-resolution datasets and from the sub-grid parameterization method are implemented into the variable infiltration capacity (VIC) mesoscale hydrological model to simulate runoff, evapotranspiration (ET), and soil moisture in the two sub-basins of the CPCRW. Simulated hydrographs based on the small-scale parameterization capture most of the peak and low flows, with similar accuracy in both sub-basins, compared to simulated hydrographs based on the coarseresolution datasets. On average, the small-scale parameterization scheme improves the total runoff simulation by up to 50 % in the LowP sub-basin and by up to 10 % in the HighP sub-basin from the large-scale parameterization. This study shows that the proposed sub-grid parameterization method can be used to improve the performance of mesoscale hydrological models in the Alaskan sub-arctic watersheds.


Introduction
The sub-arctic region of Interior Alaska lies in the transition zone between the warm temperate region to the south and the cold arctic region to the north.This region is underlain by discontinuous permafrost and is very sensitive to climate warming (Hinzman et al., 2006).Along with climate and geology, topography is a significant factor that controls the distribution of permafrost in Interior Alaska, as it has a strong control on the amount and intensity of solar radiation received at the land surface (Viereck et al., 1983;Morrissey and Strong, 1986).A difference in solar radiation be-tween north-and south-facing slopes supports the existence of permafrost on north-facing slopes and valley bottoms but not on south-facing slopes (Slaughter and Kane, 1979;Slaughter et al., 1983;Hinzman et al., 2006).Subsequently, permafrost-underlain and permafrost-free areas in this region display contrasting watershed characteristics and hydrological responses.The presence or absence of permafrost is the primary factor that creates a complex landscape mosaic of sharp spatial boundaries of contrasting vegetation cover and soil hydraulic and thermal properties, moisture dynamics, and water pathways (Bolton et al., 2000;Hinzman et al., 2002;Bolton, 2006;Petrone et al., 2006;White et al., 2008;Jones and Rinehart, 2010).Due to these small-scale complexities associated with permafrost distribution, simulating large-scale hydrological processes remains a challenge in the Interior Alaskan boreal forest.
The hydraulic conductivity of permafrost-affected soils is several orders of magnitude less than that of the overlying organic layers and the nearby permafrost-free soil (Rieger et al., 1972;Burt and Williams, 1976;Kane and Stein, 1983;Woo, 1986;Ping et al., 2005;Zhang et al., 2010a).The icerich permafrost is an impermeable layer at the permafrost surface that limits the hydraulic flow to the active layerthe thin, seasonally thawed soil layer above permafrost (Romanovsky and Osterkamp, 1995;Romanovsky et al., 2003).Streamflow in permafrost-dominated watersheds has been described as "flashy", responding rapidly to precipitation and snowmelt with storm hydrographs displaying a sharp rise and prolonged recession (Bolton et al., 2000;Quinton and Carey, 2008) and relatively low baseflow between precipitation events (Kane, 1980;Kane et al., 1981;Kane and Stein, 1983;Slaughter et al., 1983;Hinzman et al., 2002;Bolton, 2006;Petrone et al., 2006Petrone et al., , 2007)).On the other hand, in watersheds with no or relatively low areal permafrost extent, the soil hydraulic conductivity and infiltration capacity are much higher, resulting in a slower streamflow response to precipitation and snowmelt, relatively higher baseflow between storm events, and greater residence time of water in catchments (Bolton et al., 2000;Carey and Woo, 2001).
In the Alaskan sub-arctic, vegetation type, density, physiological, and structural properties such as leaf area index (LAI) and stomatal conductance display a strong variation between permafrost-dominated and permafrost-free soils (Viereck et al., 1983;Viereck and Van Cleve, 1984).These variations can lead to significant variation in the partitioning of precipitation and snowmelt into runoff and evapotranspiration (ET), and change in soil water content between permafrost-dominated and permafrost-free soils (Hinzman et al., 2002(Hinzman et al., , 2006;;Bolton et al., 2006;Cable et al., 2014;Young-Robertson et al., 2016).Permafrost-affected soils typically support coniferous vegetation that is shallowly rooted, has low physiological activity, is tolerant of cold and wet soils, and is able to survive a short growing season (Viereck et al., 1983;Viereck and Van Cleve, 1984;Morrissey and Strong, 1986;Mölders, 2011).In contrast, the well-drained, relatively warm, permafrost-free soils support deciduous vegetation that is much larger, has higher LAI and stomatal conductance, and has a deeper rooting network (Cable et al., 2014;Young-Robertson et al., 2016).This difference in vegetation between permafrost-dominated and permafrostfree soil can further influence streamflow responses (Naito and Cairns, 2011) due to the large differences in the rates of and controls on ET, particularly transpiration, canopy interception (Baldocchi et al., 2000;Ewers et al., 2005;Andreadis et al., 2009), and vegetation water storage (Young-Robertson et al., 2016).In Interior Alaska, deciduous trees have significantly higher transpiration rates and vegetation water storage compared to coniferous trees (Ewers et al., 2005;Yuan et al., 2009;Cable et al., 2014;Young-Robertson et al., 2016), which limits the availability of water for runoff in the permafrost-free watersheds.
Plot-scale and hill-slope studies have documented the differences in and relationships between soil thermal and hydrologic properties, and ecosystem vegetation composition in high-latitude cold regions (Dingman, 1973;Woo, 1976;Kirkby, 1978;Kane et al., 1981;Woo and Steer, 1983).However, larger-scale land surface parameterizations and the data products used in land surface, hydrological, and climate models do not adequately represent the complex sub-arctic watersheds with significant spatial variability in soil and vegetation dynamics.Hence, hydrological modeling using these coarse-resolution datasets cannot produce accurate estimates of the spatially variable and basinintegrated watershed responses.This is primarily because coarse-resolution datasets smoothed the non-linear smallscale watershed heterogeneities that control hydrological responses of the sub-arctic watersheds.Until the small-scale hydrologic processes, soil properties, and vegetation distributions are well represented, accurate large-scale hydrologic simulation and modeling remains extremely challenging (Walsh et al., 2005).
Several techniques, including parameter regionalization (Parajka et al., 2005;Pokhrel et al., 2008;Samaniego et al., 2010;Kumar et al., 2013), are proposed to address the limitation of large-scale hydrological simulation due to sub-grid variabilities of soil and vegetation cover properties.Standard regionalization (SR; Pokhrel et al., 2008;Troy et al., 2008) and multiscale parameter regionalization (MPR; Samaniego et al., 2010;Kumar et al., 2013) are widely used to parameterize the basin hydrological predictors.The latter approach (MPR) takes into account sub-grid heterogeneity of hydrological processes by implementing the prior predictors (the relationships between elevation, slope, vegetation characteristics, soil properties, etc.) into model parameterization (Samaniego et al., 2010;Kumar et al., 2013).However, the first approach (SR) lacks explicit representation of basin characteristics and primarily depends on calibration of subsets of the region (Troy et al., 2008).In this study, we proposed a sub-grid parameterization approach of soil property and vegetation characteristic parameterization of a mesoscale The permafrost map of CPCRW was produced by a small-scale observation across the watershed (Rieger et al., 1972;Yoshikawa et al., 2002).hydrological model.The sub-grid parameterization method proposed in this study closely fits the MPR approach.In both approaches, basin hydrological predictors are explicitly taken into account, and parameters are transferable to scales and locations other than a particular area of similar basin characteristics.
The primary objective of this study is to improve the prediction capability of hydrological models in Interior Alaska's boreal forest by implementing a small-scale parameterization scheme, which represents the spatial heterogeneity of soil properties and vegetation cover, into a mesoscale hydrological model.The study is conducted in two small subbasins of the CPCRW (LowP and HighP sub-basins of areas 5.2 and 5.7 km 2 , respectively) that are representative of the region.This study transfers the plot-and hill-slope-scale knowledge into a mesoscale distributed hydrological model (the variable infiltration capacity (VIC) model at a spatial resolution of approximately 1 km) so that its application can be extended to large basins in the region.A sub-grid parameterization method -on which the small-scale parameterization scheme is based -is derived from a high-resolution digital elevation model (DEM).Streamflow, ET and soil moisture simulations, based on the small-scale parameterization method and the coarse-resolution datasets, are presented and investigated in both sub-basins.

Study area
Located approximately 50 km northwest of Fairbanks, Alaska, the Caribou Poker Creek Research Watershed (CPCRW; centered at 65 • 10 N and 147 • 30 W; basin area ∼ 101 km 2 ) is within the zone of discontinuous permafrost and in the North American boreal forest region (Fig. 1).CPCRW is also within the Yukon-Tanana Uplands of the Northern Plateaus physiographic province (Wahrhaftig, 1965) with elevations ranging from 187 to 834 m.CPCRW has a long-term record of ecological, meteorological, and hydrological data (Bolton et al., 2000;Knudson and Hinzman, 2000;Hinzman et al., 2002Hinzman et al., , 2003;;Bolton, 2006).The climate of CPCRW is characterized as continental with large diurnal and annual temperature variations, low annual precipitation, low cloud cover, and low humidity (Haugen et al., 1982).The average annual air temperature range at CPCRW is very large (July mean temperature of 14.7 • C, January mean temperature of −18.4 • C) with a mean annual temperature of −2.1 • C (Fig. 2, Table 1).The mean annual precipitation is 417 mm, of which two-thirds occurs as rainfall (Bolton, 2006).
Seven soil series have been identified in CPCRW (Rieger et al., 1972).We group the seven series into two general categories: permafrost-dominated soils that are poorly drained with a thick organic layer, and permafrost-free soils that are well-drained with a shallow organic layer (Rieger et al., 1972).Permafrost in CPCRW is generally found along northfacing slopes and valley bottoms where the solar input is very low (Rieger et al., 1972;Haugen et al., 1982).Soils free of permafrost are generally found on south-to southwest-facing slopes.
The vegetation distribution in CPCRW displays a strong relationship with permafrost distribution.Coniferous vegetation that consists primarily of black spruce (Picea mariana) is generally found in areas underlain by permafrost.Feather moss (Hylocomium spp.), tussock tundra (mostly Carex aquatilis), and sphagnum mosses (Sphagnum sp.) are also found along the valley bottoms.Deciduous vegetation consists primarily of aspen (Populus tremuloides), birch (Betula papyrifera), alder (Alnus crispa), and sporadic patches of white spruce (Picea glauca), and is found on the welldrained, south-facing permafrost-free soils (Haugen et al., 1982).
There are several sub-basins in the CPCRW that differ in permafrost coverage and vegetation composition.The nearly permafrost-free or LowP (5.2 km 2 ) and permafrostdominated or HighP (5.7 km 2 ) sub-basins (Fig. 1) are selected for this study for their contrasting permafrost extents and vegetation cover compositions (Table 3).Unlike vegetation cover, permafrost extent, and solar input, local climate does not vary between the two sub-basins (Fig. 2, Table 1) due to the local understory convective mixing of the bulk atmosphere (Hinzman et al., 2006).The coniferous/deciduous vegetation composition, derived from Haugen et al. (1982), is approximately 30/70 % for the LowP sub-basin and 95/5 % for the HighP sub-basin (Fig. 3c, Table 3).The permafrost coverage is different between the two sub-basins, with approximately 2 and 53 % in the LowP and HighP sub-basins, respectively (Table 1; Rieger et al., 1972;Yoshikawa et al., 2002).

Sub-grid parameterization method
The sub-grid parameterization method is a simple method used to produce high-resolution soil property and vegetation cover maps for accurate representation of the watershed characteristics into the hydrological model.A fine-resolution DEM is primarily used to parameterize the sub-grid heterogeneity of the landscape.In addition, the relationships between vegetation, permafrost, slope, and aspect (Viereck et al., 1983;Morrissey and Strong, 1986;Hinzman et al., 2006) are included in the sub-grid parameterization method.
First, aspect of the landscape is calculated using the 30 m DEM (Aster, 2009) of the watershed with the ArcGIS aspect spatial analyst toolbox.The aspect toolbox uses elevation values of eight surrounding gird cells to calculate the gradient and aspect of each grid cell in the model domain (Burrough et al., 1998) resulting in an aspect map with nine classes (Fig. 3a) that are aggregated into two aspect classes: permafrost-underlain and permafrost-free aspects.The following assumptions are made in the classification of each aspect into permafrost-dominated and permafrost-free aspects: 1. North-facing slopes and valley bottoms are underlain by permafrost, and south-facing slopes are permafrost free (Rieger et al., 1972;Hinzman et al., 2002;Yoshikawa et al., 2002).
4. The hydraulic conductivity of frozen soil is 2 orders of magnitude less than the same soil in unfrozen conditions (Burt and Williams, 1976;Kane and Stein, 1983;Woo, 1986).
The small-scale observed (Rieger et al., 1972) and modeled (Yoshikawa et al., 2002) permafrost maps are used as a reference during the grouping of the nine aspect classes into permafrost-underlain and permafrost-free aspects.Finally, north, northeast, northwest, southwest, and flat aspects are classified as permafrost-underlain aspects with coniferous vegetation cover.South, southeast, east, and west are classified as permafrost-free aspects and with a deciduous vegetation cover.The vegetation cover and soil hydraulic property maps obtained from the sub-grid parameterization method are shown in Figs.3b and 4c, respectively.Since each grid cell has a unique soil property and the spatial scale or length of the observed permafrost map is on the order of few meters, a threshold by which a grid cell can be assumed 100 % permafrost or permafrost-free soil is required.After several rounds of testing, 0.5 grid cell fraction threshold is assumed to classify a grid cell as permafrost-underlain or permafrostfree soil.We used the small-scale observed (Rieger et al., 1972) permafrost map in the determination of this threshold value.As a result, a grid cell is assigned with permafrost soil property when the fraction of permafrost is greater than 0.5 (Fig. 4c).

VIC model description
The VIC model (Liang et al., 1994(Liang et al., , 1996;;Nijssen et al., 1997) is a mesoscale process-based distributed hydrological model that represents vegetation heterogeneity, multiple soil layers, variable infiltration, and non-linear baseflow.The version of the VIC model used in this study, VIC 4.1.2,contains several explicit formulations for snow accumulation and ablation, frozen soil, and evapotranspiration and solves the full energy balance or water balance at sub-daily or daily time steps to simulate the energy and water fluxes for individual grid cells.A separate routing model (Duband et al., 1993;Lohmann et al., 1996Lohmann et al., , 1998a, b) , b) is used to collect the runoff and baseflow simulations from each grid cell and simulate streamflow at the outlet of each watershed.VIC has been successfully used to simulate major hydrological and thermal processes at different spatial and temporal scales, such as basin streamflow (Abdulla et al., 1996;Abdulla and Lettenmaier, 1997a, b;Schnorbus et al., 2011;Bennett et al., 2012Bennett et al., , 2015)), evaporation and transpiration, canopy interception, soil moisture (Robock et al., 2003;Andreadis et al., 2005;Slater et al., 2007;Wu et al., 2007;Meng and Quiring, 2008;Billah and Goodall, 2012;Wu et al., 2015), surface and ground heat fluxes, ground temperature, and snowpack energy balance, including ablation and accumulation processes (Cherkauer and Lettenmaier, 1999;Cherkauer et al., 2003;Andreadis et al., 2009).In this study, a fully coupled water-energy balance mode of the VIC model, with 1/64 degree grid resolution (approximately 1 km), is used in the LowP and HighP sub-basins of the CPCRW.We utilize this higher model resolution in order to explicitly represent and simulate the contrasting hydrological responses of the valley bottoms and hill slopes that co-exist in a very short spatial distance.The frozen soil algorithm is activated for permafrost grid cells (Cherkauer and Lettenmaier, 1999;Cherkauer et al., 2003;Bowling et al., 2008) in order to solve the ground heat flux and surface energy balance.The frozen soil algorithm solves the thermal/heat flux equation through the soil column to calculate the frozen soil penetration and the ice content of each soil layer (Cherkauer and Lettenmaier, 1999).The frozen soil algorithm also solves the effect of frozen soil on moisture transport (Cherkauer and Lettenmaier, 1999).
VIC was forced with daily or sub-daily minimum and maximum air temperatures, precipitation, and wind speed.The remaining forcing data including incoming shortwave radiation, longwave radiation, atmospheric pressure, and rel-ative humidity are estimated based on the daily temperature range and precipitation as calculated by Mountain Microclimate Simulation Model (MTCLIM; Thornton and Running, 1999).MTCLIM is implemented within the VIC model (Bohn et al., 2013).The radiation estimate by the VIC model indicates an average of 15-20 Wm −2 higher shortwave radiation in the LowP sub-basin than the HighP sub-basin.In this study, we generate the daily gridded minimum and maximum temperatures, precipitation, and wind speed data from four meteorological stations located within the CPCRW (http://www.lter.uaf.edu/data.cfm)using the inverse distance weighting (IDW) interpolation method.
In addition to forcing, VIC requires parameters that describe soil properties, vegetation distribution and characteristics, and topographic information.In this study, a three-layer soil column is used with a top layer with a depth of 0.1 m and variable depth for the second and third layers based on model calibration.Soil properties are uniform within a grid cell, but these properties are allowed to vary for each layer in the grid cell.VIC requires several soil hydraulic and thermal property parameters.Due to a lack of high-spatial-resolution soil data in the region, most of the soil data -such as soil textural classes, saturated hydraulic conductivity, bulk density, and porosity -were extracted directly from the 5 arcmin (approximately 9 km) Food and Agriculture Organization digital soil map of the world (FAO, 1998).Other soil hydraulic properties -including field capacity wilting point, residual soil moisture content, water retention, and bubbling pressure -are estimated by the Brooks and Corey formulation (Rawls and Brakensiek, 1985;Saxton et al., 1986) based on soil texture classes, porosity, and bulk density information.All soil parameters are resampled to 1/64 degree using bilinear interpolation.Selected soil property values are shown in Table 2.
VIC uses a mosaic representation of vegetation coverage and subdivides each grid cell's land cover into a specified number of tiles.Each tile represents the fraction of the cell covered by a particular land cover type (coniferous, evergreen forest, grassland, etc.).The coarse-resolution vegetation cover composition data used in this study are obtained from the University of Alaska Fairbanks (UAF), Scenarios Network for Alaska and Arctic Planning (SNAP) 1 km × 1 km land cover map originating from the North American Land Change Monitoring System (NALCMS) 2005 dataset (http://www.snap.uaf.edu/data.php).Rooting depth data are extracted from the International Satellite Land Surface Climatology Project (ISLSCP) Initiative II Ecosystem Rooting Depths (Schenk and Jackson, 2009).Tree height, trunk ratio, displacement, and roughness of the primary vegetation classes are derived from field measurements by Young-Robertson et al. (2016).The remaining vegetation parameter values including LAI for each vegetation class in the region are derived from Myneni et al. (1997), Hansen et al. (2000), and Nijssen et al. (2001a, b).

VIC model parameterization methods
Based on sources of vegetation cover and soil hydraulic property data, we developed two parameterization methods: large-scale parameterization (large-scale, hereafter) and small-scale parameterization (small-scale, hereafter) schemes.The two schemes differ in vegetation cover, saturated hydraulic conductivity, and organic layer thickness.The remaining model inputs and parameters including meteorological forcing, vegetation characteristics, and soil properties are the same in both schemes.

Large-scale parameterization
The large-scale parameterization scheme is derived from the coarse-resolution SNAP vegetation cover (Fig. 3d) and FAO soil property (Fig. 4a) datasets as described in the previous section (Table 4).

Small-scale parameterization
The small-scale parameterization scheme consists of two small-scale parameterization sub-schemes: one derived from the sub-grid parameterization method (aspect parameterization, hereafter) and one derived from a permafrost map (permafrost parameterization, hereafter).In both small-scale parameterization schemes, the top layer of the permafrost cell is also assumed to be the organic layer (Ping et al., 2005;Zhang et al., 2010b) while mineral soil, as obtained from FAO dataset, is assumed for permafrost-free grid cells.

Aspect parameterization
The aspect parameterization is based on vegetation cover (Fig. 3b) and soil property (Fig. 4c) products of the sub-grid parameterization method (Table 4).Saturated hydraulic conductivity and organic layer thickness in the aspect parameterization depend on whether the grid cell is permafrost underlain or permafrost free (Fig. 4c).

Permafrost parameterization
The permafrost parameterization is primarily based on the small-scale observed (Rieger et al., 1972) and modeled (Yoshikawa et al., 2002) permafrost maps.In the permafrost parameterization, we assume the proportion of coniferous and deciduous vegetation in each grid cell is the same as the fraction of permafrost-underlain and permafrost-free soil, respectively (Table 4).Like the aspect parameterization, we classified a grid cell as containing permafrost when the fraction of permafrost in the grid cell is greater than 0.5 (Fig. 4b), as described in Sect.2.2.

Calibration and validation
Model calibration was performed by comparing streamflow simulation with observation for the period of 2001-2005 at LowP and HighP sub-basin outlets of CPCRW.The largescale parameterization scheme -FAO soil (Fig. 4a) and SNAP vegetation cover (Fig. 3d) datasets -is implemented during the calibration process.The multi-objective complex evolution (MOCOM) automated calibration approach developed by Yapo et al. (1998) was used to match the simulated and observed streamflow.MOCOM is a multi-criteria calibration approach based on random parameter sampling strategy to optimize several user-defined criteria (Yapo et al., 1998;Wagener et al., 2001).As suggested by Liang et al. (1994), the model was calibrated using baseflow generation parameters including maximum velocity of the baseflow (D smax ), infiltration parameter (bi), fraction of D smax where non-linear baseflow begins (D s ), maximum soil moisture for non-linear baseflow to occur (W s ), thickness of the second soil layer (D2), and thickness of the third soil layer (D3).Table 2 shows final calibration values of the user-defined parameters for both sub-basins.After the lumped sub-basin baseflow generation parameter values are derived by calibration, validation of the model with the large-scale and small-scale parameterization methods was conducted by comparing the observed and simulated runoff at the outlets of the LowP and HighP sub-basins from 2005 to 2008.Calibration results obtained from large-scale parameterization are directly implemented into small-scale parameterization during validation in order to examine how processes change with the new parameterizations without recalibration.We utilized two verification statistics: correlation coefficient (R 2 ; Eq. 1) and Nash-Sutcliffe efficiency (NSE; Eq. 2).
where N is equal to the number of data points (i.e., daily streamflow realizations), i is the time step (days), S is the simulated streamflow (mm day −1 ), and Q is the observed streamflow (mm day −1 ).R 2 (Eq. 1) describes the degree of linear correlation of simulated and observed runoff or goodness of fit.The value of R 2 ranges from 0.0 to 1.0, and larger values indicate better fit between simulation and observation.NSE (Eq.2) describes the relative magnitude of simulated runoff variances compared to variance in observed streamflow.NSE is also an indicator of model fit in terms of a scatter plot of simulated versus observed streamflow values, wherein a slope near the 1 : 1 line indicates a better fit.The value of NSE ranges from 1 (perfect fit) to −∞.While values larger than Table 2. Selected vegetation and soil parameters, their values and estimating methods or sources: runoff influence (RI) -runoff increases when parameter values increase (+) and runoff decreases when parameter values increase (−), coniferous vegetation (conif), deciduous vegetation (decid), permafrost-underlain soil (PF), permafrost-free soil (NPF), parameter values specifically for the LowP sub-basin (LowP), and parameter values specifically for the HighP sub-basin (HighP).

Parameters
Value ranges RI Sources

Vegetation parameters
Leaf area index: LAI 3.8-4.2(conif) and 0.11-6.0(decid) (−) (Hansen et al., 2000;Nijssen et al., 2001a 0.0 can be considered as acceptable levels of model performance (Krause et al., 2005;Schaefli and Gupta, 2007), values approaching 1.0 are more preferred depending on the study area.NSE uses of the mean observed value as a reference (Schaefli and Gupta, 2007).Hence, factors that affect the mean value of observed streamflow will have a stronger effect on the values NSE.In Interior Alaska, a lower value of NSE can be acceptable due to the large uncertainties of mean observed streamflow, which results from aufeis-related measurement errors at beginning of the snowmelt runoff season (Bolton, 2006).NSE values below zero indicates that the mean observed streamflow is better predictor than the simulated runoff (Krause et al., 2005).Percent difference (PD) is also used to indicate and visualize differences in simulated runoff from each parameterization scenario: where PD is the percent difference, Q R is reference observed streamflow, and Q i is simulated runoff.The range of PD (Eq. 3) falls between −100 and 100 % with 0 being the perfect match between observed and simulated streamflow.In addition to the three measures of success mentioned above, we compared the simulated and observed hydrographs to determine whether our sub-grid parameterization method can reproduce the strong relative contrast in runoff behavior between the two catchments.

Comparisons of vegetation cover in each parameterization scenario
Figure 3b, c, and d show the three vegetation cover representations of CPCRW derived from the sub-grid parameterization method (aspect), permafrost, and SNAP, respectively.The comparison of vegetation cover focuses on the coniferous and deciduous vegetation composition in the LowP and HighP sub-basins.Table 3 summarizes the proportion of coniferous and deciduous vegetation in the different parameterization schemes.The small-scale observed vegetation cover from Haugen et al. (1982;Fig. 3c) is used as a reference vegetation cover to evaluate the other vegetation cover scenarios.
The SNAP vegetation cover map represents the HighP sub-basin well but overestimates the coniferous proportion of the LowP sub-basin.It also shifts the dominant vegetation type from deciduous to coniferous for the LowP sub-basin (30 % in small-scale vegetation to 63 % from SNAP; see Table 3).Vegetation cover from the sub-grid parameterization method (aspect) captures the dominant vegetation type for both sub-basins (Table 3).Coniferous and deciduous distribution in the permafrost map parameterization is generally similar to permafrost-underlain and permafrost-free proportions.A 3/97 and 53/47 proportion of coniferous/deciduous vegetation is obtained from parameterization based on the permafrost map in the LowP and HighP sub-basins, respectively (Table 3)  5 summarizes the coefficient of determination (R 2 ) and NSE during the calibration of two sub-basins.Calibration results show that runoff simulation in the HighP sub-basin (Fig. 5b) is in better agreement than in the LowP sub-basin with the observed streamflow (Fig. 5a).Although the R 2 for the LowP (0.51) is greater than for the HighP (0.48), the NSE of the LowP is less than that of the HighP (0.17 for the LowP and 0.38 for the HighP sub-basin) indicating the peak flows are systematically underestimated in the LowP than in the HighP sub-basin (Table 5).The simulated streamflow for both subbasins is believed to be satisfactory given the inadequate soil and vegetation parameter representation from the coarseresolution datasets.Moreover, streamflow measurement during the early snowmelt period is mostly missing or not accurate most of the time due to the icing problem (Bolton, 2006); this contributes to the lower calibration values in this region.Figure 5c and d show the observed streamflow differences between LowP and HighP sub-basins for the calibration period of 2001-2005.The notable discharge contrast between the two sub-basins is shown by the flow duration curve in Fig. 5d.The peak flow in the HighP sub-basin is about 5 times higher than the peak flow in the LowP subbasin.Streamflow in the LowP sub-basin has lower peak flow, higher baseflow, and longer runoff response to precipitation and snowmelt than HighP.Final calibration parameter values suggest that the HighP sub-basin requires parameters that favor direct runoff while the LowP sub-basin requires parameters that favor more baseflow and infiltration (Table 2).

Hydrological fluxes under different parameterization methods
In order to test the improvement of VIC simulations in the small-scale parameterization over the coarser-resolution data products, three simulations were completed: one using largescale parameterization and two using small-scale parameterizations (aspect and permafrost).Streamflow, ET, and soil moisture were simulated in the LowP and HighP sub-basins from 2006 to 2008.Table 4 summarizes how the soil property and vegetation cover are parameterized in the three parameterization schemes.

Streamflow
The simulated hydrograph changes with vegetation cover and soil properties in both sub-basins.Figure 6a and c show the simulated versus observed streamflow of each parameterization scheme in the LowP and HighP sub-basins of CPCRW.
The percent difference (PD; Eq. 3) of each simulation from observation is also shown in Fig. 6b and d.Table 5 compares the three streamflow simulations against observations in both sub-basins.Streamflow simulation based on large-scale parameterization yields the lowest performance with a R 2 /NSE of 0.34/0.17and 0.48/0.42for the LowP and HighP subbasins, respectively (Table 5).The mean annual streamflow simulation in both sub-basins is largely underestimated under the large-scale parameterization scheme (by 45 to 68 and 27 to 52 % of the total annual runoff in the LowP and HighP sub-basins).The underestimation is higher in the LowP subbasin.Further, the simulated spring peak flow tends to occur earlier than observed peak for both sub-basins (Fig. 6).
The simulated streamflow hydrograph obtained from aspect parameterization -the small-scale parameterization scheme which is primarily derived from the sub-grid param-Table 3.Estimated percent of the sub-basins (LowP and HighP) covered by coniferous and deciduous vegetation, and percent of the subbasins underlain by permafrost.Values were obtained by applying different parameterization methods: SNAP vegetation cover (SNAP), a modified Haugen et al. (1982) vegetation cover map (modified Haugen), aspect-derived vegetation cover (aspect), the CPCRW permafrost map (permafrost), and large-scale FAO digital soil map of the world (FAO soil dataset).NA indicates that a given method was not utilized to obtain values for either the vegetation or permafrost distribution.eterization method -is closest to the observed streamflow hydrograph for both sub-basins (Fig. 6a and c) with an R 2 /NSE of 0.51/0.48and 0.53/0.52 for the LowP and HighP sub-basins, respectively (Table 5).Throughout the simulation period, the PD is also close to zero under the aspect parameterization compared to the other two parameterizations (Fig. 6b and d).Both peak and low flows in both subbasins are well simulated in the small-scale parameterization scheme (Fig. 6).Mean annual PD values are also close to zero for streamflow simulations using the aspect parameter-Table 4. Definition of the parameterization scenarios with respect to corresponding vegetation cover and soil property parameters.

Evapotranspiration
Unlike streamflow simulations, ET simulations generally do not display significant differences between the LowP and HighP sub-basins in all parameterizations.However, on an annual scale, there is a notable difference in the simulated ET between parameterization schemes in both sub-basins.
The difference in simulated mean annual ET between LowP and HighP sub-basins under aspect parameterization (356.5 and 335.2 mm in the LowP and HighP sub-basins) is almost half of the values of ET simulation when permafrost (396.6 and 352.2 mm in the LowP and HighP sub-basins) and largescale parameterization (402.8 and 362.8 mm in the LowP and HighP sub-basins) are used.
For further comparison between the LowP and HighP subbasins (Fig. 7), ET and streamflow simulations under the aspect parameterization are selected due to the better streamflow simulation performance under aspect parameterization (Fig. 6).Generally, ET is higher than runoff most of the time in both sub-basins.The HighP sub-basin, however, displays more runoff than ET during snowmelt, late summer/fall, and large storm events (Fig. 7b).

Soil moisture
Variation in soil properties and vegetation cover is generally found to be sensitive to VIC's soil moisture simulation.We also note that the rates of snowmelt, melt/refreeze, and soil moisture increase/decrease between the parameterization schemes are notably different between the two subbasins (Fig. 8).Unfrozen soil moisture content over the winter does not show any significant variation between the two sub-basins (Fig. 8).However, the change in storage in the LowP sub-basin is higher than in the HighP sub-basin.Additionally, the fluctuations of soil moisture and change in storage with the snowmelt and rainfall events are higher in the LowP sub-basin as compared to the fluctuation in the HighP sub-basin (Fig. 8).These results are consistent with previous observations (Young-Robertson et al., 2016) and modeling (Bolton, 2006) studies in the Interior Alaskan boreal forest ecosystem.These results indicate that soil moisture in the LowP sub-basin is more sensitive to climate forcing and land surface representation as compared to the wet system of the HighP sub-basin.

Discussion
One of the important findings from our small-scale subgrid parameterization method, as shown by Figs. 3 and 4, is that an acceptable spatial vegetation cover and soil property representation of the Interior Alaskan sub-arctic can be produced without any direct measurements or remote sensing methods.Large-scale vegetation cover and soil property products that are typically used in mesoscale hydrological modeling do not accurately represent spatial heterogeneities of the Interior Alaskan boreal forest ecosystem.This is because of the fact that coarse-resolution data products include regions outside the watershed area that influence the average value of the grid cell, rather than capturing the spatial heterogeneity of ecosystem properties within the watershed boundaries.This study indicates that previously documented relationships among soil hydraulic and thermal properties, vegetation cover, topography, slope, and aspect (Viereck et al., 1983;Viereck and Van Cleve, 1984;Morrissey and Strong, 1986;Hinzman et al., 2006;Mölders, 2011) can be used to formulate a methodology by which local to regional scale landscape properties can be incorporated into mesoscale hydrological models.
The soil hydraulic properties obtained from the largescale FAO soil data do not reflect any variation between permafrost-affected and permafrost-free soils (Fig. 4a).However, many surface and sub-surface soil hydraulic properties including the hydraulic conductivity and thickness of the organic layer show a large difference between permafrost-affected and permafrost-free soil (Rieger et al., 1972;Burt and Williams, 1976;Kane and Stein, 1983;Woo, 1986;Ping et al., 2005;Zhang et al., 2010a).Therefore, hydrological modeling using such data creates large uncertainty that cannot be easily corrected with model calibration, as the hydrology and landscape properties including the vegetation composition (Viereck et al., 1983;Viereck and Van Cleve, 1984;Morrissey and Strong, 1986;Hinzman et al., 2006;Mölders, 2011) are primarily controlled by the presence or absence of permafrost (Kane, 1980;Kane et al., 1981;Kane and Stein, 1983;Slaughter et al., 1983;Hinzman et al., 2002;Bolton, 2006;Petrone et al., 2006Petrone et al., , 2007)).
Our multi-objective model calibration indicates that watersheds with different permafrost proportions display consistent and contrasting parameter values between nearly permafrost-free and permafrost-dominated watersheds (Fig. 5, Table 2).As reported by previous studies (Bolton et al., 2000;Bolton, 2006;Hinzman et al., 2002), simulated runoff in permafrost-free watersheds is best fitted with observations for parameter values that favor high infiltration and baseflow and a longer recession period (Fig. 5a, Table 2).On the other hand, the best fit between simulated and observed streamflow in the permafrost-dominated watersheds (Fig. 5b) is achieved only when baseflow parameters that favor more direct runoff and inhibit infiltration (Table 2) are introduced into the model.During model calibration, especially for the heavily parameterized distributed models, it is important to consider these variabilities in the baseflow parameters of the discontinuous permafrost watersheds.It  could save a considerable amount of time and resources, especially for big watersheds and limited computational facilities.
Based on the three model evaluation criteria (Eqs.1-3) and the comparison of runoff behaviors between simulation and observation, streamflow simulations in both the permafrost-free/LowP (Fig. 6a and b) and permafrost-dominated/HighP (Fig. 6c and d) methods perform better when the small-scale sub-grid soil hydraulic and thermal properties (Fig. 4c) and vegetation cover (Fig. 3b) heterogeneity are introduced into the VIC model soil property and vegetation cover parameterizations.Comparison of the improvement between the LowP and HighP sub-basins indicates that the strongest streamflow simulation improvement (Tables 3 and 5) was observed in sub-basins that are highly unrepresented by the large-scale data products -the LowP sub-basin (Figs. 3, 4, and 6a and b).In general, our study shows that most of the peak and low flows in both sub-basins are captured well with the smallscale parameterization scheme as compared to the direct use of coarse-resolution land surface data products (Fig. 6, Table 5), except the 2007 spring flooding event (Fig. 6).The exception for the 2007 spring runoff peak can be explained by the high rainfall and mixed precipitation in the beginning of the 2006-2007 winter (nadp.isws.illinois.edu?/data?/sites?/siteDetails.aspx??id=AK01?&net?=NTN).Rainfall at the onset of winter generally freezes on the surface soil layer over the winter and limits snowmelt from infiltrating during the spring snowmelt period to generate flooding.
Although we did not detect significant differences in the ET simulations between parameterization methods, we found that sub-basin average ET in the LowP sub-basin is higher than in the HighP sub-basin, except at the beginning of snowmelt when deciduous trees of the LowP are storing snowmelt water but not yet transpiring (Young-Robertson et al., 2016;Figs. 7 and 9).Unlike other snowdominated regions in middle latitudes, where ET tends to have a strong relationship with precipitation and a change in storage (Lohmann et al., 1998b, a), ET in this region is more strongly related to temperature than precipitation or change in soil moisture storage.This apparent decoupling between ET and precipitation is likely because the deciduous trees in this region are utilizing the snowmelt water stored in their trunks rather than being directly tied to rainfall (Young-Robertson et al., 2016).This implies that temperature is the limiting factor for ET in this region.This may have important implications for the general warming trend on a broader spectrum, as in the case of permafrost degradation (Romanovsky andOsterkamp, 1995, 2000;Romanovsky et al., 2002;O'Donnell et al., 2009).
Several studies have documented that VIC soil moisture simulation is strongly sensitive to how vegetation cover (Ford and Quiring, 2013;Tesemma et al., 2015) and soil properties (Liang et al., 1996;Lohmann et al., 1998b;Lee et al., 2011;Billah and Goodall, 2012;Ford and Quiring, 2013) are represented.This study shows that vegetation cover and soil property representation are not only sensitive to the soil moisture content simulations, but they also have a strong influence on the rate of snowmelt, and snowmelt and refreeze processes in the Alaskan sub-arctic environment.Small-scale parameterization schemes are the best at capturing the expected soil moisture pattern in both the LowP and HighP subbasins.This includes high soil moisture in the permafrostdominated sub-basin due to lower hydraulic conductivity of the soil and the faster winter freeze-up of the soil column in the permafrost-dominated sub-basin than in the permafrostfree one.
We also made an effort to understand how the annual water balance is partitioned into runoff, ET, and change in soil moisture in the permafrost-free/LowP and permafrostdominated/HighP sub-basins.Both sub-basins do not differ in the percentage of each component between the large-scale and small-scale parameterizations.As shown in the previous studies (Hinzman et al., 2002;Bolton, 2006;Young-Robertson et al., 2016), our results indicate that in both sub-basins ET dominates the annual water balances (Figs. 8  and 9), while the change in storage is the lowest in both sub-basins.There are, however, variations in the percentage and pattern of the water balance components between LowP and HighP sub-basins.The runoff from LowP sub-basin is lower than that from the HighP sub-basin due to the high tree water storage and transpiration (Cable et al., 2014;Young-Robertson et al., 2016) of the deciduous trees and higher infiltration of the permafrost-free soil of the LowP sub-basin than HighP sub-basin.The mean annual change in storage in the LowP sub-basin is about 35 % higher than in the HighP sub-basin.From this result, we can conclude that the soil thermal and hydraulic properties dictate the partitioning of water into the different processes in this region.This is the first study that introduces vegetation and soil property heterogeneity from high-resolution topographic information into the mesoscale hydrological model.The results need to be verified for regional-scale basins to determine the transferability of our findings to similar areas where the land surface data products do not adequately represent the spatial heterogeneity for accurate hydrological simulations.

Conclusions
This study indicates that coarse-resolution soil and vegetation data products -data that are used extensively in land surface modeling in the middle and lower latitudes -do not adequately represent the North American boreal forest discontinuous permafrost ecosystems.Hydrological modeling with coarse-resolution products cannot adequately simulate important small-scale processes.This is because small-scale permafrost distribution and ecosystem composition primarily control the land surface processes in this region.This indicates the need for landscape modeling or sub-grid parameterization that can produce these small-scale features and incorporate them into land surface models.The strong correlative relationship between topography, vegetation, and permafrost distribution (Burt and Williams, 1976;Haugen et al., 1982;Viereck et al., 1983;Hinzman et al., 1991Hinzman et al., , 1998) ) in this region can be used to produce a sub-grid parameterization method that represents the small-scale soil property and vegetation cover heterogeneity for distributed hydrological models.
This study also shows that our sub-grid parameterization method -based primarily on this strong relationship between permafrost, topography, and vegetation composition -produced a better representation of permafrost and vegetation cover than large-scale soil and vegetation cover products.Hydrological simulations -including basin-integrated and spatially variable runoff, evapotranspiration, and soil moisture dynamics -using a small-scale parameterization scheme derived from a sub-grid parameterization method are a notable improvement over parameterizations based on coarseresolution data products.Finally, in our effort to demonstrate methodologies that can improve hydrological modeling through a small-scale parameterization method, we intend to implement and test results from this pilot study into a large river basin in the next phase of the research.

Figure 1 .
Figure 1.Location and permafrost extent (light blue shading) in the Caribou Poker Creek Research Watershed, Alaska, and the two subbasins of interest (low permafrost or LowP sub-basin and high permafrost or HighP sub-basin).The permafrost map of CPCRW was produced by a small-scale observation across the watershed(Rieger et al., 1972;Yoshikawa et al., 2002).

Figure 2 .
Figure 2. Climatology of the Caribou Poker Creek Research Watershed, Alaska, from 1970 to 2010: (a) mean annual precipitation (mm), (b) mean annual air temperature ( • C), (c) mean January air temperature ( • C), and (d) mean July air temperature ( • C).These climate datasets are extracted from the Alaska-Pacific River Forecast Center (APRFC; Bennett, 2014).Only US National Climate Data Center cooperative daily climate stations were used to generate the gridded fields of 800 m resolution using an inverse distance weighting function.Note that these data are not used to force the model in this study.In this study, data from four very close observation stations are used.

Table 1 .
Mean (1970Mean ( -2012) )  climatology of the Caribou Poker Creek Research Watershed (CPCRW), Alaska, and the sub-basins (LowP, HighP): MAP (mean annual precipitation, mm), MAT (mean annual air temperature, • C), MIN (mean minimum air temperature, • C), MAX (mean maximum air temperature, • C), MJan (mean January air temperature, • C), and MJuly (mean July air temperature, • C).Basin/sub-basin Permafrost MAP MAT MIN Figure 3. (a) Aspect map -derived from the DEM -of the CPCRW.Coniferous and deciduous vegetation composition at the CPCRW as derived from (b) model-based aspect or topography of the watershed, (c) an observation-based vegetation coverage map modified from Haugen et al. (1982), and (d) a large-scale-based Scenario Network for Alaska and Arctic Planning (SNAP) vegetation coverage map is shown.Note that panel (b) is prepared from the aspect map (a) and the sub-grid parameterization method discussed in Sect.2.2.

Figure 4 .
Figure 4. Saturated hydrologic conductivity (mm day −1 ) as derived from (a) a large-scale FAO soil dataset, (b) a permafrost map of Rieger et al. (1972), and (c) an aspect map of CPCRW is shown.Panels (b) and (c) are small-scale parameterizations obtained from altering largescale values according to the presence or absence of permafrost.Grid cells that are classified as permafrost soil are assigned hydraulic conductivity 2 orders of magnitude less than the value obtained from the large-scale (FAO) dataset.FAO soils map do not indicate any of the permafrost-underlain soil hydraulic properties which are known to be present in the CPCRW.

Figure
Figure 5a and b show daily simulated versus observed streamflow for the calibration period of 2001-2005 in the LowP and HighP sub-basins, respectively.Table5summarizes the coefficient of determination (R 2 ) and NSE during the calibration of two sub-basins.Calibration results show that runoff simulation in the HighP sub-basin (Fig.5b) is in better agreement than in the LowP sub-basin with the observed streamflow (Fig.5a).Although the R 2 for the LowP (0.51) is greater than for the HighP (0.48), the NSE of the LowP is less than that of the HighP (0.17 for the LowP and 0.38 for the HighP sub-basin) indicating the peak flows are systematically underestimated in the LowP than in the HighP

VegetationFigure 5 .
Figure 5. Observed versus simulated streamflow during the calibration period of 2001-2005: (a) at the LowP sub-basin and (b) at the HighP sub-basin of the CPCRW.The discharge contrasts between HighP and LowP sub-basins are shown in panels (c) and (d).Panel (c) denotes the streamflow hydrograph and (d) denotes the flow duration curve in both sub-basins.

Figure 6 .
Figure 6.Comparison of streamflow simulations -obtained from different parameterizations -with observations (a) in the LowP sub-basin and (c) in the HighP sub-basin.Panels (b) and (d) show the percent difference (PD) for each runoff simulation from the observed streamflow in the LowP and HighP sub-basins, respectively.

Figure 7 .
Figure 7.Comparison of evapotranspiration (ET) and runoff simulations -obtained from parameterization based on aspect parameterization (a) in the LowP sub-basin and (b) in the HighP sub-basin.

Figure 8 .
Figure 8. Partitioning of precipitation or snowmelt (P , precipitation) into evapotranspiration (ET), runoff, and change in soil moisture storage ( S) in the (a) LowP and (b) HighP sub-basins.All simulations are under aspect parameterization.

Figure 9 .
Figure 9. VIC-model-simulated (based on aspect parameterization) mean annual (2006-2008) percentage of water balance components, P (black), ET (red), and runoff (green) in the (a) LowP and (b) HighP sub-basins.The change in storage component is very small compared with the other components; hence, it is not included.

Table 5 .
Coefficient of determination (R 2 ) and Nash-Sutcliffe efficiency (NSE) values for the LowP and HighP sub-basins for calibration and validation periods.