Increased rainfall stimulates permafrost thaw across a variety of Interior Alaskan boreal ecosystems

Earth’s high latitudes are projected to experience warmer and wetter summers in the future but ramifications for soil thermal processes and permafrost thaw are poorly understood. Here we present 2750 end of summer thaw depths representing a range of vegetation characteristics in Interior Alaska measured over a 5 year period. This included the top and third wettest summers in the 91-year record and three summers with precipitation close to mean historical values. Increased rainfall led to deeper thaw across all sites with an increase of 0.7 ± 0.1 cm of thaw per cm of additional rain. Disturbed and wetland sites were the most vulnerable to rain-induced thaw with ~1 cm of surface thaw per additional 1 cm of rain. Permafrost in tussock tundra, mixed forest, and conifer forest was less sensitive to rain-induced thaw. A simple energy budget model yields seasonal thaw values smaller than the linear regression of our measurements but provides a first-order estimate of the role of rain-driven sensible heat fluxes in high-latitude terrestrial permafrost. This study demonstrates substantial permafrost thaw from the projected increasing summer precipitation across most of the Arctic region.


INTRODUCTION
Arctic amplification has caused enhanced warming in Earth's Polar Regions compared to lower latitudes 1,2 . Commensurate with this are projections for a year-round increase in precipitation [3][4][5] that is expected to be greater than at lower latitude locations 6 . Recent studies have already reported increased summer (wet) precipitation in the Arctic [7][8][9] and this is projected to continue in the future 6,10,11 . The fraction of wet (liquid) to dry (snow) precipitation is also projected to increase as the summer season lengthens 6,12 . Increases in projected Arctic summer precipitation are larger under high emission and warming scenarios than for lower emission scenarios 13 .
It is clear that increased Arctic summer precipitation will influence hydrologic processes 14 and the soil thermal regime 15 of the permafrost that underlies 24% of the northern hemisphere. Permafrost soils contain large carbon stores (twice what is currently in Earth's atmosphere) and permafrost, when kept frozen, provides stable ground surface conditions for infrastructure. Increased rainfall has recently been linked to methane release from thawing soils 16 . Changes in snow conditions and liquid water are known to affect the permafrost thermal state though some aspects of heat transfer (latent heat) are better understood than others (advective heat transport) 17 . Most of the studies linking precipitation to permafrost thaw have focused on an increasing winter snowpack that provides thermal protection against winter cold [18][19][20] . Relationships between the depth of the seasonally thawed "active layer" and increasing rainfall have been documented in a few studies 21,22 , but how these relationships vary over time, space, or a changing precipitation regime in northern ecosystems is poorly understood.
In permafrost terrains, the seasonally thawed active layer depth is controlled by local climate, vegetation and soil properties, hydrology, topography, and the timing and thermal conductivity of the previous winter season snowpack 23 . In terrains with near-surface permafrost, the seasonal thaw front moves downward throughout the summer, so end of summer season surveys represent the maximum active layer extent 24 . In many locations the permafrost is "ecosystem protected" against summer thaw by the local vegetation cover and soil and this provides resiliency against top-down thaw or disturbance 25,26 . Most of the studies linking vegetation to permafrost thermal protection have focused on the thermal conductivity or albedo of the vegetation represented. These studies report tussock vegetation 27 and areas with thick Sphagnum moss or well-developed surface organic layers 28 provide the strongest protection against seasonal thaw and permafrost degradation. Deciduous forests, soils with low organic matter or low ice content, and areas where surface vegetation has been disturbed by fire or human activities provide the least protection against surface thaw 29,30 . A key knowledge gap remains in our ability to project how (and where) these surface properties regulate the sensitivity of permafrost to changes in the type and amount of precipitation.
Here we explore the relationships between anomalously wet summers and thaw depth using 2750 end of summer active layer thickness measurements at a variety of locations in discontinuous permafrost near Fairbanks, Alaska (Fig. 1). Our measurements were taken between 2013 and 2017 at four sites along transects representing a variety of vegetation cover types (Supplementary Figs. 1 and 2). Measurements were repeated in the same exact locations every year to allow us to detect permafrost change and minimize spatial variation in active layer depth attributable to soil microform or vegetation structure heterogeneities. Our active layer measurements represent a diverse set of vegetation, including tussock tundra, black spruce forest with Sphagnum, mixed deciduous forest, wetlands, and disturbed areas. The measurements on trail crossings, clearings, and other disturbed locations allow us to examine the sensitivity of pristine versus disturbed permafrost to changing precipitation patterns. In total, the five terrain types in our study areas represent~70% of the land cover present in the boreal and taiga of the Arctic 31 . We modeled our active layer depth measurements using summer rainfall and ecosystem type as fixed effects and compared these to parameters commonly attributed to seasonal thaw depth in permafrost terrains: snowfall, summer mean temperature, summer heating degree days, the previous year's active layer depth, and the volumetric and gravimetric moisture content of active layer and permafrost soil.
Our measurement period represented two summers with wet (i.e. liquid) precipitation amounts close to the 18.5 cm per summer mean over the 91-year record, one slightly drier summer, and two summers representing the most (2014; 37 cm cumulative precipitation) and third most (2016; 33 cm cumulative precipitation) wet precipitation totals measured in the 91-year meteorological record for the area (  Table 1). The 2014 summer was unique because it included eight of the top one hundred 24 h wet precipitation totals in the 91-year record but the least cumulative precipitation in the period between May to mid-June of the 5 years we analyzed. For example, the first major rainfall in 2014 did not occur until mid-June and half of the precipitation fell during large events between mid-June and early July. Major flooding, surface runoff, and ground inundation at our sites started in mid-June and lasted throughout the summer, particularly in tussock tundra ecotypes. The other wet summer captured by our study (2016) yielded only one of the top one hundred 24 h wet precipitation totals. In contrast to 2014, in 2016 numerous rainfall events started in mid-May and continued more consistently throughout the summer. There was standing water at some field locations but not nearly as much inundation and flooding as was evident in 2014. This is likely because the 2016 summer had fewer major rain events and far more frequent smaller ones. Thus, despite a similar amount of precipitation in the summers of 2014 and 2016, the 2016 summer experienced more consistent flux of wet precipitation to our field sites over a far longer period of time. There was little to no inundation or standing water at the field sites in 2013, 2015, or 2017. Neither 2013 nor 2015 had any 24 h wet precipitation events that were in the top one hundred for the area. 2017 had two top one hundred 24 h wet precipitation values but the total for that summer was only slightly above the mean.

RESULTS
Wet precipitation and seasonal thaw of permafrost Increasing summer rainfall led to enhanced permafrost thaw across all five ecotypes (Figs. 3 and 4, Supplementary Table 2). For example, the maximum active layer depth was closest to the surface in the driest summer (2013) and was considerably deeper in the two anomalously wet years (2014 and 2016). Across the entire measurement period increasing rainfall led to an average of a 0.7 ± 0.1 cm increase in active layer depth per cm of additional precipitation. The wet summers led to a steady increase in active layer depth between 2013 and 2017. Following the extremely wet summer of 2014, thaw depths did not return to 2013 values despite summers of 2015 and 2017 being similar in air temperature and wet precipitation totals. In addition, following the anomalously wet summer of 2016, the 2017 thaw depths were only slightly shallower, despite temperature and wet precipitation values in 2017 being closer to the long term mean.
Disturbed areas and wetlands showed the greatest absolute change in active layer depth, with 0.99 and 0.89 cm of additional active layer thaw per cm of rainfall, respectively (Fig. 5). It is apparent from the linear fit y-intercept values there is a wide variation in the typical mean active layer depth for the different ecotypes. For example, the disturbed ecotype y-intercept value of 65.5 cm is the deepest active layer while the tussock tundra value of 46.2 cm is the shallowest. This is likely because disturbed areas have limited ecosystem protection while tussock tundra has been reported to have the greatest ecosystem protection properties of our ecotypes 27 . Tussock tundra provides wetter soils, which conduct surface heat downward better than dry soils due to an increase in thermal conductivity from wetting of pore spaces in the soils [32][33][34][35] . Disturbed sites commonly exhibit subsidence following permafrost thaw 29 and this promotes surface water ponding and, likely, additional thaw and associated subsidence. Though we made the thaw depth measurements in October, the sites were visited monthly during each of the summers. Both the wetland and disturbed ecotypes exhibited mostly dry surfaces during low precipitation summers but became inundated with standing water during the wet summers. These two ecotypes have limited forest cover (Supplementary Fig. 1) and thus low interception rates and potentially higher evaporation rates. There was also standing water in small localized low spots at the tussock shrub sites. However, the mixed forest and Sphagnum moss spruce forest had no visible standing water during our survey measurements. A possible explanation of why the low lying, low gradient areas exhibited the greatest rain-induced thaw is that the additional surface water that ponded in low lying areas migrated downward to warm the seasonally frozen and permafrost soils. Tussock shrubs have been found to provide the strongest ecosystem protection of the five ecotypes represented by our study areas and this is due to their thick organic layer and generally poor drainage properties 25 . At our study sites, tussock tundra yielded the shallowest active layer depths and the greatest relative increase in thaw depth in response to rainfall (0.75 cm of additional thaw per cm of an increase in rain). While wetlands and disturbed sites experienced a 25% increase in active layer depth between 2014 (the wettest summer) and 2013 (the driest summer), mean tussock tundra active layer depths were 51% deeper in 2014 compared to 2013. Tussocks grow on flat, low gradient surfaces over decadal to century time periods and provide strong ecosystem protection against summer thaw 26 . The tussocks at our study sites were surrounded by standing water during wet summer time periods, particularly in 2014 and 2016 ( Supplementary Fig. 1, bottom left photo). This inundation and commensurate lateral flow of surface water provided long-term access to soil pore waters following rain events and likely increased sensitivity to rain-induced thaw at our tussock sites.
The other terrestrial environments, such as mixed forests and conifer forests with thick moss layers, had the least amount of thaw, yielding 0.56 cm and 0.25 cm of thaw per additional cm of rain, respectively. Standing water was extremely rare at these sites for two possible reasons. First, there was likely less throughfall due to interception by the thick canopy vegetation above the forest floor. Studies quantifying interception in boreal forest sites report interception rates range from 23% in spruce 36 to 30% in larch 37 . Second, surface water was likely absorbed by the thick Sphagnum layer with large water storage capacity and upward water wicking through strong external capillary action 38 . The spruce-moss combination was the most resistant vegetation complex to thaw due to canopy shading as well as the strong thermal protection of Sphagnum moss and well-developed organic against summer thaw 28 . Due to these mechanisms our results highlight that these ecotypes are the most resistant to rain-induced thaw. A grid of 121 active layer depths since 2009 is available as a subset of the Farmer's Loop site measurements. The area only represents moss spruce forest and a linear fit between active layer depth and summer mean temperature, cooling degree days, heating degree days, and total summer wet precipitation yields the greatest correlation of determination (0.53) for the relationship between active layer depth and mean summer wet precipitation.
We explored the effects of snowfall, summer mean temperatures, summer heating degree days, the previous year's active layer depth, and volumetric and gravimetric moisture content (of both active layer and permafrost soil) on seasonal thaw depths. None of these relationships were related to active layer depth within or across ecotypes from 2013 to 2017. For example, 2016 (the third wettest year in the 91-year record) was the coldest summer in our measurement period, yet there were elevated thaw measurements across all land cover types. Four of the summers in our measurement period were warmer than the 91-yr mean summer temperature (11.8°C) but only slightly (Supplementary Table 1). The lack of dependency between summer season thaw and winter precipitation, summer air temperatures, or the ice content of the seasonally thawed soils supports our interpretation that changes in active layer depths across our measurement period were driven predominantly by interannual variation in rainfall.
The projected future wetter and warmer conditions across most of the Arctic will trigger shifts in vegetation cover and disturbance regimes and this study demonstrates how land cover types govern relationships between summer precipitation and permafrost thaw. Human land use (roads, trails, and infrastructure) represents a much smaller proportion of disturbed landscape in the Arctic relative to wildfires, which combust the thermally protective surface organic mat and vegetation 29 and accelerate thermokarst and degradation of ice rich permafrost [39][40][41] . Though our disturbed sites are not fire scars they represent areas where the surface vegetation (including protective moss and organic mats) has been removed and this is similar to what a high severity fire does. Conifer forests tend to be vulnerable to burning because they are largely comprised of plants with fire-adapted traits. Our study shows that permafrost in conifer forests is less sensitive to rainfall than in other ecotypes. Increasing fire severity is triggering vegetation state-changes away from conifer-dominated stands 42 to disturbed or mixed forest ecotypes. Our results suggest increasing areas of disturbed, deciduous, or open (unforested) habitat could promote permafrost thaw under projected future wetter summer conditions. However, projecting how, where, and at what rate the changing high latitude precipitation regime will affect the thermal state of permafrost will require better physics based ecosystem models.
Estimating the thermal inputs associated with rain-induced thaw To better parameterize the relationships between summer precipitation and the thermal state of permafrost we applied a simple energy budget model to calculate the amount of seasonal thaw attributable to the added heat from wet precipitation (see "Results" for more detail). We ask the question whether our measured increases in permafrost thaw are simply due to the thermal inputs of rain. Our first-order model is applicable because the end of summer season active layer depth represents an energy increment of latent heat and precipitation provides an added thermal energy input to the ground that causes thaw.
Our approach to estimating the added heat of precipitation is highly simplified and does not consider the full surface energy budget, for example, changes in incident solar or longwave radiation, which might be expected to covary with precipitation on a seasonal basis. Nonetheless, it provides a first-order estimate of the role of rain-driven sensible heat fluxes to our high-latitude terrestrial ecosystems. Estimated energy inputs by rain of~10 MJ/ m 2 are approximately 20% of the reported magnitude of total summertime ground heat flux measured at the base of the active layer at sites in similar climates 43 . Estimated porosity (θ) values from the sites (Supplementary Table 3) are~0.73-0.90 cm 3 cm −3 . This yields an estimate of the expected additional thaw increment per unit of rainfall of 0.15-0.24 cm thaw/cm rainfall. This is up to 4 times smaller than the observed linear correlation between rainfall and active layer thickness from our measurements (0.25-1.0 cm thaw /cm rainfall; Fig. 5), suggesting much of the variation in our measurements cannot be explained by 1-D thermal inputs of rain alone. However, our calculated estimate of the relationships between precipitation and thaw is approximately the same size response we measured for the conifer forest sites. One potential explanation for the differing magnitude and greater rates of observed thaw per increment of additional rainfall in wetter ecosystems may be due to surface and shallow subsurface water convergence. Lowland areas likely receive more water from higher   in catchments and thus receive a commensurate increase in the heat transported by the water. As noted earlier, our lowland sites were flooded throughout most of the wet summers while the other ecotypes did not have evidence of the consistent and sustained presence of water at the surface. Thus, heat advection associated with lateral flow processes may dominate the vertical advection in such regions. This explanation would be consistent with our finding that the observed rates are best explained by the heat budget approach in upland (and less inundated) ecosystems.

DISCUSSION
Following most disturbances in permafrost terrains, ice (massive ice or interstitial pore ice) in the permafrost melts and the volume change leads to ground surface subsidence. This, in turn, can lead to surface water inundation or shallow subsurface flow in low gradient locations, particularly in lowland landscapes. Our results suggest this process, along with increased summer precipitation, may cause a positive feedback to permafrost thaw and degradation. When permafrost degrades and the land subsides, this provides low lying areas on which water can pond and further thaw the permafrost soils. Though the low lying areas between tussocks are generally covered by the tall grassy shoots on top of the tussocks, ponded open water can lead to increased solar absorption which would, in turn, warm the ponded water. As stated earlier, our simple energy budget model did not adequately account for the potential that this ponding of surface and shallow subsurface water can provide a greater source of heat to thaw permafrost soils.
During the exceptionally wet seasons at our sites, it is likely that thaw encompassed the transition zone, a layer of high ice content soil that typically regulates surface permafrost temperatures because of the high latent heat required for it to thaw 44 . At our sites, the shallowest active layer depths were measured in 2013 and following the extremely wet summer of 2014 thaw depths similar to what was measured in 2013 were never repeated despite two other summers with relatively low summer precipitation similar to 2013. This suggests that at some or all of our sites the protective transition zone was thawed and did not re-form. It is likely the warming associated with elevated wet precipitation during one individual summer lasts for multiple years as the increased soil water content, particularly at the base of the active layer, provides a barrier to subsequent freezing in the next winter. These thermal memory effects may prime the system for continued thaw in the following summer(s).
Recent studies have suggested active layer deepening exposes reduced soil horizons to oxidizing conditions and mobilizes trace metals and nutrients (including the vast stores of permafrost carbon) out of soils and into rivers or streams or the atmosphere 45,46 . As a consequence, the deeper active layers at our sites from rain-induced permafrost thaw have major ramifications for the oxidation and reduction processes that govern the permafrost carbon cycle and the potential for trace metal mobilization from frozen soils. It is also clear abrupt thaw processes increasingly being identified across the Arctic 47 are not well addressed by currently available ecosystem models, particularly with respect to precipitation controls on the permafrost thermal state and potential positive feedbacks between increased summer precipitation and top-down permafrost thaw or thermokarst development.

METHODS
The field measurements for this study were made in the vicinity of Fairbanks, Alaska from 2013 to 2017 (Fig. 5). The area, part of the boreal biome, has a continental climate with a mean annual temperature of −3.3°C and mean total annual precipitation of 28 cm of water equivalent 48 . Mean summer air temperatures have increased 2°C since 1930 and are projected to increase an additional 1-2°C by 2040 49 .
Typically, mean daily air temperatures at the field sites remain above 0°C from early May to early October (160 days). We made repeated thaw depth measurements each fall during the second week of October to represent the maximum depth of seasonal thaw (otherwise known as the "active layer depth"). Measurements were made every 4 meters along four linear transects. These included a 400 m transect "Farmer's Loop 1" at 64.877N, 147.674W and a 500 m transect "Farmer's Loop 2" at 64.874N, 147.677W. These two transects represent a variety of ecotypes including regions of mixed deciduous forest, wetland, tussock tundra, moss-black spruce forest, and some disturbed trail crossings. We also measured active layer depths within an 11 by 11 grid at a 4-meter spacing that was established in 2005 following previously established protocols for the Circumpolar Active Layer Monitoring (CALM) network (https://www2.gwu. edu/~calm). The site, CALM site 41, is located near the Farmer's Loop transects at 64.876N, 147.667W and is located in a black spruce forest with 10-30 cm of Sphagnum moss cover. A 500 m transect was repeatedly surveyed at the Creamer's Field migratory bird Refuge (64.868N, 147.738W) that crosses through mixed deciduous forest, tussock tundra, and some disturbed trail crossings. A 400 m transect was established above the CRREL Permafrost Tunnel in Fox (64.950N, 147.621W). Vegetation at this site is Sphagnum moss-black spruce forest with some 2-10 m wide disturbed clearings and trail crossings.
The 550 active layer depths were made each year at small numbered flags installed at each measurement location along the transects and inside the 121-point grid. The flags allow for the collection of repeat thaw depth measurements at the same exact location to ensure repeat measurements of vertically heterogeneous soil and vegetation features, particularly low lying areas, trail crossings and clearings, high centered polygons and their troughs, or tussocks. The location of each pinflag was surveyed for elevation, northing, and easting with a differential global positioning system (dGPS) at 10 cm accuracy.
We measured seasonal thaw depth using a 1 cm diameter graduated metal rod ("frost probe") that extended to as much as 2.5 m in length. At each measurement location the frost probe was pushed downward into the ground to refusal to establish the distance between the ground surface vegetation and the ice-bonded base of the active layer/top of permafrost 24 .
An 8 cm diameter SIPRE corer was used to collect 2-3 m deep cores for volumetric ice content measurements at locations representing the major ecotypes at each site following established methods 50 . A Geoprobe 7822 Direct Push Technology track mounted drill rig was used to collect additional cores for gravimetric ice content measurements. Volumetric and gravimetric moisture contents from the cores are included in Table S1. In the 113 SIPRE core samples collected at our sites the soil porosities we calculate (using a particle density of 2.66) range from 76 to 90% (Supplementary Table 3). Using these porosities, the amount of wet precipitation for each summer, and the daily mean air temperature for each day that a rainfall event occurred in each summer we calculate values for the amount of additional thaw attributable to wet precipitation (Supplementary Table 4).
Statistical analyses were performed with the fit model application within the software program JMP version 13.0 (SAS Institute, Phoenix, Arizona). We modeled active layer depth using summer rainfall and ecosystem type as fixed effects and transect ID (location of each measurement) as a random effect to account for repeated measures. Results of this mixed effect model showed that the repeated measurements explained little variance in thaw depth and was not significant (Wald p = 0.66). Thus, for the remaining analyses, we utilized a general linear model (GLM) that included the fixed effects only. Results from the GLM show that the precipitation-thaw depth relationships varied by an ecotype interaction (Fig. 5). Within ecotypes, relationships varied in significance from the moss spruce forest (p < 0.001, t = −5.23), tussock tundra (p = 0.007, t = 2.69), wetland (p < 0.001, t = 11.51), disturbed (p = 0.09, t = −1.70), and mixed forest (p = 65, t = 0.45) ecotypes.
An estimate for the maximum total extra heat, H rain , delivered to the base of the active layer, per unit of rainfall, over a given patch of ground is (in kJ/cm rainfall): Where Cp is the heat capacity of water at constant pressure (4.1855 kJ kg −1°C−1 ), ΔT is the difference in temperature (°C) between the rainfall and the soil at the base of the active layer, ρ is the density of water (kg/m 3 ), and a is the area considered (m 2 ). We assume soil at the top of the seasonal thaw layer at any given depth or time is isothermal (0°C) and thus thawing the soil requires overcoming the latent heat of fusion and that the temperature of precipitation is equivalent to the mean air temperature T.A. Douglas et al.
during individual precipitation events. This ranged from 10.6 to 13.7°C. There was no statistically significant relationship between air temperature and amount of precipitation within a summer season (Supplementary  Table 1). We used the mean air temperature during all precipitation events for a given year for each summer of the five years we studied. Since the base of the active layer is, by definition at the freezing point, we assume it is 0°C.
The extra heat required to deepen the active layer below a given patch of ground, H ALD , by a given amount is (in kJ/cm thaw): Where θ is the soil porosity (cm 3 cm −3 ), Q is the latent heat of fusion of water (334 kJ/kg), and ρ and a are as described above.
If we assume that all of the energy delivered from the rainfall is balanced by thawing the soil, and that soils are always saturated at the base of the active layer, then we can calculate the extra thawing of frozen soil per unit rainfall as the ratio of the two terms, H rain /H ALD , which simplify to (in cm thaw/cm rainfall): Supplementary Table 4 provides results from modeling the changes in thaw depth attributable to the thermal inputs of rainfall for each summer between 2013 and 2017. We assumed a 10% error associated with all measurement inputs to this calculation and propagated a cumulative error associated with this expected thaw depth value assuming Gaussian distributions.

DATA AVAILABILITY
The datasets generated during the current study are available from the corresponding author on reasonable request.