Interannual fires as a source for subarctic summer decadal climate

Abstract


Introduction
The Earth's climate trajectory is determined by a combination of internally-generated fluctuations and external forcings such as greenhouse gases, aerosols, land-use changes, and other factors.Uncertainties in aerosol forcings continue to pose challenges to climate hindcasts and projections 1,2 .To reduce the prevailing aerosol forcing uncertainties, newer generations of climate models have incorporated an increasing number of observational constraints on historical aerosol emissions [3][4][5] .As part of this effort, global biomass burning aerosol (BBA) emissions used in the Coupled Model Intercomparison Project Phase 6 (CMIP6) 3 were updated with a merged dataset of satellite observations 6 , fire proxies, and fire model output.In the merged BBA dataset, the incorporation of the satellite-based data has introduced an abrupt onset of large interannual variability over 1997-2014, with substantially weaker interannual variability before 3,7,8 (Supplementary Fig. 1b).Although this approach utilizes best estimates of fire emissions based on available observational products, it implicitly assumes that climate impacts of aerosol forcing will be equivalent as long as the net aerosol fluxes remain the same over an extended period of time (multiple years), regardless of the magnitude of interannual variability.The invocation of this assumption reflects the fact that very little is known about the impacts of interannual variations in BBA emissions on climate.
The discontinuity in interannual variability of BBA forcing in CMIP6 between the satellite observation-based period (1997-2014) and the periods before and after is especially pronounced over boreal North America and Siberia, where interannual variations of fire occurrences and related aerosol emissions are strong (Fig. 1b).Boreal fires are largely affected by mean climatic conditions that are further modulated by different modes of climate variability such as the El Niño-Southern Oscillation, the Arctic Oscillation, or the Pacific Decadal Oscillation [9][10][11][12] .Human intervention is also important through changes in land-use, ignition, fire control, and other perturbations [13][14][15] .
With amplified warming in the northern high latitudes, boreal wildfires have been occurring more frequently and with higher intensity, and the trend is expected to increase under sustained increases in anthropogenic forcings [16][17][18][19][20][21] .Earth system models predict that changes in mean climate alter natural variability across a broad range of timescales 22 , which would in turn affect the variability of fire occurrences.As such, it is timely to ask whether and how changes in year-to-year fluctuations in boreal BBA emissions influence regional and global climate.
Recent studies have documented a northern high latitude warming response to interannually varying BBA emissions in model simulations 7,23 (Supplementary Fig. 1a).
These studies have suggested that nonlinearities in aerosol-cloud processes generate the mean warming response in high BBA variability simulations relative to low BBA variability simulations.However, nonlinearities in atmospheric processes, including aerosol impacts and associated feedbacks, alone cannot account for the sustained warming over decadal timescales (Supplementary Fig. 1c), as the typical lifetime of aerosols from fire emissions is in the order of just a few days.With such a predominantly atmospheric mechanism, atmospheric responses could only fluctuate on interannual timescales and synchronously with the enhanced BBA variability.One may then expect an ocean role in longer timescale responses 23 .Ocean heat storage and release could be considered a viable candidate mechanism for early winter warming over the Arctic Ocean under high BBA variability as the Arctic Ocean warming response maximizes in November (Supplementary Fig. 2) despite near zero BBA emissions during the cold season (Fig. 1a).This early winter warming is similar to seasonally-delayed Arctic warming through ocean heat uptake during the sea ice melting season and heat release during early winter [24][25][26] .Interestingly, however, decadal subarctic land warming (50°N-70°N) (Supplementary Fig. 2) peaks in summer, when no Arctic ocean response is evident.This finding is thereby inconsistent with the idea that Arctic ocean processes are the main driver for the terrestrial responses.
Here, our main scientific objective is to identify the physical mechanism responsible for the decadal summer subarctic land warming response to enhanced interannually varying BBA emissions simulated by the Community Earth System Model version 2 27 Large Ensemble (CESM2-LE 22 ).In CESM2-LE, 50 ensemble members follow the CMIP6 protocols for BBAs (BBA_CMIP6) and a separate group of 50 members are forced by a temporally smoothed version of the CMIP6 BBAs (BBA_smooth) 22 (Supplementary Fig. 1b).The difference in aerosol emissions between the two ensemble groups only exists for interannual scales while conserving the net emissions, and variability differences are large only between 1997-2014 (Methods).This design of the CESM2-LE provides not only a means to characterize the effect of interannual fluctuations in fire emissions on the climate system, but also to identify underlying mechanisms for the apparent decadal response.Our principal new finding is that soil water processes in permafrost provide memory allowing rectification of interannual variability in aerosol forcing to sustain decadal variability in subarctic surface temperature.Additional targeted simulations with prescribed soil moisture confirm that soil water and ice changes in permafrost can modulate subarctic summer temperature.

Regionally and seasonally varying surface temperature response
The long-term annual mean temperature difference (BBA_CMIP6 -BBA_smooth) reveals a typical Arctic amplification pattern 28,29 with more warming at higher latitudes (Supplementary Fig. 1a).This annual mean temperature difference, however, is due to regionally distinct seasonal changes.Whereas Arctic (70°N-90°N) temperatures show maximum warming in November, subarctic regions (50°N-70°N) exhibit strongest warming in July (Supplementary Fig. 2).In addition to interannual fluctuations caused by BBA, the high latitude surface temperature response is also characterized by decadalscale climate shifts (Fig. 2e and Supplementary Fig. 1c).Boreal fires typically occur during summer, which explains why the summertime BBA emission anomalies between BBA_CMIP6 and BBA_smooth also peak in summer (Fig. 1a).Due to the short lifetime (~ a few days) of aerosols in the atmosphere, the anomalous summertime input of BBAs into the atmosphere in BBA_CMIP6 relative to BBA_smooth quickly disappears and does not remain over the following summer season, as shown in the atmospheric BBA burden difference (Supplementary Fig. 3).This implies that the direct and indirect BBA forcing difference would also be active only during the BBA emission season and on interannual timescales.In the two groups of 50-member ensemble simulations, the net time-integrated aerosol differences are negligible on decadal timescales.Nevertheless, a pronounced decadal warming anomaly is apparent (Fig. 2e and Supplementary Fig. 1c), which suggests the existence of nonlinear climate rectification processes that translate the interannual forcing into a decadal climate signal.

The trigger: interannual direct and indirect effects of aerosols
We first investigate how summertime BBA emissions affect concomitant atmospheric conditions.To isolate the aerosol impacts while excluding decadal scale responses, composites of BBA_CMIP6 fields for 4 summers with high BBA emissions are subtracted from composites of BBA_CMIP6 fields for 8 summers with low BBA emissions (the mean for 8 low emission summers minus the mean for 4 high emission summers) (Figs.1c-f and Supplementary Fig. 4).The high and low emission years in BBA_CMIP6 are selected based on the BBA differences relative to BBA_smooth during summer from May through September (the 8 lowest and 4 highest were chosen).These  4b).Decreased relative humidity (Supplementary Fig. 4c) suggests that the reduced cloud cover during low BBA emission years is caused not only by reduced CCNs, but also through atmospheric feedbacks.We note that while emissions of light-absorbing aerosols such as black carbon and brown carbon from fires can impose a warming effect 30,31 , it has been estimated by observation-based studies that the net combined direct radiative forcing by all types of BBAs is negative 2,32,33 and that the net forcing is also negative [34][35][36] , which is consistent with our results.

The rectifier: nonlinear responses of clouds and soil water
The strong decadal warming over land in summer suggests that soil-related mechanisms may provide the memory and nonlinearity to rectify the interannual aerosol forcing to the decadal warming.We can rule out a role of ocean for the summer subarctic warming intensification, as the ocean surface rather absorbs heat from the atmosphere (increased heat uptake) in summer in BBA_CMIP6 relative to BBA_smooth (Supplementary Fig. 5), demonstrating that the decadal summer temperature increase over land is not caused by heat transport from nearby oceans.The time-depth evolution of BBA_CMIP6 minus BBA_smooth fields of soil ice, soil liquid, and their sum as net soil moisture over 50°N-70°N (Fig. 2) reveals that thawing of upper soil permafrost, triggered by summers of low aerosol emissions, propagates through the deeper soil layer over decadal timescales.The increased liquid water in deeper soil indicates that water drains easily through more porous soil associated with permafrost thawing 37 , contributing to the decadal drying tendency of the upper soil layer.The decadal decrease in soil moisture and retreat of permafrost in response to interannual summertime BBA fluctuations indicates that thawing and freezing of soil are not symmetric between negative and positive BBA conditions.Clearly, the upper permafrost thawing and soil moisture loss during low BBA emission years exceed the upper permafrost freezing and soil moisture increase during high emission years, thereby creating a hysteresis effect, which accumulates in time, generating decadal-scale climate responses.
Several processes are involved in the nonlinear response of soil moisture to changes in BBAs.Similar to the nonlinear relationship between anthropogenic aerosols and their indirect forcing 38 , BBA forcing in our simulations exhibits stronger sensitivity under a cleaner environment (reduced aerosol conditions).Scatter plots of CCN versus BBA reveal negative curvature, indicating that CCN formation is more sensitive to negative BBA anomalies than to positive anomalies (Fig. 3a).The same argument holds for cloud droplet sensitivity to CCN (Fig. 3b) and to the inverse relationship between surface shortwave fluxes and cloud droplet concentration (Fig. 3c).As a consequence of these nonlinearities, negative BBA anomalies can generate a stronger surface shortwave response than positive BBA anomalies (Fig. 3d).
Nonlinearities also exist in other processes related to phase changes of water such as sensitivity of soil ice melting/freezing to changes in radiation or temperature, although it may be challenging to identify its isolated effect in our simulations as the phase change processes involving water occur over the full range of surface to deeper soils, with gradual heat and water transfer over multiple seasons (response time lag as a function of depth).The largest soil ice differences between the two ensemble groups are found at subsurface depth of 0.9-1.9m (Supplementary Fig. 6a).We have found a lagged relationship between the net surface shortwave radiation to ice changes (tendency) at this subsurface depth of 0.9-1.9m where the maximum ice differences are found.The correlation between the surface shortwave flux and soil ice melting tendency at this layer is maximized with a time lag of 2 months (correlation coefficient of 0.82).The scatter diagram in Fig. 3e demonstrates that the subsurface soil ice melting is more sensitive to higher surface shortwave fluxes than to lower shortwave fluxes.This is largely because ice melts only at temperatures > 0°C, providing strong nonlinearity.
Albedo feedbacks involving snow cover changes can also add nonlinearity to the response 39 .This, however, is not the main source of nonlinearity in summer as both ensemble groups become nearly snow-free over the most subarctic regions by July (Supplementary Fig. 7).Thus snow depth differences and albedo differences between the two ensemble groups are mostly not significant except for small areas adjacent to the Arctic (Supplementary Fig. 9a-b).The surface temperature response to shortwave fluxes is almost linear on average over land between 50°N-70°N, although small snow covered areas still exhibit a nonlinear behavior in surface temperature responses to shortwave flux variations in summer (Supplementary Figs.8c-d).
Overall, the combined effects of nonlinearities from aerosol-cloud-radiation interactions and liquid-ice phase changes in soil result in a nonlinear response of soil ice to BBA emissions (Fig. 3f).This promotes net thawing of permafrost on decadal scales in BBA_CMIP6.As a consequence of permafrost thawing, an accelerated drainage of water gradually deprives the upper soil of moisture, providing the long-term climate memory (or reddening of the temporal spectrum).

The climate memory: decadal summer subarctic warming
The hydrological processes occurring in the upper soil layer have several types of impacts on surface temperature.First, the overall drying of the upper soil layer in BBA_CMIP6 compared to BBA_smooth (Fig. 4b) lowers the heat capacity of the soil surface, inducing higher surface temperatures and associated positive feedbacks to nearsurface atmospheric warming in summer.The heat capacity of soil is determined by the volumetric fraction of minerals, organic matters, and water content [40][41][42] .Porous soil with low soil wetness is filled with air, which has lower heat capacity than liquid or ice water.
The lower heat capacity of the drier soil more readily facilitates surface warming in summer.Second, changes in evaporative fluxes modify surface temperature by absorbing or releasing latent heat.The Bowen ratio, defined as the ratio of sensible to latent heat flux, is a good indicator of this effect 43 .For example, a high Bowen ratio indicates that the surface releases heat to the atmosphere more through sensible than latent heat fluxes due to limited surface moisture to evaporate, raising surface temperature.Although not all areas over the subarctic land have the additional warming effect due to reduced evaporation, Supplementary Fig. 9h suggests that temperature damping effects by evaporation are suppressed over the regions where surface moisture is strongly depleted in BBA_CMIP6 relative to BBA_smooth (drier soil moisture in Fig. 4b and higher positive Bowen ratio differences Supplementary Fig. 9e).However, we note that, except for those extremely dry areas, the evaporative effect does not seem to be the dominant contributor to the overall decadal warming (surface temperatures are similar both over areas with positive Bowen ratio changes and with negative changes in Supplementary Figs.9f-g).This is likely because moisture can still be supplied by thawed surface soil during summer, as inferred by the distributions of sensible and latent heat fluxes in BBA_smooth and their change in BBA_CMIP6 (Supplementary Figs.9a-d).Another impact of the hydrological processes on soil temperatures is through latent heat uptake and release by melting and freezing of soil ice.The decreased upper soil ice content causes less melting from spring to summer (Supplementary Fig. 10), absorbing less latent heat of fusion, and thus sustaining less cooling (or more warming) at the surface.In other words, the presence of smaller amounts of ice through the transition from spring to summer contributes to additional surface heating in July.
The decadal warming in BBA_CMIP6 relative to BBA_smooth mediated by these soil moisture processes is further enhanced by immediate cloud feedbacks.The atmospheric warming induced by soil interactions over the subarctic land domain is more intensified towards the surface (Supplementary Fig. 11).Although specific humidity also increases due to enhanced surface evaporation, the increase in atmospheric moisture is nevertheless modest relative to the temperature increase, thereby lowering relative humidity.This in turn reduces cloudiness throughout the mid-to lower-troposphere, with more cloud reduction occurring in the lower reaches.Fewer clouds allow more incoming solar radiation to reach the surface, thereby heating the surface and lower atmosphere further (Fig. 4 and Supplementary Fig. 11).The presence of this atmospheric feedback is supported by the positive shift in anomalous surface shortwave fluxes even for moderate positive anomalous BBA emissions (see the dotted lines in Fig. 3d).

Prescribed soil moisture experiment
The soil moisture feedback on summer surface climate is demonstrated by conducting additional two sets of simulations with prescribed soil moisture calculated from CESM2-LE.The purpose of this experiment is to quantify to what extent surface climate differences in BBA_CMIP6 and BBA_smooth can be explained by soil moisture differences.The two simulation sets were all restarted from the year of 2000 and 2010 from original 40 BBA_smooth members, and ran for a year under the same forcings including BBA emissions.(The choice of the years is to utilize existing restart files.)The only difference between the two sets is that for one set we prescribed daily mean soil moisture fields (liquid water and ice) obtained from CESM2-LE climatology across the 50 BBA_CMIP6 ensemble members and across 15 years (2000-2014), while for the other set from the 50 BBA_smooth members.
The surface temperature difference between the two sets in Fig. 5 demonstrates that reduced soil moisture induces substantial surface warming.As soil moisture is more depleted at higher latitudes in BBA_CMIP6 compared to BBA_smooth (Fig. 4b), the land temperature response to soil moisture changes becomes larger closer to the Arctic (Fig. 5).The contribution of surface warming solely by soil drying (Fig. 5) relative to all contributions (Fig. 4a) is 22% over 50°N-60°N and 43% over 60°N-70°N.The soil drying-induced warming is especially severe over Siberia (Fig. 5).This prescribed soil moisture experiment confirms that the overall surface warming in BBA_CMIP6 relative to BBA_smooth in Fig. 4a is partly due to asymmetric response to reduced and enhanced BBA fluctuations on interannual scales, and partly due to permafrost thawing-induced soil moisture drying on decadal scales.

Discussion
Climate models are highly sensitive to a variety of aerosol forcings, including those from biomass burning.However, the sensitivity of the climate system to changes in the interannual variance of aerosols has not previously been considered in experimental designs for coordinated simulations such as CMIP6.Given that wildfires vary seasonally and from year-to-year, we have chosen to address this important question over the satellite observation period of fire emissions (1997-2014) using two groups of large ensemble simulations with 50 members each.Our large ensemble approach using a single Earth system model, comprised of 50 members forced by pronounced interannual BBA variations during the satellite era following the CMIP6 protocol, and 50 members in which these interannual BBA variations are smoothed out 22 , enables us to identify the forced response signal that can otherwise be easily obscured by strong internal variability.The 100 ensemble members for the historical period spanning 1850-1990     forced with identical BBA fluxes facilitate quantification of residuals from internal variability that are not completely neutralized even after averaging 50 members (Methods), providing higher confidence in identifying a signal above natural variability noise for the time interval of interest.
We found here that there is a strong decadal rectified response to interannual BBA variability over the NH subarctic land regions, with a net summer warming in the presence of elevated interannual variability of biomass burning emissions.This mean state response is mediated through interactions with subarctic soil moisture, resulting in distinct summertime footprints over the NH high latitude land at the decadal timescale.
Although recent studies 7,8,23 have identified a similar warming response in the annual mean to interannual fluctuations of BBA forcing, their proposed mechanistic pathways for the climate response are less clear and distinct from the rectifier mechanisms identified here.The earlier studies have hypothesized that asymptotic behaviors of aerosol effects on clouds with skewness in the emissions lead to a net reduction in low cloud amount, which increases incoming shortwave radiation over northern high latitudes 7 .However, their proposed mechanism cannot explain the decadal scale climate response, as aerosol forcings along with the invoked atmospheric processes cannot sustain longerterm rectified responses (see the schematic timeseries of the nonlinearity-only response in Fig. 6).Our study elucidates that indirect aerosol effects are only one of the triggers of decadal permafrost thawing, through which the decadal summer land warming is mediated (Fig. 6).With a simple conceptual framework for the combined effects of nonlinearity (by atmospheric processes and soil ice melting) and reddening (by soil wetness serving as a climate memory), we can explain a decadal warming response during the high BBA variability period of 1997-2014 and the recovery of the climate system afterward (Fig. 6 and Supplementary Fig. 12).
The rectified responses identified for surface variables such as temperature and soil moisture are largely limited to the time interval of the perturbation (aerosol variance modulations) itself, but longer-term responses do persist for the subsurface and deeper soil layers in permafrost.The total loss of soil ice in BBA_CMIP6 is significant and its impact lasts for decades (Supplementary Fig. 13).Depending on the season, BBA_CMIP6 loses soil ice by a few percent up to ten percent relative to BBA_smooth (Supplementary Fig. 14).In reality such a permafrost perturbation could potentially contribute to the release of methane, promoting additional warming.The permafrost thawing-induced warming is also expected to modulate wildfire activity that can have feedback effects through ensuing carbon and aerosol releases.It is our hope that the rectifier mechanism identified in this study motivates improved representation of processes that impact BBA emissions, including through the inclusion of interactive fires in models.
Our analyses have revealed that modulations of variance in aerosols and shortwave radiation can be rectified by inherent nonlinearities and slow responses (reddening) in the soil system (Fig. 6).This new conceptual framework may have further implications for understanding the generation of long-term climate variability in response e.g. to year-toyear changes in fires, volcanic activity or even randomly-occurring interannual changes in Arctic and subarctic cloudiness.The soil system serves as a memory for subarctic summer climate in this study, but a similar conceptual model with Arctic sea ice and/or ocean as a memory may also be applied to other components of the climate system.Our results, which demonstrate mechanistically that spurious modulations of variability in aerosol fluxes can lead to biases in the forced response, should be considered in the design of aerosol forcing protocols for the upcoming Coupled Model Intercomparison Project Phase 7.

CESM2-LE simulations
The 100 ensemble members comprising the CESM2-LE are divided into two subgroups that are distinguished by their respective BBA emissions 22 .The first group of 50 members follows the CMIP6 protocols for BBA and thereby includes strong interannual variations during 1997-2014 due to the use of satellite-based Global Fire Emissions Database (GFED) observations 3,6 (Supplementary Fig. 1b).The second group with 50 members has much smaller interannual variability, as the prescribed BBA fields are smoothed by applying 11-year running means at each grid point from the CMIP6 BBAs.
The 11-year running mean is performed for each month of the year separately to retain the annual cycle of BBA emissions.These two groups are referred to as BBA_CMIP6 and BBA_smooth, respectively, in this study.Net aerosol emissions are nearly conserved between the two groups of ensembles (difference < 0.35 % over 30°N-90°N and < 0.1% globally for 1997-2014), and the aerosol difference only exists at the interannual scale.
Detailed descriptions of the dataset can be found in 22,27 .To focus on responses over the high latitude NH, the BBA emissions over 30°N-90°N is considered in this study.We note that the strongest contribution is emitted over 50°N-70°N as shown in Fig. 1b, and results are not sensitive to the choice of a BBA emission latitude band between 30°N-90°N and 50°N-70°N.

Estimating significance of composites for low -high emission summers
A bootstrap method is used for estimating uncertainty ranges for the composites of low emission summers minus high emission summers in BBA_CMIP6 in Fig. 1 and Supplementary Fig.   (1997, 1999, 2000, 2001, 2005, 2007, 2009, 2011) minus composites of 4 high BBA emission years (1998, 2002, 2003, 2012)    CCN and cloud droplet concentrations are vertically integrated throughout the troposphere.The BBA flux is integrated over 30°N-90°N, and all other simulated values are averaged over the land domain between 50°N-70°N.In (a-f), each scatter dot corresponds to a 1-month value for June (blue), July (orange), and August (green).In (e), the ice tendency is calculated as a change in subsurface ice during 2 months following surface shortwave flux for a given month.In (f), each scatter dot corresponds to June-July-August (JJA) BBA flux versus a soil ice change from a previous year.Since the soil ice response to summer BBA emissions propagates to deeper soil over following months as depicted in Fig. 2b, mean soil ice for summer to the following winter (June to next February) was taken for each given year (e.g., the dot for 1998 corresponds to the BBA flux in JJA 1998 versus the mean soil ice for June 1998 -February 1999 minus the mean soil ice for June 1997 -February 1998).Spring values are excluded in the soil ice calculation given the ambiguity as to whether spring soil ice changes are influenced by the previous summer or the same year's spring emissions (Note that spring BBA emissions are irregular as shown in Fig. 1a).The black curve is a least-square seconddegree polynomial fit from filled dots.Unfilled dots are considered as outliers whose squared error is greater than the three-standard deviation of all squared errors.
composites represent interannual responses to decreased prescribed BBA fluxes.Aerosols modify the temperature structure of the atmosphere and surface by changing radiation through direct and indirect effects.For decreased BBA conditions, the associated reductions in aerosol optical depth (AOD) and cloud condensation nuclei (CCN), whose patterns match the main BBA source regions over Siberia and northern Canada, enhance surface shortwave fluxes, and thereby heat the surface (Figs.1b-f).The composites in Supplementary Figs.4a-b further illustrate that over subarctic land regions the direct effect by reduced aerosols (clear-sky shortwave flux in Supplementary Fig. 4a) is comparable to the indirect effects due to reduced clouds (shortwave flux change by clouds in Supplementary Fig. 4. A bootstrapped field is calculated by subtracting the mean of BBA_CMIP6 fields for 4 randomly selected years from the mean of BBA_CMIP6 fields for 8 randomly selected years.Random years are selected between 1997-2014 without replacement, and all calculations are for May to September (MJJAS).Values within one

FiguresFig. 1 .
Figures in BBA_CMIP6 for (c) aerosol optical depth (AOD), (d) cloud condensation nuclei (CCN), (e) the net shortwave flux at the surface (positive downward), and (f) surface temperature.The 8 low and 4 high emission years in BBA_CMIP6 are selected based on the BBA emission differences, shown in (a), for summer months for May to September (MJJAS).The composite difference maps in (c-f) are for the same summer months (MJJAS).Stippling denotes non-significant areas (see Methods).

Fig. 2 .
Fig. 2. Evolution of soil ice and surface water and temperature in response to BBA fluctuations over 50°N-70°N land.(a) Monthly time series of BBA flux difference over 30°N-90°N (same as Fig. 1a).Monthly differences of (b) soil ice, (c) soil liquid, and (d) net soil water (ice + liquid) show that melting of ice propagates into deeper soil and that upper layer soil moisture is diminished in BBA_CMIP6.(e) Monthly time series of differences of surface temperature and (f) upper 10 cm soil moisture show decadal scale changes.Dot colors in time series represent months (same as in Fig. 1a) with enlarged pink dots for July.Unit is kg m -2 for (b-d) and (f).Stippling in (b-d) denotes nonsignificant areas whose values are within one standard deviation of an internal variability range (see Methods).All values correspond to differences between the two ensemble groups (BBA_CMIP6 -BBA_smooth).

Fig. 3 .
Fig. 3. Nonlinear sensitivity of BBA-cloud-radiation and BBA-soil ice.Relationship between anomalies of (a) prescribed BBA flux and CCN, (b) CCN and cloud droplet concentration, (c) cloud droplet concentration and net surface shortwave flux by clouds (positive downward), (d) BBA and net surface shortwave flux, (e) net surface shortwave flux and subsurface (0.9-1.9 m depth) soil ice tendency, and (f) BBA flux and total soil ice tendency.All values are taken from BBA_CMIP6 -BBA_smooth over 1997 -2014.

Fig. 4 .
Fig. 4. July responses in BBA_CMIP6 -BBA_smooth.Difference maps of (a) temperature, (b) upper 10 cm soil moisture, (c) net surface shortwave flux (positive downward), (d) low-level (> 700 hPa) relative humidity, and (e) low-level cloud fraction in July over 2000-2014.Stippling denotes non-significant areas whose values are within one standard deviation of an internal variability range (see Methods).

Fig. 5 .
Fig. 5. Surface temperature response from the prescribed soil moisture experiment using CESM2.Color shows surface temperature difference in July between the two prescribed soil moisture cases of BBA_CMIP6 and BBA_smooth.Except for prescribing different daily soil liquid water and ice fields, all other simulation conditions have remained the same for the two simulation sets.See simulation details in Results.

Fig. 6 .
Fig. 6.Schematic of atmosphere-land-permafrost coupling and generation of decadal variability from interannual forcing.For each given year, reduced (enhanced) aerosol emissions from biomass burning induces net warming (cooling) of the atmosphere and surface through direct and indirect aerosol impacts.As a consequence of nonlinear sensitivities in aerosol-cloud-radiation processes and soil ice melting, interannual changes in aerosols result in a decadal scale permafrost loss, promoting drainage of upper soil water to the deeper soil layer.The decadal decrease in upper soil moisture by the melt-water drainage induces decadal warming of the surface and lower troposphere that is intensified by cloud feedbacks through changes in relative humidity.The bottom panels illustrate three conceptual response models when (upper left) symmetric interannual fluctuation of aerosols is present, mimicking summertime aerosol