MODELS, DATASETS, AND METHODS

MESH modelling framework

The model utilized here is the Modélisation Environmentale Surface et Hydrologie model (MESH: Pietroniro et al., 2007). MESH is a physically-based, semi-distributed modelling system with three main components: (1) the vertical processes of moisture and heat flux land-atmosphere transfers, represented either by the Canadian Land Surface Scheme (CLASS: Verseghy, 1991, 2000) or the Soil, Vegetation and Snow scheme (SVS: Husain et al. , 2016), (2) the lateral movement of surface (overland) and subsurface (interflow) flows to the drainage system, represented by the WATROF (Soulis et al. , 2000) or PDMROF (Mekonnen et al. , 2014) algorithms, (3) the hydrological routing between river-network grids, represented by the WATROUTE component of the WATFLOOD hydrologic model (Kouwen et al. , 1993b).
To represent the landscape, MESH is based on a model grid that is subdivided into grouped response units (GRU: Kouwen et al. , 1993a) based on land cover, soil type, slope/aspect. Water and energy fluxes are computed at the tile-level (GRUs mapped onto grids) and then aggregated to the grid-scale using weighted averaging based on the areal fractions of GRUs in each grid. MESH runs at a sub-daily time-step forced by seven meteorological variables, namely air temperature, barometric pressure, incoming longwave and shortwave radiation, precipitation, specific humidity, and wind speed. MESH has been widely utilized to simulate land surface-hydrology processes in cold regions (e.g. Razavi et al. , 2010; Haghnegahdar et al. , 2014; Davison et al. , 2016; Yassin et al. , 2017; Elshamyet al. , 2020).
In this study, CLASS is used as the underlying land surface model. CLASS simulates the coupled water and energy balances for a user-defined soil layering that is generalized across the modelled watershed. Above ground, CLASS encompasses four plant functional types, needle-leaf forest, broadleaf forest, grassland and cropland. Below ground, soil parameters (defined by Sand, Clay, and Organic matter percentages) implicitly link soil thermal (i.e . heat capacity and thermal conductivity), and soil hydraulic properties (e.g . porosity and saturated hydraulic conductivity); the latter being defined using Cosbyet al. (1984). During runtime, each soil layer’s temperature and moisture content can evolve and update the associated thermal properties. This occurs down to the depth of bedrock, or the soil permeable depth (SDEP), below which no moisture migration is allowed; only heat can transfer vertically between the soil layers in this region.
CLASS incorporates the Neumann-type boundary condition at the bottom of the soil column (constant geothermal flow) by which the user can replicate the presence of an upward geothermal flux. Flux exchanges with the atmosphere determine the upper boundary condition through the solution of the surface energy balance. Initial conditions include prognostic variables for each soil layer, such as temperature and volumetric moisture content, in addition to other surface state variables.
Fully organic soils can be handled by CLASS using three predefined organic peat types (i.e . fibric, hemic and sapric) based on the work of Letts et al. (2000). Compared to mineral soil, peats have higher porosities (0.93, 0.88,0.83 for the three sub-types respectively compared to 0.49 for clay), higher retention capacity (0.275, 0.62, 0.705), higher residual water content (0.04, 0.15, 0.22), higher heat capacity (2.5 x 106Jm-3K-1), and lower thermal conductivity (0.25 Wm-1K-1) than mineral soils. Thermal properties are taken to be the same for all organic sub-types. Further details are provided in Verseghy (2012).
MESH/CLASS is usually run at a 30-min time-step, and different permafrost characteristics can be output from the simulated temperature profiles. For the present study, the following related aspects of permafrost are considered: temperature envelopes (Tmax and Tmin), mean annual ground temperature profile at the top of the permafrost (MAGTp), active layer thickness (ALT), depth of the zero-annual amplitude point (DZAA), depth to the base of permafrost (BP), thermal offset, surface offset, and date of maximum thaw (ALT-DOY) (see Fig. 1 andTable 1 ).
Possible Position of Fig. 1 .
Possible Position of Table 1 .

Study area and data

The experimental sites selected for this study are near the outlet of the Liard River Basin, Northwest Territories, Canada (Fig. 2 ). The area is located along the divide between sporadic and discontinuous permafrost regions based on the permafrost Map of Canada (Hegginbottomet al ., 1995). More than half of the basin is underlain by sporadic permafrost, mainly in the south, discontinuous permafrost underlays the northern third of the basin, and the rest of the basin is underlain by patchy permafrost (Fig. 2 ). The climate is characterized as subarctic according to the Köppen-Geiger classification (Peel et al. , 2007). The basin is underlain by warm-permafrost (near 0°C), where a rapid reduction in the extent of permafrost in the Canadian sub-arctic has been observed due to climate change (DeBeer et al. , 2016; Connon et al. , 2018). The Liard River basin plays a central role in the sub-continental Makenzie River Basin’s hydrology, as it has the highest runoff coefficient and contributes the largest mean annual flow to the Mackenzie River at the outlet (Woo, 2012).
The availability of soil temperature data is limited to a few experimental sites (e.g.Scotty-Creek (Quinton and Marsh, 1999)), and measurements made during/after the construction of infrastructure for maintenance and monitoring purposes (e.g. Norman Wells-Zama pipeline (Smithet al. , 2004); Yukon Alaska Highway (Oldenborger et al. , 2015)). In this study, two representative permafrost sites were selected due to the availability of soil temperature data at multiple depths and corresponding borehole logs (Fig. 3 ). These were Jean Marie Creek (borehole 85-12B), underlain by sporadic permafrost, and Wrigley Highway (borehole 99TC03), underlain by discontinuous permafrost (Table 2 ), initially installed to monitor the Norman Wells-Zama pipeline’s impact on permafrost. Several Geological Survey of Canada (GSC) reports have been used to extract the thermal and geological data for the current study (Smithet al. , 2004, 2009, 2010, 2016; Ednie et al. , 2012; Chartrand et al. , 2014).
The Jean Marie Creek (JMC) site is dominated by boreal forest (mainly needleleaf) and scattered ericaceous shrubs on peat plateaux where permafrost is warm (Mean Annual Ground Temperature (MAGT) of -0.1°C) and of limited thickness (~ 4m) and the active layer is shallow (~ 1.5m). The data span 1986 to 2000, with no records available in the 21st century. The Wrigley Highway (WH) site is dominated by shrubs with a small black spruce thicket and moss. At this site, permafrost is also warm (MAGT of -0.2 °C) but has a larger thickness (~ 10m) than JMC while the active layer is slightly deeper (~ 2m). Since the two sites are relatively close (~60 km apart), they have similar climatic conditions with an average annual daily air temperature over the 1979-2017 period of -2.5°C and -1.9°C, and average annual precipitation of 430 mm yr-1 and 420 mm yr-1 for the JMC and WH sites respectively.
MESH requires seven climatic variables at a sub-daily resolution to drive CLASS, as mentioned in Section ‎2.1 . Selection of a forcing dataset is constrained by the quality of atmospheric data and the availability of permafrost data, which spans 1986 - 2000 for JMC and 2007-2015 for WH (Table 2 ). A few forcing datasets start prior to the 1980s such as WFD (WATer and global CHange (WATCH) Forcing Data) (Weedon et al. , 2011), Princeton (Sheffield et al. , 2006), and WFDEI (WFD with the ERA-Interim analysis) (Weedon et al. , 2014). However, the WFD and Princeton datasets end in 2001 and 2012 respectively. The combined Global Environmental Model (GEM; Côtéet al. , 1998) atmospheric forecasts and the Canadian Precipitation Analysis (CaPA; Mahfouf et al. , 2007) have been found to compare well with ground observations (Wong et al. , 2017), but GEM-CaPA is not available prior to 2002. Since the WFDEI dataset provides reasonable estimates of climate fields, as shown by Wong et al. (2017) for precipitation, and is available from 1979, covering the duration of the permafrost records, WFDEI is used for driving CLASS for the period 1979-2016, at 3-hour resolution.