The Permafrost and Organic LayEr module for Forest Models (POLE-FM) 1.0

. Climate change and increased fire are eroding the resilience of boreal forests. This is problematic because boreal vegetation and the cold soils underneath store approximately 30% of all terrestrial carbon. Society urgently needs projections of where, when, and why boreal forests are likely to change. Permafrost (i.e., subsurface material that remains frozen for at least two consecutive years) and the thick soil-surface organic layers (SOLs) that insulate permafrost are important controls of boreal forest dynamics and carbon cycling. However, both are rarely included in process-based vegetation models used to 15 simulate future ecosystem trajectories. To address this challenge, we developed a computationally efficient permafrost and SOL module named the Permafrost and Organic LayEr module for Forest Models (POLE-FM) that operates at fine spatial (1 ha) and temporal (daily) resolutions. The module mechanistically simulates daily changes in depth to permafrost, annual SOL accumulation, and their complex effects on boreal forest structure and functions. We coupled the module to an established forest landscape model, iLand, and benchmarked the model in interior Alaska at spatial scales of stands (1 ha) to landscapes 20 (61,000 ha) and over temporal scales of days to centuries. The coupled model could generated intra-and inter-annual patterns of snow accumulation and active layer depth (portion of soil column that thaws throughout the year) generally consistent with independent observations in 17 instrumented forest stands. The model was also skilled at representeding the distribution of near-surface permafrost presence in a topographically complex landscape. We simulated 394.36% of forested area in the landscape as underlain by permafrost; compared a close match to the estimated 33.4% from the benchmarking product. We 25 further determined that the model could accurately simulate moss biomass, SOL accumulation, fire activity, tree-species composition, and stand structure at the landscape scale. Modular and flexible representations of key biophysical processes that underpin 21 st -century ecological change are an essential next step in vegetation simulation to reduce uncertainty in future projections and to support innovative environmental decision making. We show that coupling a new permafrost and SOL module to an existing forest landscape model increases the model’s utility for projecting forest futures at high latitudes.


Introduction
The boreal forest is warming at a rate at least twice the global average (IPCC, 2021;Chylek et al., 2022), which can reduce fuel moisture and causeing climate-sensitive disturbances, like forest fire, to increase Walker et al., 35 2020). Together, pronounced warming and larger, more severe fires are initiating abrupt changes in forest cover, structure, functions, and tree-species composition (Johnstone et al., 2010a;Alexander and Mack, 2016;Walker et al., 2019;Mack et al., 2021;Baltzer et al., 2021); trends that will likely continue for at least the next several decades (Mekonnen et al., 2019;Foster et al., 2019Foster et al., , 2022. This is important because biophysical properties of the boreal forest underpin feedbacks to regional climate (Foley et al., 1994;Chapin et al., 2008;Rogers et al., 2013;Potter et al., 2020), and ~30% percent of all terrestrial organic 40 carbon stocks are stored in the biome (Lorenz and Lal, 2010;Schurr et al., 2018). Some portion of those stocks could be released to the atmosphere and further accelerate warming (Anderegg et al., 2022). Thus, society urgently needs projections of where, when, and why the boreal forest will change.
Ecological legacies are the organismal adaptations (i.e., information), physical materials, and energy that persist in ecosystems through multiple disturbances (Ogle et al., 2015). Legacies will underpin how the boreal forest responds to climate 45 change and fire (Turetsky et al., 2016;Johnstone et al., 2016). For example, adaptive traits, like cone serotiny (cones that stay closed for many years until heated by fire) and asexual resprouting, are information legacies that facilitate postfire forest recovery (Johnstone et al., 2009(Johnstone et al., , 2010a. Thick moss-dominated soil-surface organic layers (SOL) form over decades of postfire forest development, and a portion often escapes burning in the subsequent fire, leading to accumulation of SOL over multiple fire cycles (Walker et al., 2018). This serves as a physical legacy that preserves permafrost (subsurface material that 50 remains frozen for at least two consecutive years) Jorgenson et al., 2010) and shapes tree species composition by controlling seedling germination and establishment . In conjunction with insulative physical legacies, energy legacies of past temperature regimes also maintain permafrost underneath forests where current air temperature would otherwise not support it .
Physical and energy legacies underpin spatio-temporal patterns of permafrost at multiple scales. At the biome scaleIn 55 the boreal forest of North America, permafrost is continuous in the north, becomes discontinuous, sporadic, and is then eventually absent from the southern boreal forest in the south (Obu et al., 2019). Within the discontinuous zone, the permafrost distribution is heterogeneous, varying on fine spatial scales with topography, dominant forest type, and fire history (Brown et al., 2016;Gibson et al., 2018). Permafrost dynamics are particularly important for shaping boreal forest structure and function as well as hydrology (Turetsky et al., 2010;Baltzer et al., 2014;Dearborn and Baltzer, 2021). Within permafrost-affected soils, 60 a portion of the soil column termed the "active layer" undergoes an annual cycle of freezing and thawing. The annual maximum active -layer depth can vary from a few centimeters to several meters (Smith et al., 2022). This freezing and thawing determines the seasonality, vertical distribution, and amount of plant-available soil water and influences nutrient availability (Abbott and Jones, 2015;Young-Robertson et al., 2017).
In response to continued warming, annual maximum active-layer depth is predicted to increase, and the distribution 65 of permafrost will likely contract, with large hydrologic and biogeochemical consequences (Pastick et al., 2015;Schuur and Mack, 2018). Increasing wildfire (Veraverbeke et al. 2017, Phillips et al. 2022) will also impact permafrost by combusting SOLs and altering treeforest regeneration pathways (Baltzer et al. 2021, Johnstone et al. 2010). However, permafrost and the legacies that affect its dynamics, are rarely considered in forest models. In fact, just a handful of models explicitly simulate permafrost (Foster et al., 2019;Gustafson et al., 2020;Kruse et al., 2022), and those that do often operate at relatively coarse 70 spatial (≥ 25 ha grid cells) and/or temporal (≥ monthly) resolutions (but see Kruse et al. 2022, which describes a permafrost module that runs with a five minute temporal resolution). This makes it difficult to capture the fine-scale spatial heterogeneity of permafrost distributions and the effects of daily temperature variability on plant water availability during short, but critical shoulder seasons. Further, most existing permafrost algorithms rely on computationally intensive numerical methods (Sitch et al., 2003;Beer et al., 2007;Karra et al., 2014;Perreault et al., 2021;Yokohata et al., 2020;Westermann et al., 2016), limiting 75 the spatio-temporal resolutions at which they can be applied, particularly across broad domains.
To address this challenge, we present the Permafrost and Organic LayEr module for Forest Models (POLE-FM) that was designed to mechanistically simulate daily changes in active layer depth, annual SOL accumulation, and the associated ecological effects on boreal forests and fire at a fine spatial resolution (i.e., grain of ~1-ha) in a computationally efficient manner (Fig. 1). When paired with a state-of-the-art forest model, such as iLand, the module allows for simulation of complex 80 feedbacks among forests, fire, and permafrost dynamics in topographically complex landscapes under historical and future conditions. In this paper, we describe the module and benchmark its ability to represent permafrost and SOLs in forest stands to landscapes of interior Alaska across days to centuries.

Permafrost and SOL module 85
The module represents daily changes in active layer depth and long-term trends (years to decades) in permafrost presence. Permafrost is represented based on physical principles of heat transport through vegetation and soil media with varying thermal resistances affected by soil moisture content. We incorporate the insulating effects of snow and deep SOLs and capture transient shifts between permafrost regimes (e.g., a transition from temporally continuous to sporadic permafrost due to climate change). Moreover, we aimed for a computationally efficient approach that operates well within the runtime 90 and memory constraints of forest models. The module tracks the energy fluxes that thaw and freeze water at the edge of the active layer (zero isoline, or the depth at which soil temperature is 0 °C), requires only a few state variables, and provides daily values of active layer depth with little computational overhead by avoiding iterative numerical approximations of differential equations.
To capture daily changes in active layer depth, we first estimate the thermal resistances [m² W -1 K -1 ] of snow (when 95 present), SOL, and the mineral soil layer (Eq. 1).
Snow depth is represented as a function of the precipitation that falls during days with mean air temperature below 0°C and the density of snow pack (set at 190 kg m -3 ) (Bonan, 1991;Bennett et al., 2019). We set snow thermal conductivity, 100 . , at 0.3 w m -1 K -1 (Cook et al., 2008). SOL depth is estimated based on the mass of live and dead mosses and litter pools in each grid cell. SOL thermal conductivity, . is set at 0.09 w m -1 K -1 (Hinzman et al., 1991;O'Donnell et al., 2009).
Characteristics of the mineral soil layer that determine its conductivity are explicitly considered. We allow mineral soil thermal conductivity, .
. , to vary with soil texture and soil moisture. We derive mineral soil conductivity following 105 the approach of Farouki (1981) as described in Bonan (2019) (Eq. 2).
. is determined by linearly ramping between saturated conductivity, .
. . , and dry conductivity, . . . , based on a factor, , that varies with relative soil moisture and soil texture, represented separately for unfrozen (Eq. 3) and frozen (Eq. 4) soils. 110 Where is the Kersten number, and is the volumetric soil water content (VWC) relative to the volumetric soil water content at saturation (VWC.sat). .
. . = 0.135 +64.7 is estimated as a function of the conductivity of solids, water, and ice in the matrix, modeled separately for unfrozen (Eq. 6) and frozen (Eq. 7) soils. Using the total thermal resistance from Eq. 1, we can then estimate the daily sum of energy flow ( ; MJ day -1 ) that reaches the zero isoline from the atmosphere above (Eq. 9). Where . is the daily mean air temperature , . . = 0°C, and the constant converts from J s -1 to MJ day -1 . is then used to calculate the daily sum of water that thaws or freezes at the zero isoline based on the enthalpy (or latent heat) of fusion ( ℎ ; 0.33 MJ /liter -1 water). Eq. 9 is also used to estimate the daily energy flux from soil below the active layer by replacing .
with the temperature of the soil below. Deep soil temperatures, here set at 5m, is 130 assumed to be at equilibrium with mean annual air temperature of the previous decade (Riseborough, 2004).
We then model the daily amount of water that thaws or freezes, . . (Eq. 10).
Where . . is constrained to values between -10 and +10 mm (only 10mm of water is allowed to freeze or thaw each day in order to avoid numerical instabilities close to the soil surface). Finally, the corresponding depth of soil that freezes 135 or thaws each day in m, . . , is calculated (Eq. 11).
Since frozen soil (and the water captured therein) is not accessible for plants, the actual water holding capacity of the soil is dynamically modified each day. If soil thaws in a given day, that freshly melted water is added to the soil water pool and the capacity for soil to hold water increases. The approach described here also works for estimating seasonal thawing and 140 freezing of soils in areas not underlain by permafrost.
The SOL component was adapted from Bonan and Korzuhin (1989) and Foster et al. (2019) and represents SOL depth as a function of annual moss net primary production, biomass accumulation, respiration, and turnover. It adds live and dead moss to the fuels for forest fires, and the depth of the SOL influences post-fire tree regeneration. Annual moss productivity is simulated as a function of environmental scalars that represent effects of light attenuation through the forest canopy and moss 145 layer and growth inhibition from fresh deciduous litter. The amount of light that reaches moss for photosynthesis attenuates with increasing forest canopy cover and with increasing moss biomass. Effects of light attenuation are represented by first calculating the amount of light available for photosynthesis in year t as ℎ .
is the leaf area index (m 2 leaf area 150 m -2 ground) of tree cover in year t.
Where is the light saturation point, or the amount of light, relative to the light level above the canopy, above which, an increase in light does not increase moss GPP; set at 0.05. is the light compensation point, or the amount of light, relative to light level above the forest canopy, beyond which moss begins to photosynthesize; set at 0.01.
Field experiments show that fresh leaf litter from deciduous broadleaf tree species strongly inhibits moss productivity (Jean et al., 2020). Such inhibitory effects, , are modeled as (Eq. 14). 160 When . > 0 or 1 when . = 0. Where . −1 is the fresh (previous year's) forest floor deciduous litter biomass in Mg ha -1 . , annual assimilation by moss in year t (kg biomass m -2 leaf area) is then computed (Eq. 15).
Where , the maximum moss productivity per unit leaf area, is 0.3 kg m -2 year -1 (Foster et al., 2019). We estimate effective 165 assimilation in year t, .
Moss productivity in year t, , in kg m -2 biomass then depends on turnover, , and respiration, , in year t (Eq.
Where . −1 is the previous year's moss biomass in kg m -2 and and are empirical parameters set at 0.136 and 0.12, respectively (Foster et al., 2019). The moss biomass pool is updated (Eq. 20).
Note that the biomass pool can shrink if Pt becomes negative, e.g., due to a closing canopy.
Thickness of the live moss layer is calculated as biomass divided by a bulk density of 31 kg m -3 calculated from field observations described in Walker et al. (2020). Dead moss and forest floor litter layer thickness is calculated as biomass divided by bulk density, set at 91 kg m -3 (Walker et al., 2020).
The permafrost and SOL module is implemented in C++ for computational efficiency and is relatively compact (< 180 1,000 lines of code). It is compatible with PC, Linux, or Mac and full full source code and documentation is available under a GNU General Public License (GNU GPL www.gnu.org/licenses/gpl-3.0.html) (See code availability section).While the design is modular, we note that the complex feedbacks between vegetation, permafrost dynamics, and SOL accumulation may require some adaptations and code modifications when integrating our work in different forest models. Below, we detail the integration into the individual-based forest landscape and disturbance model iLand (Seidl et al., 2012a). 185

Coupling the permafrost and SOL with iLand
iLand simulates the growth and mortality of individual trees in spatially explicit landscapes as a function of canopy light interception, climate, nutrient availability, and disturbance (Seidl et al., 2012a, b). The model was originally designed to study effects of natural disturbances, like forest fire, on forest landscapes in the context of climate change (Seidl et al., 2012a).
Thus, iLand emphasizes representation of disturbances and the processes that underpin forest responses to disturbance, 190 including tree-seed production and dispersal, abiotic filters of tree-seedling establishment, and multiple pathways of tree mortality (Seidl et al., 2012a, b;Hansen et al., 2018Hansen et al., , 2020. For an exhaustive technical description of iLand, including carbon cycling and simulation of forest fire, see Appendix A and https://iland-model.org/, which includes full model source code. The proportion of moss biomass that turns over (dies) each year in the new module is fed into the litter layer of iLand's decomposition module. iLand simulates decomposition as a function of climate and pool-specific carbon to nitrogen ratios 195 (Seidl et al., 2012b). The C:N ratio of moss litter is set at 30 (Melvin et al., 2015). Together, live moss, dead moss and forest floor litter layers comprise the SOL in iLand. Wildfire ignition, spread, and severity are partially contingent on downed fuel availability in iLand (Seidl et al., 2014a), and we now include live and dead moss as available fuel in the fire module. When a grid cell burns, the combusted forest floor litter, dead moss, and live moss pools are subtracted from SOL depth.
The tree species that establish in years following fire shape multi-decadal successional trajectories (Seidl and Turner, 200 2022). The depth of burning in the SOL is an important determinant of seedling establishment success because the SOL is often dry and seedlings must expand their roots into mineral soil to access water (Johnstone and Chapin, 2006;Brown and Johnstone, 2012). We therefore included the effect of deep SOL as an additional limiting factor when calculating tree-seedling establishment in iLand. For each 1-ha iLand cell, the probability of establishment is scaled with a negative exponential function following Trugman et al. (2016) (Eq. 21). 205 Where . is a multiplicative factor reducing the abiotic establishment probability in year t, .
ℎ is the depth of the SOL (cm) in year t, and c is a species-specific shape parameter, set at 0.50 for trembling aspen (Populus tremuloides

Model benchmarking
We used a pattern-oriented modeling framework (Grimm et al., 2005) to evaluate the new module by simulating forests of interior Alaska at stand and landscape scales over days to centuries. Pattern-oriented modeling is an approach to benchmarking where patterns of many variables operating at multiple temporal and spatial scales are compared to observational datasets. We chose interior Alaska because it is located in the discontinuous permafrost zone where permafrost presence, moss 215 production, and SOL accumulation vary with dominant forest type, disturbance history, and topography. For example, areas dominated by mature black spruce in lowland valley bottoms and north facing slopes are generally underlain by permafrost and support a relatively productive forest-floor moss layer and thick SOLs. Upland and south-facing slopes are dominated by deciduous trembling aspen and Alaskan birch, which are often not underlain by permafrost, and moss is far less prevalent.
White spruce also inhabits upland positions, on its own, or mixed with black spruce, and contains SOLs of intermediate 220 thickness (Van Cleve and Viereck, 1981). The multiple interacting biotic and abiotic drivers of permafrost and moss productivity create complex landscape mosaics (Johnstone et al., 2010a) that we wanted to ensure the module could produce.
We first evaluated whether the module could generate reasonably realistic daily patterns of snow accumulation/melting and active layer thawing/freezing at the stand level. We then simulated a ~61,000 ha forested landscape to test whether the approach could realistically generate complex mosaics of near-surface permafrost presence, moss 225 productivity, and SOL accumulation consistent with observations.. To ensure robust simulations, we updated an existing iLand tree-species parameter set for interior Alaska (Hansen et al., 2021) (Table B1) and parameterized the iLand carbon cycle (Table   B2) using values derived from the literature.

Temporal patterns of snow and active layer depth
To evaluate whether the module could generate realistic intra-and inter-annual patterns of snow accumulation and 230 active layer depth, we selected 1717 forested sites in interior Alaska that span approximately 700 km. the southernmost site sits along the Alaskan highway at the border between Canada and Alaska. The northernmost site is just south of the Brooks mountain range along the Dalton highway. that wereEach site was instrumented with temperature probes to measure daily soil temperature at depths of zero to six m between 2014 and 2018 (https://permafrost.gi.alaska.edu/sites_list). Seven of the sites were recorded as having an annual maximum active layer depth of less than 2 m (permafrost present). Ten of the sites had an 235 annual maximum active layer deeper than 2 m (permafrost absent). We used the 2 m depth cutoff because it is the maximum effective soil depth assumed in iLand. The sites were initialized from field inventories covering the same domain selected to match the species composition recorded in the soil temperature database (Walker and Johnstone, 2014;. Soil information used to initialize iLand were was extracted from the global SoilGrids250m V. 1.0 (for effective soil depth) and 2.0 (for % sand, silt, and clay) (Hengl et al., 2017). Relative soil fertility, expressed as plant available nitrogen, was set to 240 45 kg ha -1 yr -1 (Hansen et al., 2021). Depth of the SOL was not recorded in the soil temperature database for the 17 sites. Thus, we used photos from the instrumented sites and information on the dominant forest type to assign initial SOL depths to the iLand stands. Sites where researchers recorded dominance of deciduous trees, or where SOLs appeared absent or shallow in photographs were assigned a depth of 0 or 0.07 m to match independent field estimates of SOL depths in deciduous forests located in the Tanana Valley near Fairbanks (Melvin et al., 2015). Sites dominated by black spruce, or where photographs 245 suggested a deep SOL, were assigned a depth of 0.25 m based on field surveys of black spruce stands (Johnstone et al., 2010a).
Stands dominated by white spruce were assigned an intermediate depth of 0.16 m.
Stands were simulated in iLand with 2001-2018 daily climate (minimum and maximum daily temperature, precipitation, shortwave solar radiation, and vapor pressure deficit) from the 1-km Daymet product (Thornton et al., 2021).
We benchmarked simulated maximum annual snow depth and timing of snow melt for the period 2001-2017 (the period when 250 snow observations were available) using a gridded snow product (Yi et al., 2020). This product was developed by integrating downscaled reanalysis data with satellite imagery to provide a continuous estimate of snow depth at 1-km spatial grain. When compared with a meteorological station network (SNOTEL), the gridded observational product had a RMSE of 0.32 m with a bias of -0.09 m in mid-elevations (400-800m) where 70% many of our forested sites were located, and a bias of 0.01m at low elevations (< 400m) where the rest of our sites were located (Yi et al., 2020). 255 We compared simulated and observed maximum annualdaily changes in active layer depth for 2014-2018, the period where soil temperature observations were available, at the seven permafrost sites and maximum annual freezing depth for the 10 non-permafrost sites. with root mean squared error (RMSE). We converted observed daily soil temperatures at depths of 0.03, 0.5, 1, 1.5, 2, 4, and 6 m to active layer depth by identifying the zero isoline with linear interpolation. We also compared the day of year when maximum active layer depth and freezing depth were reached in simulations and observations. 260

Landscape heterogeneity in near-surface permafrost presence, moss productivity, and SOL accumulation
We evaluated whether the module, coupled with iLand, could simulate landscape-scale mosaics of near-surface permafrost (≤ 1 m deep), moss production, and SOL accumulation in a large forested area (~61,000 ha of land area). We initialized the model with a tree-species composition map based on a remotely sensed plant functional type (PFT) product for Alaska and western Canada that classified vegetation as spruce, deciduous, mixed forest, or non-forest  and 265 reflected fire history. We further decomposed PFTs into black spruce, white spruce, trembling aspen, Alaskan birch, mixed forest, potential forest (i.e., areas currently unforested that could support forest in the future), and nonforest using rules based on aspect, elevation, and a permafrost map (Table B3). While this approach allowed us to disaggregate PFTs to the species level, we lack robust datasets to evaluate the accuracy of the species composition map. This is a challenge as dominant tree species determines SOL accumulation and permafrost distribution. In the future, well validated remotely sensed tree-species 270 composition maps would markedly reduce initial condition uncertainty of forest simulations in interior Alaska (Hermosilla et al., 2022).
Initial stand densities, tree sizes, and forest-floor carbon pools (litter, coarse wood, live and dead moss; Table B4) for the appropriate tree species were initialized in the model as early postfire (11 years old) forest based on field inventories described earlier (Walker and Johnstone, 2014;. Because the forest landscape was initialized as entirely 275 early postfire, it did not reflect variation in forest stand age. Thus, we ran a 200-year spin up as a function of historical climate (climate years 1950-2005 recycled randomly with replacement) and simulated fire dynamically to generate spatial heterogeneity consistent with internal model logic, following protocols established in previous iLand studies Turner et al., 2022). We then simulated forests for another 100 years and used this period in all analyses.
We want to eventually conduct simulations with future 21 st century climate. Thus, we used daily meteorological data 280 from the historical period of the CMIP5 generation CCSM4 General Circulation Model (GCM) (Gent et al., 2011) to force landscape-level simulations instead of DAYMET (as was used in the stand-level experiment). This GCM corresponds closely with observed historical climate in Alaska (Walsh et al., 2018), and we statistically downscaled it to a 1-km spatial resolution using quantile matching with Daymet as the observational grid (Hansen et al., 2021).We extracted soils data from the same sources as the stand-level experiment that geographically corresponded to the 1-ha grid-cells in our simulated landscape. 285 Because fire is stochastic in iLand, and an important determinant of permafrost dynamics, SOL depth, tree-species composition, and stand structure, we ran ten replicates and analyzed output from the run with the smallest difference between modeled and observed mean annual burned patch size and annual probability of a fire event.
We compared fire from simulation years 201-300 to observations in the Alaska Large Fire Database from the period 1980 -2021. This database contains perimeters for larger fires (size threshold for inclusion has varied over time, ranging from 290 10-1,000 ha) and point locations for smaller fires in Alaska. We chose years 201 -300 for evaluation because a century aligns with the historical mean fire return interval in Alaska (Johnstone et al., 2010b). We combined these datasets to ensure comprehensive coverage and assumed a circular shape for the smaller fires when perimeters were unavailable. Fire is a stochastic process in iLand, so, we did not expect perfect correspondence between modeled and observed individual fire sizes and locations. Instead, we aimed for the model to generate fire characteristics (i.e., frequency, patch size, annual area burned, 295 and severity) that were generally consistent with the observational record. We took two approaches for benchmarking. First, we compared simulated and observed annual probability of fire occurrence and mean annual burned patch size, as well as the proportion of stems and basal area killed by fire. Second, we compared simulated and observed fire characteristics from the landscape with observed fire characteristics in all of forests of interior Alaska broken into 625 ~ 61,000ha landscapes. This allowed us to determine how the dynamic fire module in iLand performed for our landscape, specifically, and how the model 300 performed relative to the spatial variation in fire regimes across interior Alaska.
We compared the proportion of the landscape underlain by near-surface permafrost in the last 40 years of simulation (years 261-300) to a remotely sensed product of near-surface permafrost presence (Pastick et al., 2015). Fourty years was chosen because we wanted to evaluate permafrost over a multi decadal period and because it aligned with the period used to evaluate postfire SOL combustion and tree seedling density (see below). This product was created by integrating satellite 305 records and other geospatial datasets to predict the probability of near-surface permafrost presence at a 30m spatial resolution with machine learning. Because iLand operates at 1-ha spatial resolution for permafrost, we aggregated the remotely sensed data from 30-m to 1-ha grid cells by calculating the mean probability of near-surface permafrost presence in each 1-ha grid cell. We then used a ≥ 50% probability of permafrost presence, the same cutoff used in the original analysis (Pastick et al., 2015), to map the permafrost distribution. In iLand, near-surface permafrost was considered present in any grid cell where the 310 annual maximum active layer depth was ≤ 1 m in 15 (3850%) of years in the last 430 years of simulation. This cutoff ensured we only included areas that were underlain by consistently frozen ground most years. We compared the total proportion of the landscape underlain by near-surface permafrost and how permafrost presence varied as a function of aspect in simulations and the benchmarking product. We also evaluated how permafrost presence varied as a function of simulated dominant tree species, but did not compare to the benchmarking product because we lack tree species composition maps in interior Alaska. 315 We compared SOL carbon in simulation year 300 separated by forest type to field inventories (Alexander and Mack, 2016;Walker et al., 2020). While benchmarking data was unavailable, we also evaluated landscape variability in total SOL and live moss depth. We assessed SOL combustion by fire in different forest types for model years 26610-300, and as compared model output to the two extensive sets of postfire field plots (Walker and Johnstone, 2014;Walker et al., 2020) also used for initialization. The period of analysis was selected to ensure a sufficient number of fires while balancing 320 the computational intensity of these calculations.
Because near-surface permafrost presence and moss productivity are affected by and feedback to influence forest dynamics, we determined whether the model could realistically represent landscape-level patterns of tree-species composition and stand structure. We explored how landscape patterns of dominant forest type shifted through 300 years of simulation and compared simulated stand density and basal area of each forest type from the end of the simulation with two field inventories. 325 The first was a regional network of permanent plots in interior Alaska collected by the Bonanza Creek Long Term Ecological Research Network site (Ruess et al., 2021). The second inventory was the Cooperative Alaska Forest Inventory, which is a set of permanent plots covering interior Alaska, south-central Alaska, and the Kenai Peninsula (Malone et al., 2009). We reran the 300-year simulation with the SOL and permafrost module turned off to evaluate how the module shaped landscape distributions of tree species composition. 330 We also compared simulated aboveground live tree biomass from the end of the simulation with remotely-sensed estimates of aboveground live woody biomass for interior Alaska and western Canada the same landscape (Wang et al., 2021).
This dataset is a 30-m product that characterizes annual live woody biomass for the years 1984-2014. We aggregated 2014 biomass estimates to the 1-ha spatial resolution of iLand using bilinear interpolation. We further benchmarked snag and coarse wood carbon pools in model year 300 with published field observations Mack 2015, Melvin et al. 2015). 335 To quantify the underpinning drivers of landscape variability in tree-species composition and aboveground live and dead biomass, we compared simulated variation in postfire tree-seedling density by species and SOL depth from years 261-300, with field observations (Walker and Johnstone, 2014; using the same fires that were analysed for postfire SOL combustion. from years 260-300 by species and SOL depth with field observations (Walker and Johnstone, 2014;. Finally, we analyzed the computational efficiency of the module by simulating the landscape with and 340 without the permafrost module turned on to quantify its memory requirement and run time. Dominant forest type was determined using species importance values (IV), a measure of stand dominance based on the relative proportions of species density and basal area. It ranges from zero to two . We considered stands dominated by a particular species if their IV was greater than one. Stands were considered mixed-spruce or mixeddeciduous forest if black sprue and white spruce or aspen and birch IVs summed to greater than one, respectively. Averages 345 in the text are presented as medians and inter-quartile ranges (IQRs) (25 th -75 th percentiles). When comparing simulated and observed datasets, parametric statistics were not used because sample sizes can be increased with simulations to artificially inflate statistical significance. Benchmarking analyses were conducted in R statistical software V. 4.0.4 (R Core Team, 2021) using the packages tidyverse (Wickham et al., 2019) and terra (Hijmans, 2021). When forced with 2014-2018 climate, simulated daily patterns of active layer depth (Fig. 3) and, simulated median annual maximum annual active -layer depth was 1.6 (1.3 -1.8) m, and observed median annual maximum active layer depth was 1.4 (1.0 -1.5) m s closely matched observations fromin seven forest stands underlain by permafrost (RMSE of 0.37 m) 360 ( Fig. 2C). Simulated daily patterns of active layer depth also corresponded well with observations (Fig. 3). On average, maximum annual active -layer depth occurred 20 days later 15 days later in iLand than in observations with a IQR of 10 days earlier to 39 days later. The model also reasonably recreated maximum annual freezing depths at 10 forest stands not underlain by permafrost Simulated and observed median annual (maximum freezing depths were 2.0 (1.9-2.0) m and 1.9 (1.9-2.0) m, respectively RMSE of 0.44 m) (Fig. 2D). On average, the maximum annual freeze depth was reached 10 1.3 days earlier in 365 simulations than in observations with an IQR of 28 days earlier to 7 days later than observations.

Landscape-level fire characteristics
Mean annual burned patch size was 3,628 ha and annual probability of a fire event was 11% in the best of the ten replicate landscape simulations and differed from observed values by only 5% and 8%, respectively (Fig. 4A, 4B). However, among all ten replicates, burned patch size and probability of a fire event differed by as much as 44% and 42%, highlighting 370 the stochastic nature of fire. Both observed and simulated fire metrics for the landscape were also representative of observed fire characteristics in all 625 sampled 61,000 ha landscapes across the boreal domain of Alaska (Fig. 4C). On average, 73 (70 -76) % of stems and 52 (46 -60) % of basal area was killed by fire in the model (Fig. B1).

Landscape near-surface permafrost, moss, and SOL depth
The model simulated 39.3% of forested area in the landscape as underlain by permafrost between years 261-300; 375 compared to the estimated of 33.4% of forested area from the benchmarking product (Fig. 5A). Aspect was an important determinant of permafrost presence in the model and in observations (Fig. 5B). Simulated permafrost was over represented on north-facing slopes, as compared to the benchmarking product, but corresponded well on all other aspects. Near-surface permafrost presence also varied with dominant tree species in iLand. Seventy-onethree percent of simulated black spruce forest area was underlain by near surface permafrost, followed by 51% of white spruce forest, 14% of aspen dominated stands , 11% 380 of mixed spruce, 9% of white spruce forest, 2% of aspen dominated stands, 16% of mixed deciduous forest, and 0.21% of birch dominated forest.

Landscape-level tree species composition and forest structure
Between simulation year 0 and 300, forest cover increased from 48,811 ha to 60,629 ha, as trees colonized areas initialized as potential forest. The model was initialized with black spruce forest comprising 41% of the land area, followed 395 by white spruce (22%), aspen (7%), birch (6%), and mixed forest (5%) (Fig. 7A). By year 300, the land area dominated by black spruce remained high at 40% (Fig 7B). However, white spruce-dominated forest area declined markedly to 2% because black spruce trees colonized white-spruce stands, as is commonly found in interior Alaska (Van Cleve and Viereck, 1981;Burns and Honkala, 1990). At the end of the simulation, mixed spruce stands comprised 42 % of land area. Aspen and birch also intermixed by year 300, with mixed-deciduous forest covering 11% of the landscape. The land area dominated by aspen 400 in year 300 declined to 2%, and birch-dominated forest declined to 3% of the landscape. Rerunning simulations with the permafrost and SOL module turned off led to markedly different tree species composition compared to initial conditions and after 300 years of simulation where permafrost and SOL were dynamically represented (Fig. 8).
Stand density and basal area in the model corresponded well with multiple field observation datasets in year 300 (Fig.   98). Aspen and birch stands were most dense, followed by black spruce and white spruce dominated stands. Deciduous 405 dominated stands also had the greatest basal area, followed by white spruce and black spruce stands (Fig. 98B).

Computation efficiency
The memory footprint of the permafrost and SOL module was approximately 15 MB (~ 0.1% of total memory footprint), and it increased overall run time by 1%.

DiscussionConclusions
Ecological legacies will determine how forests are affected by climate change and increasingly prevalent 420 disturbances, like fire (Turetsky et al., 2016;Kannenberg et al., 2020;Hansen et al., 2022). However, some legacies uniquely important to the structure and functioning of boreal forests (e.g., permafrost and SOLs) are rarely considered in models used to project 21 st -century ecological change. Here, we present a new permafrost and SOL module that operates at fine temporal (daily) and spatial (1-ha) scales and is computationally efficient. The module simulates daily changes in active layer depth, moss production, and annual SOL accumulation (Fig. 1). When coupled to a forest model, it also represents the complex 425 ecological effects of permafrost and SOLs on boreal forests and fire. With some exceptions discussed below, bBenchmarking results demonstrate the model recreates temporal and spatial patterns consistent with observations at stand to landscape scales over days to centuries. Our model will contribute to improving 21 st -century projections of boreal forest change.
Process-based simulation models are powerful tools for assessing how forests will change (Seidl, 2017;Albrich et al., 2020;Fisher and Koven, 2020). Forests often respond slowly to stressors relative to other ecological systems (Hughes et 430 al., 2013;Turner et al., 2022). As a result, models must capture dynamic feedbacks among variables and represent the key legacies that accumulate over decades to centuries in order to project future trajectories of forests . Our objective was to mechanistically represent permafrost and SOLs and capture effects of daily variability in weather as well as the feedbacks that arise among forest dynamics, fires, and permafrost in topographically complex landscapes. The model was skilled at capturing daily patterns of freezing and thawing as well as the inter-annual variability in maximum thaw depth but and the timing of maximum thait generally occurred later in simulations than in observations. There are a number of potential reasons for this. First, the model does not track the moisture content of the SOL separately from the mineral soil layer. In reality, the low bulk density of SOL relative to mineral soils leads to more variable moisture content, and thus, a greater range of thermal conductivities, which could lead to slower thawing in simulations if simulated SOL moisture was lower during the spring thaw (Fisher et al., 2016). Another potential explanation is that forest structure (density and leaf area index) and tree 440 species composition have been shown to strongly modulate microclimate and permafrost thaw in complex ways (Stuenzi et al., 2021). Effects of forests on microclimate are also not yet included in the model. Finally, snow depth and melt plays a critical role in active layer dynamics. OurThe model also reasonably recreated snow accumulation patterns for most years, but overestimated depth in years where snow pack was unusually deep., This is likely because iLand takes a relatively simple approach to simulating snow derived from Running and Coughlan (1988), which is not particularly mechanistic.. For example, 445 we used a single snow-density parameter value, which ignores compaction. In reality, snow density varies tremendously across landscapes and over time. iLand takes a relatively simple approach to simulating snow derived from Running and Coughlan (1988). In the future, our approach would benefit from separately tracking moisture content of the SOL, a representation of forest structure and composition effects on microclimate, and a more advanced snow model could be added that includes key processes affecting snow depth and conductive properties, including the representation of variation in snow density, freeze 450 thaw cycles, and sublimation (Bormann et al., 2013;Jafarov et al., 2014).
At landscape scales, the model generally captured mosaics of near-surface permafrost, moss production, and SOL accumulation in a large forested area (~61,000 ha of land area), but it did overestimate the area underlain by permafrost on north facing slopes. This may have occurred for two reasons. First, the climate data used to force iLand was statistically downscaled to a 1-km resolution from a global general circulation model, and thus, does not perfectly capture variability in 455 climate as a function of fine scale variation in aspect and topography. Further downscaling using lapse rates might help improve simulations. Further, dominant forest type varies strongly with aspect in interior Alaska, and in turn, shapes SOL thickness and permafrost distributions. However, we lack landscape level maps of individual tree species distributions to initialize the model. In the future, well validated remotely sensed tree-species composition maps would markedly reduce initial condition uncertainty of forest simulations in interior Alaska (Hermosilla et al., 2022) and could improve landscape level simulations of 460 permafrost distribution.
The module was designed to represent permafrost and SOL effects on forest dynamics and fire. In particular, it determines the water available to plants and accumulation of forest floor biomass, which serves as fuels for fire and influences postfire tree regeneration. When coupled with iLand, the model reproduced common secondary successional trajectories found in interior Alaska, including self-replacement and disturbance-induced abrupt transitions in forest types (Johnstone et al., 465 2010a(Johnstone et al., 465 , 2016. For example, when thick SOLs remained after fire in black-spruce stands, self-replacement was common, leading to recovery of forests functionally and structurally similar to the prefire stands (Anderson et al., 2003;Johnstone and Kasischke, 2005). In contrast, when fires combusted most of the SOL in black spruce stands, abrupt transitions from spruce-to deciduous-dominated forest (mixtures of aspen and birch) occurred, consistent with regional trends documented in the last few decades (Johnstone et al., 2010a. 470

Conclusions
The boreal forest biome is warming at least two times faster than the global average (IPCC, 2021), causing climate sensitive disturbances, like fire, to increase in frequency and severity Walker et al., 2020). Our permafrost and SOL module will help process-based modelers produce more accurate projections of how forests in the biome are likely to change over the next century. Better projections will resolve a number of important uncertainties, including 1) where 475 increased burning due to climate change may reduce boreal fuel loads such that fire-self limitation emerges (Héon et al., 2014;Buma et al., 2022); 2) when shifts in postfire successional trajectories will initiate biophysical feedbacks that further alter regional climate; and 3) how climate change, fire, and permafrost thaw will interact to reshape boreal carbon cycling (Schurr et al., 2018;Schuur and Mack, 2018;Mack et al., 2021). Because boreal forests have disproportionate impacts on the climate system through biogeochemical and biophysical pathways, such information is essential to inform innovative and effective 480 global climate mitigation and adaptation strategies.

Carbon cycling in iLand 485
iLand dynamically models carbon in live foliage, branch, stem, and root compartments, and in standing snag, forestfloor litter, downed coarse wood, and mineral soil organic material pools (Seidl et al., 2012b). Primary production is simulated with a radiation use efficiency approach. Carbon fixed by trees is then allocated to different tree compartments based on allometric equations, representing functional balance. Influxes of carbon from live compartments to dead organic matter pools are calculated based on leaf turnover rates, tree mortality, and snag dynamics. Snag fall occurs over time based on a species-490 specific half-life. When snags fall, they are added to the downed coarse wood pool. Decomposition of dead organic matter pools is represented with a pool-and species-specific optimal decomposition rate (10°C, no water limitation) that is then modified by prevailing temperature and precipitation.

Forest fire in iLand
The model also includes robust representations of several natural disturbances, including forest fire (Seidl et al., 495 2014a, b;Hansen et al., 2020). Fire occurrence and spread are dynamically simulated at a 20-m resolution as a function of 20 th -century fire probability and size distributions, landscape topography, model-generated wind speed and direction, and the proportion of total downed litter and coarse wood pools that are burnable, which is determined by fuel moisture (as quantified by the Keetch Byram Drought Index; KBDI). For every 20-m grid cell that burns, the available fuels are assumed combusted.
Percent crown kill of live trees is estimated as a function of tree size, available fuel loads and aridity. For the portion of live 500 tree canopies that are killed, we assumed 90% of foliage, 50% of branch, and 30% of the burned stem biomass is combusted.
Tree mortality from fire is simulated probabilistically based on tree size, percent crown kill, and bark thickness; a model parameter that varies by tree species. If a tree dies, the non-combusted foliage and branches are added to the downed litter and coarse wood pools. Portions of killed tree stems that were not combusted enter the standing snag pool. 505 510 Appendix B

Species Rule
Black spruce  PFT is spruce and aspect is north  PFT is spruce, aspect is flat, and permafrost is present  PFT is woodland White spruce  PFT is spruce and aspect is not north  PFT is spruce, aspect is flat, and permafrost is not present Trembling aspen  PFT is deciduous and aspect is south