Supplementary Information forInterannual fires as a source for subarctic summer decadal climate
variability mediated by permafrost thawing
Ji-Eun Kim, Ryohei Yamaguchi, Keith B. Rodgers, Axel Timmermann, Sun-Seon Lee,Karl Stein, Gokhan Danabasoglu, Jean-Francois Lamarque, John T. Fasullo, Clara Deser,Nan Rosenbloom, Jim Edwards, Malte F. Stuecker

Corresponding author: Ji-Eun Kim
Email: jieunkim@pusan.ac.kr

npj Climate and Atmospheric Science1Supplementary Figure 1. Decadal mean surface temperature response tointerannually varying BBA emissions. (a) Mean annual temperature difference(BBA_CMIP6 – BBA_smooth) for 2000-2014. Stippling denotes non-significant areas
whose values are within one standard deviation of the internal variability range (seeMethods). (b) Time series of 30°N-90°N annual BBAs prescribed in BBA_CMIP6(orange) and BBA_smooth (blue) simulations. (c) Time series of annual temperature overthe land domain between 50°N-70°N for the two ensemble groups. Shading represents
one standard deviation from 50 ensemble members for each group.

2Supplementary Figure 2. Regionally and seasonally distinct decadal temperatureresponses. Temperature differences between the two ensemble groups (BBA_CMIP6 –BBA_smooth) over 2000-2014 for each month. The cyan contour corresponds to 2°C.
Stippling denotes non-significant areas (see Methods).

Supplementary Figure 3. Monthly evolution of differences (BBA_CMIP6 –BBA_smooth) in atmospheric BBA burden. Monthly time series of BBA column
burden differences over the land domain between 50°N-70°N. Dot colors denote months
as in Fig. 1a.
3

Supplementary Figure 4. Aerosol direct and indirect effects. Composites of 8 low
BBA emission years (1997, 1999, 2000, 2001, 2005, 2007, 2009, 2011) minus
composites of 4 high emission years (1998, 2002, 2003, 2012) in BBA_CMIP6 for (a)clear-sky shortwave flux at the surface, (b) changes in surface shortwave flux by clouds,(c) and low-level (> 700 hPa) relative humidity. The downward direction is positive for(a) and (b). The selected composite years are the same as in Fig. 1c-f, and all values arefor MJJAS. Stippling denotes non-significant areas (see Methods).

Supplementary Figure 5. July difference (BBA_CMIP6 – BBA_smooth) map of net
heat uptake at the surface averaged over 2000-2014. The net heat uptake is estimatedas residuals of net surface shortwave, longwave, sensible heat, and latent heat fluxes.
Positive values represent anomalous heat entering (downward) the surface, and stippling
denotes non-significant areas (see Methods).
4

Supplementary Figure 6. Soil ice melting in response to increases in surface
shortwave flux. (a) Anomalous soil ice in BBA_CMIP6 relative to BBA_smooth as afunction of depth. The orange dots correspond to soil model layers in CLM5. The largest
ice melting occurs at depth between 0.9-1.9 m (shading). (b) Correlation coefficients
between net surface shortwave flux and soil ice tendency at soil layers of 0.9-1.9 m as afunction of lag month.

Supplementary Figure 7. Monthly mean snow cover for 2000-2014 in BBA_smooth.

5Supplementary Figure 8. Snow cover and albedo feedback. July difference(BBA_CMIP6 – BBA_smooth) maps of (a) snow depth and (b) Albedo. Stippling
denotes non-significant areas (see Methods). Relationship between net surface shortwave
flux and surface temperature averaged over (c) 50°N-70°N land and (d) areas wherealbedo difference < -0.01. For scatters and the polynomial fit in (c-d), the same methods
are used as for Fig. 3a-d.

Supplementary Figure 9. Sensible and latent heat fluxes, and their impact on
surface temperature. The mean (a) sensible and (b) latent heat flux in BBA_smooth,and differences (BBA_CMIP6 – BBA_smooth) of (c) sensible and (d) latent heat flux,and (e) Bowen ratio (sensible to latent heating). Monthly time series of surface
temperature over areas where (f) Bowen ratio > 0, and (g) Bowen < 0, and (h) Bowen
ratio > 0.07 over 50°N-70°N land. Dot colors denote months.

6

Supplementary Figure 10. Annual cycle of upper 12 cm soil ice in BBA_CMIP6(orange) and BBA_smooth (blue), and their difference (black) over 50°N-70°N for
2000-2014. The grey shading represents ranges of non-significant values (see Methods).

Supplementary Figure 11. Atmospheric responses in July over the land domain
between 50°N-70°N for 2000-2014. Vertical profiles (left to right) of temperature,
specific humidity, relative humidity, and cloud fraction for BBA_CMIP6 (orange) andBBA_smooth (blue), and their differences (black). The grey shading represents ranges of
non-significant values (see Methods).
7

Supplementary Figure 12. Three conceptual model responses to (upper left) aerosol
fluctuations in BBA_CMIP6 – BBA_smooth. (lower left) A response varies at the sameinterannual timescale as the emission timescale when there is only nonlinearity inatmospheric processes. (upper right) A response fluctuates interannually with a gradualdecay of a signal (reddening) after a symmetric linear response to an aerosol pulse whena climate system has a memory effect. The decay timescale of 2 years is used in thisexample. (lower right) When both nonlinearity and reddening are present, a response issimilar to the temperature response shown Fig. 2e. The upper left aerosol anomalies are
estimated from annual mean BBA differences between BBA_CMIP6 and BBA_smooth,and the time-integrated value is set to zero. Units are arbitrary.

Supplementary Figure 13. Same as Fig. 2b, but for an extended period. Unit is kg m2.

8

Supplementary Figure 14. Loss of soil ice. Percentage change of vertically-integrated
soil ice over the entire model layer (up to 8.6 m depth (1)) in BBA_CMIP6 relative toBBA_smooth over 50°N-70°N, with colors distinguishing months.

Supplementary Figure 15. Monthly time series of surface temperature difference(BBA_CMIP6 – BBA_smooth) over the land domain between 50°N-70°N for theperiod of 1850-2100. Surface temperature differences between the two 50-member
ensemble groups nevertheless exhibit residuals even before 1990 due to non-zerocancellation of internal variability. Colors denote different months.

9