Using field observations to inform thermal hydrology models of permafrost dynamics with ATS ( v 0 . 83 )

Climate change is profoundly transforming the carbon-rich Arctic tundra landscape, potentially moving it from a carbon sink to a carbon source by increasing the thickness of soil that thaws on a seasonal basis. However, the modeling capability and precise parameterizations of the physical characteristics needed to estimate projected active layer thickness (ALT) are limited in Earth system models (ESMs). In particular, discrepancies in spatial scale between field measurements and Earth system models challenge validation and parameterization of hydrothermal models. A recently developed surface–subsurface model for permafrost thermal hydrology, the Advanced Terrestrial Simulator (ATS), is used in combination with field measurements to achieve the goals of constructing a process-rich model based on plausible parameters and to identify fine-scale controls of ALT in ice-wedge polygon tundra in Barrow, Alaska. An iterative model refinement procedure that cycles between borehole temperature and snow cover measurements and simulations functions to evaluate and parameterize different model processes necessary to simulate freeze–thaw processes and ALT formation. After model refinement and calibration, reasonable matches between simulated and measured soil temperatures are obtained, with the largest errors occurring during early summer above ice wedges (e.g., troughs). The results suggest that properly constructed and calibrated onedimensional thermal hydrology models have the potential to provide reasonable representation of the subsurface thermal response and can be used to infer model input parameters and process representations. The models for soil thermal conductivity and snow distribution were found to be the most sensitive process representations. However, information on lateral flow and snowpack evolution might be needed to constrain model representations of surface hydrology and snow depth.


Introduction
In Arctic tundra, the thickness of the soil layer that reaches above 0 • C, defined as the active layer thickness (ALT), largely determines the volume of carbon stores available for decomposition.Predicting ALT is therefore critical when characterizing potential climate feedbacks due to greenhouse gas release into the atmosphere from decomposition of organic soil carbon (McGuire et al., 2009;Koven et al., 2011;Schneider von Deimling et al., 2012).Current long-term predictions of ALT generally use large-scale Earth system models (ESMs) with simplified representations of the hydrothermal processes, and are thus producing results with significant uncertainty (Schaefer et al., 2009;Slater and Lawrence, 2014;Koven et al., 2013).The freeze-thaw dynamics that determines the ALT functions on a vertical scale of centimeters and vary horizontally on a scale of meters across the characteristic microtopography of polygonal tundra (Painter et al., 2013).Freeze-thaw dynamics is also strongly controlled by a local inundation state (Muster et al., 2012), which can vary over a horizontal extent of meters to hundreds of meters.These local-scale to intermediate-scale processes are under-resolved or completely missing in ESMs.Improved fine-scale simulation capabilities can inform the representation of soil thermal processes in regional to global-scale models by identifying appropriate representations of key processes governing ALT, and by providing calibrated model parameterization.
Previous efforts have been made to characterize ALT using field, lab, and numerical experiments (e.g., Osterkamp and Romanovsky, 1996;Romanovsky and Osterkamp, 1997).Site-specific properties of Arctic soils, such as porosity, bulk thermal conductivity, and water retention characteristics, have been measured in lab settings from samples taken in the field (Hinzman et al., 1991;Letts et al., 2000).Those field and lab measured properties were then used in ESMs in order to predict future ALT and permafrost conditions (Beringer et al., 2001;Lawrence and Slater, 2008;Subin et al., 2013).However, such regional-and global-scale projections are difficult to constrain by measurements of soil properties made at vastly smaller scales of observation.This scale gap between the governing fine-scale physical processes and large-scale simulations impedes direct model validation against measurements, which has motivated development of fine-to intermediate-scale hydrothermal models (e.g., Hinzman et al., 1998;Hansson et al., 2004;Daanen et al., 2007;Mckenzie et al., 2007;Painter, 2011;Karra et al., 2014;Endrizzi et al., 2014;Yi et al., 2014); for a review see Kurylyk and Watanabe (2013).Numerical experiments using high-resolution coupled hydrothermal models, which are calibrated against fine-scale measurements, can play a fundamental role in understanding the governing physical processes of ALT formation.
Simulating thermal hydrology in polygonal tundra systems is a challenging endeavor that requires simultaneous representation of multiple physical processes including phase change and highly nonlinear constitutive relationships (e.g., Painter, 2011).Soil thermal conductivity alone depends on volumetric water content, mineral composition, porosity, density, and temperature (Farouki, 1981).In soils experiencing freeze-thaw cycles, the phase of water also affects bulk thermal conduction (e.g., Johansen, 1977;Peters-Lidard et al., 1998).Latent heat of fusion and evaporation impart further control on the propagation of the freezing front and therefore thermal conduction.Thermally driven vapor transport can slowly change ice content and thus thermal conduction in partially and fully frozen soils (Grimm and Painter, 2009;Karra et al., 2014).Characterizing subsurface properties for modeling is further complicated due to variability in microtopography and cryoturbated soil that create a heterogeneous surface and subsurface in polygonal tundra systems.In addition, coupling of the soil to the atmosphere involves a balance among multiple energy transfer processes, which occur across interfaces of snow, water, ice and exposed ground.All of the above attributes describing soil structure, surface energy balances, and processes of phase change result in a tightly coupled hydrothermal system.Therefore, numerical experiments using high-fidelity representations of fine-scale processes require calibrated parameters that are able to effectively link dependent processes.
Despite the model gains of calibrating thermal properties (Romanovsky and Osterkamp, 1997;Nicolsky et al., 2009), relatively few hydrothermal modeling studies of Arctic systems have documented calibration procedures, with the noted exception of Tang and Zhuang (2011) and Jiang et al. (2012).Additionally, correct model structure representation, capable of representing the system based on known physical relationships while using plausible model parameters, is typically not known a priori.Calibration of a model with an inadequate model structure may result in over-fitting and unreliable forward simulations that incorrectly predict system behavior based on faulty process representation (e.g., Beven, 2005;Gupta et al., 2012).Therefore, when dealing with a coupled system of complex processes, it is imperative that the conceptual model is refined during the calibration process to increase model structure adequacy (Gupta et al., 2012).
Iterative modeling approaches that use repeated model runs with different combinations of parameters, governing mechanisms, or process representation can help fundamental system understanding (Clark et al., 2008;Kavetski and Fenicia, 2011;Fenicia et al., 2011;Larsen et al., 2014).Here we use an iterative procedure that integrates finely resolved models with field observations and measurements to develop a process-rich model with physical mechanisms and parameters consistent with measurements from the Department of Energy Office of Science -Next Generation Ecosystem Experiment (NGEE-Arctic) site, Barrow Environmental Observatory (BEO), Barrow, Alaska (Fig. 1).The iterative process of using field observations to inform model development and subsequent simulations to inform new data needs is referred to here as the model-observation/experiment or ModEx cycle (Fig. 2).Clearly, there is no unique way to approach iterative modeling procedures (Larsen et al., 2014), which is intrinsically subjective and highly dependent on expert knowledge.Well-documented examples of successful applications of model refinement are thus invaluable for building the required experience base.We use repeated calibration of model parameters against site-specific field measurements and iterative model adjustments of the model structure to reduce mismatch between model predictions and measurements in order to attain a viable model of thermal hydrological conditions.
In this paper we summarize our ModEx experience involving the detailed use of subsurface temperature and snow cover field data to develop and test process-rich simulations of ALT dynamics, such that observational data and necessary physical dynamics are incorporated into the model.In order to calibrate and refine model structure in a tractable fashion, the model development first focuses on a series of subsurface-only calibrations in Sect. 3 before moving onto a series of coupled surface energy balance and subsurface calibrations in Sect. 4. The end result is a set of calibrated thermal and hydrological parameters for moss, peat, and mineral soil layers, along with a consistent model structure,  employed for various microtopographic positions characteristic of polygonal tundra.We demonstrate how the detailed calibration and model development effort informs understanding of the key processes that define the ALT in polygonal ground.We further complete the ModEx cycle by discussing how future data needs can reduce system uncertainty and refine our understanding of process behavior.

ModEx process applied to thermal hydrology processes in permafrost
Our variant of the ModEx approach is shown schematically in Fig. 2. Starting with site identification and characterization, field observations and measurements begin to form the modeling activity by providing model parameter inputs and targets for the model calibration process.Standard model calibration -denoted by the inner loop -aims to match simulations to field measurements by varying parameters while keeping the model structure fixed.Here the ModEx procedure moves beyond the standard calibration by assuming the model itself is uncertain, but can be further constrained through successive comparison to observations (outer loop in Fig. 2).These improved model runs then inform the observation process by specifying the data needs, either through further calibration or through informal numerical experimentation.Such model refinement is not a unique process, and can be achieved through multiple avenues.For example, flexible modeling approaches have been used in understand structural errors by combining functional aspects of several models (Clark et al., 2008;Kavetski and Fenicia, 2011;Fenicia et al., 2011).We implement ModEx model refinement by evaluating the plausibility of calibrated parameters in addition to the mismatch between field measurements and simulated responses.
The calibration process uses a multi-dimensional response surface to evaluate the plausibility of parameters and the degree of mismatch between simulated results and observed data.Sets of parameter values are mapped to the response surface with the respective mismatch between simulated results and field observations/measurements, quantified by the root mean squared error (RMSE), which determines the shape of the responses surface.RMSE is given by where θ is a vector comprised of a combination of parameter values, Ti (θ ) is the ith simulated temperature given θ , T i is the ith calibration measured temperature target, and N is the number of calibration targets.Simulations with a poor fit to data have high RMSE and a corresponding high value on the response surface.Conversely, simulations with a good fit to data have a low RMSE and therefore a low value on the response surface and may constitute a minimum in the response surface.A minimum in the response surface indicates that a possible calibration has been achieved.However, in the case of a complex model with high dimensionality, multiple local minima may exist, which causes gradient-based calibrations to find non-unique solutions (Beven, 2006).Model structure error can also cause the response surface to slope to a parameter boundary, indicating that over-fitting is necessary to calibrate to observed data (Beven, 2005).Therefore, it is important to extend calibration boundaries beyond the acceptable parameter range to allow the optimization algorithm to travel into the infeasible range when the response surface dictates an implausible combination of parameter values, indicating an inadequate model.By altering the model itself, and not just model parameters, the ModEx process can work to reduce model structure error and reshape the response surface such that the simulated system matches the observed data and calibrated parameters are realistic.The ModEx process is facilitated by two software components.First, for calibrating a given model to determine an optimal match to measurements we use the parameter estimation software, PEST (Doherty, 2004), which implements the Levenberg-Marquardt algorithm (Marquardt, 1963).This method uses gradient descent to determine (from a highdimensional space of calibration parameters) a set of parameters that (in a local sense) minimize the forward model's error in predicting observed data.Second, the ModEx process requires iterative exchange, comparison, and addition of process models, which is greatly facilitated by a dynamically configured model with many process options.Therefore, a framework that manages complexity and allows for rapid development of new physical representations is critical.To this end, we have implemented the Advanced Terrestrial Simulator (ATS), version 0.83, as a collection of physics modules managed by the Arcos multiphysics framework (Coon et al., 2015b).At runtime, Arcos dynamically forms a dependency graph where each variable identifies its data requirements, allowing the automation of model evaluation.Process kernels (i.e., a single PDE (partial differential equation), such as mass balance) are coupled to form complex systems of equations in which each term or component can easily be replaced.The ease of swapping and adding processes makes model verification and evaluation more tractable, and facilitates the ModEx process by allowing the model structure to be easily changed and extended.

Site description and initial conceptual model set-up
The lowland, cold continuous permafrost tundra at BEO was established as the end-member of the NGEE-Arctic sites, which follow a bioclimatic gradient that extends to the warm discontinuous permafrost, shrub tundra environment of the Seward Peninsula.The site supports the NGEE-Arctic goal to improve climate model predictions through advanced understanding of coupled processes in Arctic terrestrial ecosystems.NGEE-Arctic scientists are collecting multiscale in situ field measurements and remote-sensing observations of polygonal tundra.A range of polygon types including low center polygons, which are surrounded by rims and, in some areas, shallow troughs, and high center polygons with deep troughs as a result of ice-wedge degradation.The focus of the model development chronicled here is NGEE-Arctic site "area C" (Fig. 1), which is characterized by ∼ 50 cm deep troughs, rims and shallow low centers.The site was chosen because it serves as a representative state that polygonal tundra may develop into as permafrost degrades.Three one-dimensional (1-D) model domains represent the main ice-wedge polygon sub-features: center, rim, and trough.Each domain includes a unique model structure and parameterization (Figs. 1  and 3).Nine soil temperature sensors (0.1 to 1.5 m depth) from three soil profiles representing center, rim, and trough were used to compare simulated to measured soil temperatures (http://lapland.gi.alaska.edu/vdv/vdv_historical.php?station_id=20&page_id=-1&direct=1).The shallowest soil temperature sensor (2 cm depth), located just under a layer of green moss, provided the subsurface model with an upperboundary condition.Each column had unique near-surface soil temperature forcing, measurements for calibration, and assigned peat layer thicknesses typical of the microtopographical features.The center, rim and trough columns had an organic peat layer of 10, 6 and 14 cm respectively.The underlying mineral soil was a silty loam to a total depth of 50 m.A far field bottom boundary condition was held constant at −6 • C to represent the average deep permafrost temperature in the North Slope of Alaska (Romanovsky et al., 2010).All columns were initialized by first freezing the entire column from the bottom with a no flux upper-boundary condition and then spunup to a cyclical steady state using a "decadal average" year of daily values looped for 20 simulation years.The decadal average year was made by averaging the daily mean temperature from 1 October 1998 to 30 September 2009 at Barrow, AK, for each day of the year to produce forcing data that represented seasonal trends.Each calibration parameter combination was then simulated for an additional year using the same decadal average year before the in situ soil temperature forcing data at 2 cm depth were applied.

Model description
The ATS solves water and energy flow in variably saturated soils at temperatures above and below freezing using the conservation equations described by Karra et al. (2014) (see also Painter, 2011 andCoon et al., 2015a).Liquid and ice partitioning is represented by the model of Painter and Karra (2014).In this model liquid water can coexist with ice below 0 • C, as is well known (e.g., Miller, 1980;Williams and Smith, 1991), which occurs due to soil-surface forces and pore geometry.Ice-water partitioning is related to the soil water characteristic curve under unfrozen conditions.Thus, soil moisture characteristic curve parameters directly contribute to thermal conduction regimes when the soil is saturated and frozen.Two variations of a three-phase thermal conductivity model (Painter, 2011), both an extension of Johansen (1977), were used to relate bulk thermal conductivity to ice and liquid contents.The three-phase thermal conductivity model is described in detail in Appendix A. The first thermal conductivity model variant is a simplification of the Johansen method and is referred to as the bulk phase component (BPC) model.The BPC model has porosity and the bulk-phase unfrozen saturated thermal conductivity (K sat,uf ) and bulk-phase dry thermal conductivity (K dry ) as input parameters to be calibrated (Eq.A3 in Appendix A).The third bulk-phase component, saturated frozen thermal conductivity (K sat,f ) (Eq.A3) is then calculated based on an empirical relationship with K sat,uf shown in Eq. (A8) in Appendix A. The second option for thermal conductivity is denoted the material component (MC) model.The MC model has porosity and the solid material thermal conductivity K soil as an input parameter; K sat,uf and K dry are then calculated using functional relationships shown in Eqs. ( A6) and (A11), respectively.Material components ice, water, and gas are fixed material thermal conductivities in the MC model.Switching from the BPC model to the MC model reduces the dimensionality of parameter space by one.Perhaps more importantly, the MC model calculates all bulk-phase components as a function of soil porosity; thus, porosity is more correlated to thermal conductivity in the MC model as compared to the BPC model.

Parameter starting values and ranges from literature
Parameter value ranges for moss, peat, and mineral soils of Arctic tundra systems were drawn from literature and field observations at the NGEE-Arctic site (NGEE-Arctic data portal; http://ngee-arctic.ornl.gov;see references in Appendix C).Estimates of reasonable calibration ranges are listed in Table 1.Depending on the thermal model being calibrated, seven to eight parameters for both peat and mineral soil were calibrated creating a 14-16-dimensional parameter space.Based on the literature and assigning greater weight to study sites with characteristics and proximity to Barrow, AK, a probable parameter guess was selected as one starting point of the calibration process, along with seven additional starting calibration parameter sets located near the boundary of parameter space.Together the eight starting calibration parameter sets determined the dependence of calibration results on starting location (i.e., the degree of non-uniqueness in the calibration results).

ModEx applied to the subsurface system
Our experience with the ModEx cycle applied to the coupled subsurface hydrothermal system at the BEO is shown in process flow form in Fig. 4. In this cycle the ATS model only included subsurface processes, and the shallowest measurement of temperature (2 cm depth) was used as a timedependent upper-boundary condition to force the model.Measurements at deeper locations (from 0.1 to 1.5 m) (Fig. 3) represented the calibration targets.In the initial iteration, calibration was performed using the BPC model for thermal conductivity and assumed full saturation of the soil column.That calibration resulted in parameters being out of range.
In the second iteration, the thermal conductivity model was changed to an alternative model (the MC model), which 1.3-2.821.3-1.90.1-0.33Residual VWC [-] 0.02-0.180.04-0.220.05-0.18resulted in improved parameter values but inferior match to measured soil temperatures.In the final iteration, surface pressure was calibrated at the borehole locations, which determines liquid saturation that affects near-surface thermal conductivity.The iteration to calibrate surface pressure resulted in a calibration that was judged to be adequate for continuation of a coupled surface energy balance-subsurface calibration and model development (see Sect. 4).Details of the subsurface calibration and model development are discussed in the remainder of this section.

Subsurface BPC vs. MC thermal model
The first subsurface calibration attempt used the BPC model (Fig. 4) and resulted in unrealistic parameters sets.The response surface of the center and rim columns resulted in calibrated peat porosities to move to the lower parameter boundary (Fig. 5).With a few exceptions, the thermal conductivities for peat in the center, rim, and trough calibrated outside the acceptable parameter range to the lower boundary for peat.The first calibration iteration produced unrealistic parameter values and indicated that the BPC model is not an adequate calibration tool for subsurface hydrothermal modeling.
In the second iteration of our model-data integration cycle, subsurface thermal conductivity was simulated using the MC model instead of the BPC model, which reshaped the calibration response surface such that calibrated porosities spread out across parameter space and away from the parameter boundary.Calibrating with the MC model generally kept the porosity parameters within the acceptable range and improved the thermal conductivity parameters; however, RMSE increased for all columns (Table 2).Yet, the MC model was selected for the remainder of the paper because calibrated parameters were reasonable.

Simultaneous calibration of center, rim, and trough
Upscaled parameters for larger-scale models were calibrated by coupling all three columns to find a single set of peat and mineral soil hydrothermal parameters.The calibration was coupled by combining objective function results from each microtopographical feature in the PEST Levenberg-Marquardt algorithm to inform the next parameter update that is then applied to all 1-D columns.The initial application of the coupled calibration resulted in unrealistic parameter values and motivated a reformulation of the conceptual model to include near-surface unsaturated conditions necessary for center and trough simulations.The saturated condition response surface decreased the K e for the peat layer, and maintained or increased heat conduction for mineral soil.Peat porosity and peat K sat,uf calibrated to the lower calibration boundary of 0.59 and 0.33 W m −1 K, respectively, and mineral porosity calibrated to a higher value (0.65) than the peat porosity, while the mineral K sat,uf calibrated to 1.04 W m −1 K.An unsaturated near surface could conversely result in a reduced thermal conductivity for the

Variably saturated versus unsaturated soils
The fourth iteration of the ModEx cycle allowed the surface pressure to be a calibration parameter for the center and trough columns, which were previously assumed fully saturated for the duration of the year.A surface pressure less than atmospheric results in an unsaturated condition at the top of the soil column, and introduces air with low thermal conduction, creating a gradient of increasing K e with depth.
The surface pressure in the rim, which did not manifest the issues described above, was still fixed at 25 % gas saturation.It is important to note that calibrating a top pressure for this set of subsurface calibrations does not allow the nearsurface saturation to vary throughout the year and; therefore, the saturation state is only a function of pressure and ice content.Figure 6 illustrates how K e of peat decreases with lower surface pressure.Decreasing surface pressure results in decreased K e , but the effect is especially large during the winter.Ice has a large thermal conductivity compared to either water or gas; any variation in the amount of ice in the domain will cause a large change in K e .
The eight calibration starting locations for the uncoupled column calibration were then re-tested for the center and trough by calibrating surface pressures (Fig. 4).Here we only tested unsaturated conditions using the MC thermal model rather than posthumously retesting prior model structural decisions, as the MC model was thought to be more physically accurate.The new conceptual model with unsaturated conditions at the soil surface became the second model refinement, which resulted in a reshaped parameter response surface.More calibrated center porosity values were within the acceptable parameter range when surface pressures were calibrated, but more trough peat porosities calibrated to the upper peat boundary.Both the center and trough had more calibrated K dry, material within the realistic range.The increase in calibrations resulting in porosities outside their acceptable range for the trough may be indicative of the trough being more saturated than the center, or being fully saturated.However, unsaturated conditions reduced the RMSE for both the center and trough indicating a better model fit (Table 2).The increased model fit with more realistic parameters suggests that it is necessary to capture characteristic saturation states of the dominant topographical features (center, rim, and trough) to constrain model calibration.Furthermore, the single coupled center-rim-trough calibration, where surface pressures were calibrated, also resulted in realistic parameters with surface pressures at 95 440.9 and 97 638.2 Pa for the center and trough, respectively (Table 5).Moreover, the revised coupled calibration found a low RMSE of 0.554 • C and the temperature time-series results fit measured data near the point of the active layer depth (Fig. 7).

Surface methods
After the calibration of subsurface thermal properties, a 2 cm moss layer was added to each of the three columns and a surface energy balance model was used to calibrate both the thermal properties of the moss layer and parameter values for the surface energy balance in a second set of ModEx iterations (Fig. 8).Parameters from the subsurface calibration were used in the coupled snow-surface energy balance-subsurface simulation.The ranges of hydrothermal parameters for moss are listed in Table 1.The surface energy balance, described in detail in Appendix B, is implicitly coupled with subsurface thermal hydrology and is based on the work of Hinzman et al. (1998) and Ling and Zhang (2004).Simulated snow deformation and snow density changes described by Eqs. ( B6) and (B7) in Appendix B are applied on a single layer snowpack.The center, rim, trough columns had unique maximum head boundary conditions of 8, 0.7, and 15 cm, respectively, where water spills off each column at or above the specified head heights.The maximum head boundary conditions were selected according to relative elevation differences observed in polygonal tundra.
For the surface energy balance calibration each column was spunup over a 10-year loop using decadally averaged air temperature along with shortwave radiation, relative humidity, and wind speed data from 1 October 1998 to 30 September 2009 at Barrow, AK, where meteorological data from each day in the 10 years were averaged together.After spinup, daily meteorological data from 2010 to 2013 were used to drive the model.This forcing data were compiled from several sources; the incoming solar radiation is from the Atmospheric Radiation Measurement (ARM) Climate Research Facility (ARM, 1993(ARM, , 1996)); rainfall and snowfall is from Barrow airport (Station GHND:USW00027502 National Weather Service, National Atmospheric and Oceanic Administration); air temperature, relative humidity and wind speed are from individual research projects at the BEO (Liljedahl et al., 2011;Zona et al., 2014); landscape-averaged end-of-winter snow depth is from the Circumpolar Active Layer Monitoring (CLAM) Program (Shiklomanov et al., 2012).Daily rain-and snowfall were adjusted for undercatch according to Yang et al. (1998).A second adjustment was applied to the snowfall where the average ratio between the 1997-2006 CALM observations and the undercatch-adjusted and measured soil temperature profile (black).The biggest difference between initial temperature profiles and the calibrated profiles is the wintertime temperature for each column, and is a result of distributing snow on the center, rim and trough and depth hoar representation.Snow distribution also had the greatest control in the ALT (Table 4).
National Weather Service (NWS) snow accumulation was applied to respective daily precipitation events.The simulation results from 2013 were then compared with measured subsurface temperature data, at a 2 cm depth below the moss layer.The runtime increased when including the surface energy balance component model such that automated calibration algorithms could no longer be employed.Manual calibration was used with 2 cm soil temperature borehole measurements and observed ALT, as calibration targets.

ModEx applied to the coupled surface energy balance system
The second set of ModEx cycle iterations is presented in Fig. 8 in process flow form.The focus of the second set of ModEx cycles is process identification and calibration of the moss layer and surface energy balance parameters.The first iteration of the cycle coupled the surface energy balance model and 2 cm moss layer to the previously calibrated and refined subsurface model.The initial iteration matched surface temperatures well in all three columns; however, soil temperatures were generally under simulated for center and trough columns, especially during winter (Fig. 9).The second iteration added a microtopography-informed snow depth from measurements between universal transverse mercator (UTM) coordinates: Northing 7910330-7910350 and Easting 585900-585930, which encompasses the borehole temperature locations.temperatures substantially improved, which also resulted in late summer ALT to be in or near the observed ALT range.However, near-surface winter rim temperatures were colder than measured because microtopography-informed snow distribution produces less snow on rims and results in less snow cover insulation.The third iteration of the ModEx cycle added a depth hoar representation in the snowpack, which resulted in a better representation of winter rim soil temperatures and caused the rim ALT to be within the range of observed ALT.In the final ModEx iteration hydrothermal properties of moss and surface energy balance parameters were hand calibrated within the plausible range of parameters space, which resulted in only slight improvements of near-surface temperature simulations.Details of how each iteration of the ModEx cycle (for the coupled surface energy balance -subsurface model) informed both model development and future data needs are presented below.

Importance of surface energy balance governing saturation time series
Forcing the subsurface thermal propagation through a surface energy balance in the second set of ModEx cycles attempts to capture variable surface thermal conductivities due to changing surface saturation states as pulses of precipitation enter the subsurface and subsequently dry from evaporation.Modeling studies that do not explicitly model surface energy balance processes may not adequately capture near-surface saturation states and have reported the greatest error during the summer when highly variable soil moisture states occur (Romanovsky and Osterkamp, 1997;Jiang et al., 2012).It is known that soil moisture influences soil temperature in addition to meteorological controls, by governing the amount of latent heat of fusion necessary to freeze-thaw and evaporate water from soils (Johansen, 1977;Farouki, 1981;Peters-Lidard et al., 1998;Subin et al., 2013).Consequently, the timing of the precipitation pulses and subsequent drying may have a significant impact on ALT because the highly variable saturation states coincide with summer soil warming.Therefore, the second set of ModEx cycles starts with a more detailed representation of transient soil moisture conditions, which is the third major model refinement.Simulation results showed that it is important to capture the freeze-up timing with the highly variable fall saturation state in order to set up near-surface ice content and thermal conductivity during winter (Fig. 10a).Properly representing the freezeup with transient soil moisture is especially important giving that winter has the largest range of possible thermal conductivity values (Fig. 6) and therefore is highly variable from year to year.Simulating the surface energy balance for each column resulted in varied model fits to the measured 2 cm soil temperature time series.For example, the simulated center and trough 2 cm soil temperature during the summer is consistently lower than the measured 2 cm temperature (Fig. 9, center and trough plots), especially for the early summer, which in turn lowers the simulated soil temperature at depth.However, simulated 2 cm deep soil temperatures for the rim matched measured soil temperatures.The ability for the model to match measured summer surface temperatures for the rim versus the center and trough is most likely attributed to either the spatial differences and local microtopography of the three columns and/or the surface saturation state.The rim is higher and therefore drier than the center and trough columns (Fig. 3).To mimic microtopographical differences in the three columns, unique maximum ponded water depths were assigned to each column; the rim had a negligible maximum ponded depth with effectively no standing water from snowmelt compared to the center and trough columns.Unfortunately, limitations to our surrogate 1-D model exist and inherently contribute to model structural error.For example, the largest deviation of surface temperature for the trough occurred during the fall as the temperature dropped below freezing.The measured surface temperature at 2 cm depth had a longer duration of the zero curtain, where soil temperatures are at 0 • C as water freezes, compared to the simulated surface temperature (Fig. 9).One possible explanation for this difference is that there is greater soil moisture in the trough than was simulated, as added soil moisture will extend the time to freeze a block of soil.A possible reason for the underestimated soil moisture is that the 1-D surrogate model neglected lateral surface and subsurface flow that could be flowing on to the column, especially for troughs that are connected to an extensive trough-network.Monitoring of lateral flow in polygonal tundra systems could help to constrain the conceptual model needed to understand soil moisture dynamics.

Snow model refinement
The largest gains from calibrating the surface energy balance portion of the model came from the fourth model refinement, which resulted from two additional ModEx iterations (1) updating the conceptual and numerical model to add snow depth variation informed by microtopography and (2) including a depth hoar representation in the snowpack model.The snowpack at Barrow, AK, is scoured relatively flat due to strong winds (Benson and Sturm, 1993;Zhang et al., 1996) resulting in deeper snow in depressions such as troughs and low centers.To match measured snow depths of the three topographical features (Table 3), snowfall was in-creased for the center and trough columns by 30 % (3.6 cm) and 82.5 % (9.9 cm), respectively, and reduced for the rim to 87 % (10.4 cm) of the total adjusted snowfall (12 cm) for the snow year of 2012-2013.Although manually distributing snow does not fully capture snowpack dynamics, especially year-to-year snowpack variation, simulated nearsurface (2 cm) winter temperature more accurately matched the measured temperatures (Fig. 9, center and trough plots).Summer ALT increased for both the center and trough, which improved the model prediction to be within the observed ALT range for the trough and closer to the observed ALT range for the center column (Table 4).Conversely, the decreased snow depth over the rim cooled the winter surface soil temperature below the measured soil temperatures.Including a depth hoar layer in the model counteracted the reduced insulation of a shallower snowpack on the rim.The combination of reduced snow depth and depth hoar representation on the rim translated to a slightly shallower ALT, resulting in the rim ALT to be within the observed ALT range.Without snow re-distribution or depth hoar representation the snowpack evolved to a density of 410 to 440 kg m −3 by mid-May and early June as determined from Eq. (B26).At first, this seemed reasonable because the surface of tundra snow forms a wind slab layer due to the wind scouring affect with densities between 400 and 500 kg m −3 (Benson and Sturm, 1993;Dominé et al., 2002).Having a snowpack surface with high densities is required to accurately capture snow surface albedo.However, underneath the wind slab layer, a hoar layer forms during the winter with a density between 100 and 250 kg m −3 (Benson and Sturm, 1993;Zhang et al., 1996;Zhang, 2005), which reduces the thermal conductivity of the snowpack.The single layer snow model did not include the formation of a depth hoar layer and would overestimate the thermal conduction of the snowpack and therefore increase winter cooling of the ground surface.The iterative ModEx process, however, encouraged us to formulate a way of both representing snowpack top densities in order to properly simulate surface albedo, and capture a depth hoar layer to account for lower snowpack thermal conduction.The new formulation, similar to the snow classes used by Schaefer et al. (2009) and Sturm et al. (1995), employed in the model runs plotted in Fig. 9 conduction by assuming a depth hoar layer forms for 15 % of the snowpack with a calibrated density.Then a harmonic mean snow density is taken between the depth hoar layer and rest of the snowpack in order to calculate an adjusted thermal conductivity of the snowpack.Because this process applies only to calculating the snowpack thermal conduction, the simulation of snow albedo is unaffected.Center and rim depth hoar densities calibrated to 110 kg m −3 and the trough depth hoar density calibrated to 190 kg m −3 .The addition of the depth hoar also reduced end of winter (2 May) snowpack densities from above 400 kg m −3 to between 320 and 370 kg m −3 (Table 3), which is closer to the measured endof-winter average snowpack density of 326 kg m −3 .Adjusting the snow accumulation due to topographically informed snow distribution and including a depth hoar representation increased the insulative effect of the snowpack and had a clear impact on winter near-surface temperatures (Fig. 9).In addition snow distribution and depth hoar representation improved summertime ALT predictions (Table 4).Summertime changes in ALT due to winter conditions highlights a memory trait of the system and the necessity to capture dominant winter processes in order to simulate transient thermal conditions in physically based models.Research by Hinkel and Hurd (2006) showed that large snow drifts cause long-term deepening of the ALT, due in part from the additional insulation for the snow and the loss of cold thermal propagation into the subsurface.Timing of snowpack accumulation and thickness has also been shown to govern permafrost formation (Zhang, 2005).However, at the scale of microtopographical relief, where trough to rim vertical relief changes by 40 cm within a horizontal distance of a meter, questions regarding how snow thickness and associated meltwater inputs affect ALT formation remain.Results for this work show that topographically informed snow distribution will change the spring and early summer surface saturation state (Fig. 10d) due to distributed snow water equivalence amounts (Table 3).The change in early summer surface saturation state then affects the thermal conduction for early sum-mer as well as adding greater water mass that then requires a greater amount of energy to heat up (Hinkel and Hurd, 2006).Moreover, studies have found that the depth hoar layer can be as thick as 50 % of the snowpack height in arctic conditions (Sturm et al., 1995;Schaefer et al., 2009).However, due to continuous wind slab and depth hoar formation significant snowpack heterogeneities develop within and across topographical features (Sturm and Benson, 2004;Sturm et al., 2004).Therefore, spatially distributed snow depth measurements and snowpack density profiles that characterize local snowpack variability and over microtopographical features can help constrain both modeled snowpack thermal conduction representation, and surface water inputs.

Surface energy balance calibration
In the final ModEx iteration and model refinement, attempts to increase the simulated summer surface (2 cm) temperature were made (Fig. 8).Special attention was paid to the early summer wet conditions found in the center and trough for the Julian dates between 150 and 200 (Fig. 10b and d), where the biggest error in surface temperatures is found (Fig. 9 center and rim plots).It was thought that by calibrating parameters that control the amount of energy entering the subsurface under wet conditions, such as the albedo of standing water (see Appendix B for details), the surface temperature of the center and trough, which are wet, will increase without affecting the relatively dry rim surface temperature.However, variables specific to the surface energy balance and moss properties had little effect of simulated soil temperature during the snow free summer.The range of accepted albedo values for tundra varied from 0.12 to 0.17 based on wet or dry conditions (Grenfell and Perovich, 2004), and the albedo range for standing water values ranged from 0.11 to 0.20 for the months of May through September for latitude of 70 • near Barrow, AK (Cogley, 1979).Only slight gains in simulated surface temperature were observed by decreasing albedo of standing water from 0.14 to 0.11 and tundra from 0.15 to 0.12.This iteration of the ModEx cycle shows that adjusted standing water albedo and roughness length within the perceived parameter range did not substantially improve model fit, which suggest that the model is lacking either a necessary process representation or the calibration parameter range is not correct.One possible improvement would be a distributed surface albedo representation that provides a unique albedo for centers, rims, and troughs.Local-scale tundra albedo measurements can inform models of spatially distributed albedo conditions.Another possible explanation is how atmospheric mixing coefficients such as roughness length (noted as z 0 in Eq.B12 in Appendix B) could change over microtopographical features.Specific exchange coefficients for each microtopographical feature would then produce unique sensible and latent heat fluxes.For example, rim surface temperatures were well matched under current roughness lengths.But topographically protected troughs and centers could have a different roughness length, which may result in changes to latent and sensible heat exchanges and higher surface temperatures.Observations of how microtopography affect near-surface wind and associated atmospheric mixing could support an improved conceptualization of sensible and latent heat exchanges.

Summary and conclusions
1-D thermal hydrology models of transient saturation and frozen states combined with a surface energy balance model were used to represent active layer dynamics in polygonal tundra at the Barrow Environmental Observatory.In the coupled model, surface water was allowed to pond to a specified maximum height but any additional water was removed (spill over condition).The surface model also includes a surface energy balance model for bare, snow-, ice-or watercovered ground.The model was used in combination with borehole temperature and snowpack field measurements in an iterative model-data integration (ModEx) framework to produce calibrated model parameters and refine constitutive models and process representations.The particular variant of the ModEx approach combined calibration with iterative refinement of the model structure; parameter feasibility and model-observation mismatch were used as metrics to achieve the objective of model development and identification of viable representations of key thermal hydrological processes.
The results demonstrate the effectiveness of using borehole temperature measurements to effectively develop and refine the model structure for hydrothermal models of permafrost-affected landscapes.Results also suggest that properly constructed and calibrated 1-D models coupled to a surface energy balance may be adequate for representing thermal response at a given location provided the maximum ponded depth (spill point) is known for that location.This suggests a multiscale modeling strategy that uses overland flow models to establish the spill point (maximum ponded depth) at each surface location in conjunction with a set of thermal hydrology simulations.Further evaluations of the 1-D representations against 3-D model representations are needed to identify additional process representation and the appropriate level of model complexity to capture scale dependencies of thermal dynamics.In addition, it is important to note that the largest discrepancy between model and field measurements occurred during early summer in the troughs and that mismatch is likely indicating model structural error with inflow of water from upstream locations and/or unique surface energy balance conditions.Observations of water fluxes such as evapotranspiration, lateral flow, and snowmelt at the sub-polygon scale would help model representation and, in particular, the role of advective lateral heat transport.However, the temperature mismatch was brief and confined to the trough location, and is thus not expected to have large consequences for integrated results such as thaw depth.
The model refinement process identified the representation of thermal conductivity -specifically the dependence of bulk thermal conductivity on porosity, water content and ice content -as a constitutive model that affects model performance.Thus, field and laboratory work to better constrain hydrothermal representation and the governing model parameters would help reduce uncertainty in model projections.Further modeling efforts focusing on uncertainty analysis and environmental parameter sensitivity can provide information regarding which parameters govern model outcome and thus inform future observational efforts.Similarly, snowpack properties and snow distribution were found to be important.Investigations similar to Benson and Sturm (1993), Zhang et al. (1996) and Tape et al. (2010) that better define the relationship between depth hoar, microtopography and wind slab formation would help reduce uncertainty in projections.For example, snowpack dynamics and density profile observations at the NGEE-Arctic site will inform models of how the snowpack develops and how snow will distribute across microtopography.
More generally, these results demonstrated the utility of one particular approach to merging observations and models in environmental applications.In this particular iterative approach, formal parameter estimation methods are used iteratively.Each calibration run -the inner loop in Fig. 2 -minimizes mismatch between data and models with fixed model structure.The "reasonableness" or feasibility of the calibrated parameters and the RMSE are performance metrics for the calibrated model.Model structural adjustment, the outer loop in Fig. 2, is initiated when calibrated parameters fall outside reasonable bounds.Although structural model adjustments were done in an ad hoc manner guided by experience and knowledge of the system being modeled, the resulting refinements have produced robust representation of system response.Such an approach combining structural model adjustments drawing from literature, field observations and for- The α of ponded water is the average α of standing water at a latitude of 70 • from May to September.During freezing and thawing of the ground surface any ponded water is subdivided into an unfrozen water fraction and a frozen water fraction in ATS.The α values for this surface are then an average of water and ice α values and are found to transition linearly between the two states (Grenfell and Perovich, 2004) based on unfrozen water fraction.Transitional α values between each type of surface can occur and are triggered when the snowpack height is less then 2 cm, or the standing water height is less then 10 cm.The transition height for ponded water is based on the penetration depth of shortwave radiation in ice ( 10 In this model, if snow is present it is always the top surface, and ponded water or surface ice will always be below snow and above the tundra surface.Therefore, the α value is set first by snow, if present, then by standing water and/or ice if present, and finally by the tundra surface. Once the incoming radiation components of the energy balance model are computed, evaporative resistance (E r ) is then calculated by where the air resistance term (R air ) is the inverse of the turbulent exchange of latent and sensible heat (D eh ) and the sta-bility function (ζ ): κ is the von Karman Constant 0.41 [-], U s is the wind speed at the reference height (z r ) of the meteorological measurement location, and z 0 is the roughness length.Due to the changing conditions of the landscape at barrow, z 0 changes from 0.005 [m] for wind swept snow (Wieringa and Rudel, 2002), to 0.04 [m] for polygonal tundra (Weller and Holmgren, 1974;Hansen, 1993).
The stability function (ζ ) accounts for both stable (ζ stable ) and unstable (ζ unstable ) atmospheric conditions (Price and Dunne, 1976) ζ unstable conditions occur when the ground surface (T s ) is warmer than the air temperature (T a ) causing more air to mix vertically.R i defines atmospheric stability; where R i is positive in the stable condition and R i is negative in an unstable condition.
g is the acceleration due to gravity.R soil [m s −1 ] is calculated following the methods used by Sakaguchi and Zeng (2009) and is only implemented during ground surface evaporation when the saturation state of the upper most subsurface cell adjacent to the domain surface is less than 1.
The exponent b in Eq. ( B16) is a Clapp and Hornberger (1978) fitting parameter for the soil water characteristic curve, assumed to be 1 for moss (Beringer et al., 2001), which covers the tundra surface and is simulated as the top subsurface layer for the tundra.
L is dry layer thickness or the length vapor must travel from the point of evaporation.Once all necessary components of the energy balance are calculated, either the snow energy balance or surface energy balance is computed.The snow energy balance, Eq. (B1), is calculated if snow height (Z s ) is more than 2 cm.The ground surface energy balance, Eq. (B2), is used if no snow is present.Between Z s of 0 and 2 cm, a transition between the snow energy balance and the ground surface energy balance is used where both surface conditions are solved.When calculating the energy balance for the transitional regime, the snow energy balance assumes a Z s of 2 cm for all components that depend on Z s and an area-weighted average is used between the ground surface and snow energy balance based on the actual Z s that is equal to or less than 2 cm.Assuming a 2 cm Z s within the snow energy balance calculation prevents unreasonable heat conduction through the snowpack (Q c ), calculated by where k s is the effective thermal conductivity of snow [W m −1 K] and is calculated from an empirical function of ρ s used by Ling and Zhang (2004), described by Goodrich (1982) k s = 2.9 × 10 −6 ρ 2 s . (B19) The snow and surface energy balance use the same formulation for Q h and Q Out lw .Q h is where ρ a is the density of air 1.275 kg m −3 , and C p is the specific heat of air (1004 J K −1 kg).Q Out lw is where ε s is the emissivity of the surface.The ε s for snow and ice 0.98 [-] is taken from Liston and Hall (1995), and the ε s for tundra is 0.92 (Ling and Zhang, 2004) and for standing water is 0.979 (Robinson and Davies, 1972).Q e is slightly different between the snow and ground surface energy balance where the porosity (ϕ s ) of the top cell in the ground surface is included for the surface energy balance calculation.where E r is the evaporation resistance as defined by Eq. (B8) and R soil is 0 in the case of snow, or condensation on the surface.L s is the latent heat of sublimation for snow (2 834 000 J kg −1 ) and L e is the latent heat of evaporation for the ground surface (2 497 848 J kg −1 ).e s is the vapor pressure of the snow or surface, and A pa is the atmospheric pressure (101.325kPa).
Once the energy balance is calculated, then the water fluxes to the ground surface are calculated.In the case of snow, if the snow surface temperature (T s ) is greater than freezing, T s is set to freezing and the snow surface energy balance is recalculated with all excess energy assigned to the melting energy (Q m ), and a melting rate (Mr) [m s −1 ] is calculated from where ρ w is the density of water and H f is the heat of fusion for melting snow 333 500 J kg −1 .Condensation or sublimation of the snow surface is also calculated from Q e , where the sublimation/condensation rate (S r ) is added to the total water flux.If T a and Z s > 0 and S r is positive, then Sublimation is removed from the snowpack when S r is positive.If only the ground surface energy balance is used then water is delivered to the ground surface as precipitation and condensation when S r is negative.Water is evaporated from the surface-subsurface when S r is positive.Snow water equivalence (SWE), Z s , and ρ s are tracked through the simulation of snowpack evolution and related by: Both Z s and ρ s are important in the snow energy balance equation for calculated Q c and snow α, and both variables evolve as the snowpack ages through snowpack deformation simulated by (Martinec, 1977) ρ settled = ρ freshsnow SP age 0.3 , (B26) where ρ freshsnow is assigned a density of 100 kg m −3 , SP age is the age of the snowpack.The total snowpack density and Z s are then calculated by a weighted average of three components: old settled snow, new snow accumulation, and any ice from condensation.The density of condensation is assigned 200 kg m −3 .

Figure 1 .
Figure 1.Lidar of site-C with the three observation locations mapped and the greater Barrow, AK, area (credit Garrett Altmann).

Figure 2 .
Figure 2. Schematic representation of a Model Observation/Experiment (ModEx) process involving traditional parameter estimation-calibration (inner loop) and model structuralconceptual refinement (outer loop).Observations inform simulation input and provide a starting point for a conceptual model.Both the conceptual and numerical model is then tested against observations.In successive ModEx iterations the model is then refined and at times re-drawn in order to elicit governing processes that shape model outcome to match observed and measured phenomena.Finally, model experiments and the identification of governing processes inform future observations as to which measurements are needed to assess the state of the system.

Figure 3 .
Figure 3. Diagram of the three 1-D columns and the associated measured soil temperature depths.

Figure 4 .
Figure 4.The ModEx cycle as applied here to subsurface thermal hydrologic system in freezing-thawing soils.

Figure 5 .
Figure 5. Panels (a), (b) and (c) show center, trough and rim, respectively, calibrated peat and mineral porosities from eight calibrations starts.Panels (d), (e) and (f) show calibrated saturated unfrozen thermal conductivities (K sat,uf ) for peat and mineral soil layers from the same eight calibrations starts.K sat,uf values from the MC calibration are calculated from Eq. (3).Blue diamonds used the BPC model for soil thermal conductivity, red squares used the MC model for soil thermal conductivity, and green triangles added surface pressures as a free calibration parameter to the MC model for soil thermal conductivity.Color-coded asterisks represent the average calibrated parameter for each model tested for the eight calibration starts, but are not actual calibrated results.Accepted parameter space delineated from literature and site observations in all cases are mapped as clear areas.Shaded areas are the calibration space outside of the acceptable parameter space.This figure shows how the calibration response surface changes as the model changed from (1) BPC to (2) MC to (3) unsaturated.

Figure 6 .
Figure 6.Thermal conductivity of peat throughout a year with different surface pressures.Percent liquid saturation is based off of summer time water liquid saturation, which changes during winter due to an increase in ice saturation.The change in thermal conductivity coincides with spring thaw, approximately Julian day 160 or early June, and fall freeze-up near Julian day 265 or late September.

Figure 7 .
Figure7.The subsurface un-calibrated and calibrated temperature time series is compared to measured soil temperature time series to showcase the improvement from the calibration process at 40 cm depth for the center, trough and rim.The initial un-calibrated parameters were selected from the literature search described in Sect.2.4 and Appendix C. Calibration fit to observation varies from the three columns, but shows marked improvement from initial un-calibrated time series and are most accurate for all three during the summer at depth where active layer thickness is delineated.

Figure 8 .
Figure 8.The ModEx cycle applied to the surface energy balance and moss parameters.

Figure 9 .
Figure 9. Temperature profiles for a 2 cm depth are shown for the center (a), rim (b) and trough (c), using the initial surface energy balance parameters (blue), calibrated surface energy balance (red)and measured soil temperature profile (black).The biggest difference between initial temperature profiles and the calibrated profiles is the wintertime temperature for each column, and is a result of distributing snow on the center, rim and trough and depth hoar representation.Snow distribution also had the greatest control in the ALT (Table4).

Figure 10 .
Figure 10.Ice and liquid saturation are shown in (a) for the simulated years of 2010-2013 at 2 cm depth along with bulk thermal conductivity for a center column.Notice that ice saturation and thermal conductivity during the winters are unique for each simulation year.(b) is a detailed view of year 2013 of ice and liquid saturation and the bulk thermal conductivity for the center.(c) and (d) show the corresponding ice and liquid saturations for the trough and rim, along with the respective thermal conductivities for the 2 cm soil depth for the year 2013.(b)-(d) have unique ice and liquid saturation and therefore bulk thermal conductivity for each column, which is a result of both the maximum ponded depth for each column and the snow distribution that mimics wind scouring of the snow surface at Barrow, AK.

−
Tran snow ] Tran tundra = [1 − Tran snow ] − Tran water , (B8) where Z is the height of water or snow and Pen is the penetration depth of shortwave radiation.The transitional α value is then calculated by α trans = α snow Tran Snow + α water Tran water + α tundra Tran tundra .(B9)

Q
e,snow = ρ a L s E r 0.622 e a − e s A pa Q e,ground_surface = φ s ρ a L e E r 0.622 e a − e s A pa , (B22)

Table 1 .
Valid parameter range for calibration sets.

Table 2 .
The calibration error from the measured values reported as the RMSE • C (phi) increased between the (1) BPC model to the (2) MC saturated model.Thus there was greater error in the model results, but the calibrated parameters were more realistic.Phi then decreased between the (2) MC saturated model and (3) the MC unsaturated model.
peat layer while maintaining thermal conduction for the mineral soil layer.

Table 3 .
Measured snow depth ranges were gathered from a compilation of 258 snow depth measurements taken on 2 May 2013 in the area encompassing all three borehole temperature measurements.universal transverse mercator (UTM) coordinates: Northing 7910330-7910350 and Easting 585900-585930.Measured snow water equivalence (SWE) ranges were calculated from measured snow depth and the measured average snowpack density of 326 kg m −3 .All simulated values were taken on simulation day 2 May 2013.

Table 4 .
The ALT for all three columns are listed for each iteration of the calibration process, also with the range of possible ALT from the observed data.The observed ALT range was made by finding the deepest borehole measurement for center rim and trough with a temperature above 0 • C for at least a day and the shallowest borehole measurement with all temperatures below 0 • C.

Table 5 .
Final calibrated parameter table (referred to throughout the text).

Table B1 .
Albedo values and parameter range.