1

Interannual fires as a source for subarctic summer decadal climate

2

variability mediated by permafrost thawing

3
4

Authors

5

Ji-Eun Kim*a,b, Ryohei Yamaguchic, Keith B. Rodgersa,b, Axel Timmermanna,b, Sun-Seon

6

Leea,b, Karl Steina,b, Gokhan Danabasoglud, Jean-Francois Lamarqued, John T. Fasullod,

7

Clara Deserd, Nan Rosenbloomd, Jim Edwardsd, Malte F. Stueckere

89Affiliations

10

a

11

b

12

c

13

Technology, Kanagawa, Japan.

14

d

15

e

16

and Earth Science and Technology, University of Hawaiʻi at Mānoa, Honolulu, HI, USA.

Center for Climate Physics, Institute for Basic Science, Busan, South Korea.
Pusan National University, Busan, South Korea.

Research Institute for Global Change, Japan Agency of Marine-Earth Science andNational Center for Atmospheric Research, Boulder, CO, USA.

Department of Oceanography and International Pacific Research Center, School of Ocean

17
18*Correspondence to: jieunkim@pusan.ac.kr

19
20
21

npj Climate and Atmospheric Science

22
23

1

24

Abstract25

Climate model simulations run under the Coupled Model Intercomparison Project Phase26

6 (CMIP6) use an inhomogeneous biomass burning aerosol (BBA) emission dataset,

27

which exhibits pronounced interannual variability after 1997 due to the infusion of28

satellite data. Here we quantify the impact of the sudden change in aerosol variance on

29

subarctic summer climate and propose a mechanism connecting the interannual forcing to30

longer timescale climate responses by using the Community Earth System Model version

31

2 Large Ensemble (CESM2-LE). Considering 50 ensemble members with a 199732

transition to high BBA variability (following the standard CMIP6 protocol) and 50

33

members with consistently lower fluctuations, interannual aerosol variations between

34

1997-2014 generate interannual atmospheric responses synchronous with BBA

35

fluctuations as well as pronounced decadal-scale subarctic land warming. Under strong

36

interannually varying aerosol amounts, soil ice melting during reduced aerosol years

37

exceeds freezing during enhanced aerosol years due to nonlinearities in aerosol-cloud

38

interactions and in soil water phase change processes, resulting in a decadal accelerated

39

reduction of permafrost and increased soil water drainage. By conducting idealized

40

simulations with prescribed soil moisture, we further demonstrate that soil moisture

41

changes in permafrost feedback on summer surface temperatures. In more general terms

42

our proposed mechanism reveals that climate variability at shorter timescales can induce

43

climate responses on longer timescales if nonlinearity and memory effects coexist.

44
45

Introduction

46

The Earth’s climate trajectory is determined by a combination of internally-generated

47

fluctuations and external forcings such as greenhouse gases, aerosols, land-use changes,

48

and other factors. Uncertainties in aerosol forcings continue to pose challenges to climate

49

hindcasts and projections 1,2. To reduce the prevailing aerosol forcing uncertainties,

50

newer generations of climate models have incorporated an increasing number of51

observational constraints on historical aerosol emissions 3-5. As part of this effort, global

52

biomass burning aerosol (BBA) emissions used in the Coupled Model Intercomparison
2

53

Project Phase 6 (CMIP6) 3 were updated with a merged dataset of satellite observations 6,

54

fire proxies, and fire model output. In the merged BBA dataset, the incorporation of the55

satellite-based data has introduced an abrupt onset of large interannual variability over56

1997-2014, with substantially weaker interannual variability before 3,7,8 (Supplementary

57

Fig. 1b). Although this approach utilizes best estimates of fire emissions based on

58

available observational products, it implicitly assumes that climate impacts of aerosol

59

forcing will be equivalent as long as the net aerosol fluxes remain the same over an

60

extended period of time (multiple years), regardless of the magnitude of interannual

61

variability. The invocation of this assumption reflects the fact that very little is known

62

about the impacts of interannual variations in BBA emissions on climate.

63

The discontinuity in interannual variability of BBA forcing in CMIP6 between the64satellite observation-based period (1997-2014) and the periods before and after is

65

especially pronounced over boreal North America and Siberia, where interannual

66

variations of fire occurrences and related aerosol emissions are strong (Fig. 1b). Boreal

67

fires are largely affected by mean climatic conditions that are further modulated by68

different modes of climate variability such as the El Niño-Southern Oscillation, the69

Arctic Oscillation, or the Pacific Decadal Oscillation 9-12. Human intervention is also

70

important through changes in land-use, ignition, fire control, and other perturbations 13-15.

71

With amplified warming in the northern high latitudes, boreal wildfires have been

72

occurring more frequently and with higher intensity, and the trend is expected to increase73

under sustained increases in anthropogenic forcings 16-21. Earth system models predict

74

that changes in mean climate alter natural variability across a broad range of timescales

75

22

76

ask whether and how changes in year-to-year fluctuations in boreal BBA emissions77

influence regional and global climate.

78

Recent studies have documented a northern high latitude warming response to79

interannually varying BBA emissions in model simulations 7,23 (Supplementary Fig. 1a).

80

These studies have suggested that nonlinearities in aerosol-cloud processes generate the81mean warming response in high BBA variability simulations relative to low BBA, which would in turn affect the variability of fire occurrences. As such, it is timely to3

82variability simulations. However, nonlinearities in atmospheric processes, including

83

aerosol impacts and associated feedbacks, alone cannot account for the sustained

84

warming over decadal timescales (Supplementary Fig. 1c), as the typical lifetime of85

aerosols from fire emissions is in the order of just a few days. With such a predominantly

86

atmospheric mechanism, atmospheric responses could only fluctuate on interannual

87

timescales and synchronously with the enhanced BBA variability. One may then expect88

an ocean role in longer timescale responses 23. Ocean heat storage and release could be89

considered a viable candidate mechanism for early winter warming over the Arctic Ocean

90

under high BBA variability as the Arctic Ocean warming response maximizes in91

November (Supplementary Fig. 2) despite near zero BBA emissions during the cold

92

season (Fig. 1a). This early winter warming is similar to seasonally-delayed Arctic

93

warming through ocean heat uptake during the sea ice melting season and heat release94

during early winter 24-26. Interestingly, however, decadal subarctic land warming (50°N-

9570°N) (Supplementary Fig. 2) peaks in summer, when no Arctic ocean response is96

evident. This finding is thereby inconsistent with the idea that Arctic ocean processes are97

the main driver for the terrestrial responses.

98

Here, our main scientific objective is to identify the physical mechanism responsible for99

the decadal summer subarctic land warming response to enhanced interannually varying

100

BBA emissions simulated by the Community Earth System Model version 2 27 Large

101

Ensemble (CESM2-LE 22). In CESM2-LE, 50 ensemble members follow the CMIP6

102

protocols for BBAs (BBA_CMIP6) and a separate group of 50 members are forced by a103

temporally smoothed version of the CMIP6 BBAs (BBA_smooth) 22 (Supplementary Fig.

1041b). The difference in aerosol emissions between the two ensemble groups only exists for

105

interannual scales while conserving the net emissions, and variability differences are

106large only between 1997-2014 (Methods). This design of the CESM2-LE provides not107

only a means to characterize the effect of interannual fluctuations in fire emissions on the

108

climate system, but also to identify underlying mechanisms for the apparent decadal

109

response. Our principal new finding is that soil water processes in permafrost provide

110

memory allowing rectification of interannual variability in aerosol forcing to sustain4

111

decadal variability in subarctic surface temperature. Additional targeted simulations with

112

prescribed soil moisture confirm that soil water and ice changes in permafrost can

113

modulate subarctic summer temperature.

114
115

Results

116

Regionally and seasonally varying surface temperature response

117

The long-term annual mean temperature difference (BBA_CMIP6 – BBA_smooth)118

reveals a typical Arctic amplification pattern 28,29 with more warming at higher latitudes

119(Supplementary Fig. 1a). This annual mean temperature difference, however, is due to120

regionally distinct seasonal changes. Whereas Arctic (70°N-90°N) temperatures show121

maximum warming in November, subarctic regions (50°N-70°N) exhibit strongest

122

warming in July (Supplementary Fig. 2). In addition to interannual fluctuations caused by

123

BBA, the high latitude surface temperature response is also characterized by decadal-

124

scale climate shifts (Fig. 2e and Supplementary Fig. 1c). Boreal fires typically occur

125

during summer, which explains why the summertime BBA emission anomalies between

126

BBA_CMIP6 and BBA_smooth also peak in summer (Fig. 1a). Due to the short lifetime127(~ a few days) of aerosols in the atmosphere, the anomalous summertime input of BBAs

128

into the atmosphere in BBA_CMIP6 relative to BBA_smooth quickly disappears and129

does not remain over the following summer season, as shown in the atmospheric BBA

130burden difference (Supplementary Fig. 3). This implies that the direct and indirect BBA

131

forcing difference would also be active only during the BBA emission season and on

132

interannual timescales. In the two groups of 50-member ensemble simulations, the net

133

time-integrated aerosol differences are negligible on decadal timescales. Nevertheless, a134

pronounced decadal warming anomaly is apparent (Fig. 2e and Supplementary Fig. 1c),135

which suggests the existence of nonlinear climate rectification processes that translate the136interannual forcing into a decadal climate signal.

137
5

138

The trigger: interannual direct and indirect effects of aerosols

139

We first investigate how summertime BBA emissions affect concomitant atmospheric140

conditions. To isolate the aerosol impacts while excluding decadal scale responses,

141

composites of BBA_CMIP6 fields for 4 summers with high BBA emissions are

142

subtracted from composites of BBA_CMIP6 fields for 8 summers with low BBA

143emissions (the mean for 8 low emission summers minus the mean for 4 high emission

144

summers) (Figs. 1c-f and Supplementary Fig. 4). The high and low emission years in

145

BBA_CMIP6 are selected based on the BBA differences relative to BBA_smooth during

146

summer from May through September (the 8 lowest and 4 highest were chosen). These

147

composites represent interannual responses to decreased prescribed BBA fluxes. Aerosols

148

modify the temperature structure of the atmosphere and surface by changing radiation

149

through direct and indirect effects. For decreased BBA conditions, the associated150

reductions in aerosol optical depth (AOD) and cloud condensation nuclei (CCN), whose151

patterns match the main BBA source regions over Siberia and northern Canada, enhance152

surface shortwave fluxes, and thereby heat the surface (Figs. 1b-f). The composites in

153

Supplementary Figs. 4a-b further illustrate that over subarctic land regions the direct154

effect by reduced aerosols (clear-sky shortwave flux in Supplementary Fig. 4a) is

155

comparable to the indirect effects due to reduced clouds (shortwave flux change by

156

clouds in Supplementary Fig. 4b). Decreased relative humidity (Supplementary Fig. 4c)

157

suggests that the reduced cloud cover during low BBA emission years is caused not only158

by reduced CCNs, but also through atmospheric feedbacks. We note that while emissions

159

of light-absorbing aerosols such as black carbon and brown carbon from fires can impose

160

a warming effect 30,31, it has been estimated by observation-based studies that the net

161

combined direct radiative forcing by all types of BBAs is negative 2,32,33 and that the net

162

forcing is also negative 34-36, which is consistent with our results.

163
164

The rectifier: nonlinear responses of clouds and soil water

165

The strong decadal warming over land in summer suggests that soil-related mechanisms

166

may provide the memory and nonlinearity to rectify the interannual aerosol forcing to the6

167

decadal warming. We can rule out a role of ocean for the summer subarctic warming

168

intensification, as the ocean surface rather absorbs heat from the atmosphere (increased

169

heat uptake) in summer in BBA_CMIP6 relative to BBA_smooth (Supplementary Fig.

1705), demonstrating that the decadal summer temperature increase over land is not caused

171

by heat transport from nearby oceans. The time-depth evolution of BBA_CMIP6 minus

172

BBA_smooth fields of soil ice, soil liquid, and their sum as net soil moisture over 50°N-

173

70°N (Fig. 2) reveals that thawing of upper soil permafrost, triggered by summers of low

174

aerosol emissions, propagates through the deeper soil layer over decadal timescales. The

175

increased liquid water in deeper soil indicates that water drains easily through more176

porous soil associated with permafrost thawing 37, contributing to the decadal drying

177

tendency of the upper soil layer. The decadal decrease in soil moisture and retreat of178

permafrost in response to interannual summertime BBA fluctuations indicates that179

thawing and freezing of soil are not symmetric between negative and positive BBA

180

conditions. Clearly, the upper permafrost thawing and soil moisture loss during low BBA

181

emission years exceed the upper permafrost freezing and soil moisture increase during

182

high emission years, thereby creating a hysteresis effect, which accumulates in time,

183

generating decadal-scale climate responses.

184
185

Several processes are involved in the nonlinear response of soil moisture to changes in

186

BBAs. Similar to the nonlinear relationship between anthropogenic aerosols and their187

indirect forcing 38, BBA forcing in our simulations exhibits stronger sensitivity under a188

cleaner environment (reduced aerosol conditions). Scatter plots of CCN versus BBA

189

reveal negative curvature, indicating that CCN formation is more sensitive to negative

190

BBA anomalies than to positive anomalies (Fig. 3a). The same argument holds for cloud191

droplet sensitivity to CCN (Fig. 3b) and to the inverse relationship between surface

192

shortwave fluxes and cloud droplet concentration (Fig. 3c). As a consequence of these193

nonlinearities, negative BBA anomalies can generate a stronger surface shortwave

194

response than positive BBA anomalies (Fig. 3d).

195
7

196

Nonlinearities also exist in other processes related to phase changes of water such as

197

sensitivity of soil ice melting/freezing to changes in radiation or temperature, although it

198

may be challenging to identify its isolated effect in our simulations as the phase change

199

processes involving water occur over the full range of surface to deeper soils, with

200

gradual heat and water transfer over multiple seasons (response time lag as a function of201

depth). The largest soil ice differences between the two ensemble groups are found at

202

subsurface depth of 0.9-1.9 m (Supplementary Fig. 6a). We have found a lagged

203

relationship between the net surface shortwave radiation to ice changes (tendency) at this204

subsurface depth of 0.9-1.9 m where the maximum ice differences are found. The

205

correlation between the surface shortwave flux and soil ice melting tendency at this layer

206

is maximized with a time lag of 2 months (correlation coefficient of 0.82). The scatter

207

diagram in Fig. 3e demonstrates that the subsurface soil ice melting is more sensitive to208higher surface shortwave fluxes than to lower shortwave fluxes. This is largely because

209

ice melts only at temperatures > 0°C, providing strong nonlinearity.

210
211

Albedo feedbacks involving snow cover changes can also add nonlinearity to the212

response 39. This, however, is not the main source of nonlinearity in summer as both213

ensemble groups become nearly snow-free over the most subarctic regions by July214

(Supplementary Fig. 7). Thus snow depth differences and albedo differences between the215

two ensemble groups are mostly not significant except for small areas adjacent to the216

Arctic (Supplementary Fig. 9a-b). The surface temperature response to shortwave fluxes

217

is almost linear on average over land between 50°N-70°N, although small snow covered

218

areas still exhibit a nonlinear behavior in surface temperature responses to shortwave flux

219

variations in summer (Supplementary Figs. 8c-d).

220
221

Overall, the combined effects of nonlinearities from aerosol-cloud-radiation interactions

222

and liquid-ice phase changes in soil result in a nonlinear response of soil ice to BBA

223

emissions (Fig. 3f). This promotes net thawing of permafrost on decadal scales in
8

224

BBA_CMIP6. As a consequence of permafrost thawing, an accelerated drainage of water225

gradually deprives the upper soil of moisture, providing the long-term climate memory

226(or reddening of the temporal spectrum).

227
228

The climate memory: decadal summer subarctic warming

229

The hydrological processes occurring in the upper soil layer have several types of230

impacts on surface temperature. First, the overall drying of the upper soil layer in

231

BBA_CMIP6 compared to BBA_smooth (Fig. 4b) lowers the heat capacity of the soil232

surface, inducing higher surface temperatures and associated positive feedbacks to near-

233

surface atmospheric warming in summer. The heat capacity of soil is determined by the234volumetric fraction of minerals, organic matters, and water content 40-42. Porous soil with235

low soil wetness is filled with air, which has lower heat capacity than liquid or ice water.

236

The lower heat capacity of the drier soil more readily facilitates surface warming in237

summer. Second, changes in evaporative fluxes modify surface temperature by absorbing

238

or releasing latent heat. The Bowen ratio, defined as the ratio of sensible to latent heat

239

flux, is a good indicator of this effect 43. For example, a high Bowen ratio indicates that240

the surface releases heat to the atmosphere more through sensible than latent heat fluxes

241

due to limited surface moisture to evaporate, raising surface temperature. Although not242

all areas over the subarctic land have the additional warming effect due to reduced243

evaporation, Supplementary Fig. 9h suggests that temperature damping effects by244

evaporation are suppressed over the regions where surface moisture is strongly depleted

245

in BBA_CMIP6 relative to BBA_smooth (drier soil moisture in Fig. 4b and higher

246

positive Bowen ratio differences Supplementary Fig. 9e). However, we note that, except

247

for those extremely dry areas, the evaporative effect does not seem to be the dominant

248

contributor to the overall decadal warming (surface temperatures are similar both over249

areas with positive Bowen ratio changes and with negative changes in Supplementary

250

Figs. 9f-g). This is likely because moisture can still be supplied by thawed surface soil

251

during summer, as inferred by the distributions of sensible and latent heat fluxes in

252

BBA_smooth and their change in BBA_CMIP6 (Supplementary Figs. 9a-d). Another9

253

impact of the hydrological processes on soil temperatures is through latent heat uptake

254

and release by melting and freezing of soil ice. The decreased upper soil ice content

255

causes less melting from spring to summer (Supplementary Fig. 10), absorbing less latent

256

heat of fusion, and thus sustaining less cooling (or more warming) at the surface. In other257

words, the presence of smaller amounts of ice through the transition from spring to258

summer contributes to additional surface heating in July.

259
260

The decadal warming in BBA_CMIP6 relative to BBA_smooth mediated by these soil

261

moisture processes is further enhanced by immediate cloud feedbacks. The atmospheric262

warming induced by soil interactions over the subarctic land domain is more intensified

263

towards the surface (Supplementary Fig. 11). Although specific humidity also increases

264

due to enhanced surface evaporation, the increase in atmospheric moisture is nevertheless265

modest relative to the temperature increase, thereby lowering relative humidity. This in266

turn reduces cloudiness throughout the mid- to lower-troposphere, with more cloud267

reduction occurring in the lower reaches. Fewer clouds allow more incoming solar

268

radiation to reach the surface, thereby heating the surface and lower atmosphere further

269(Fig. 4 and Supplementary Fig. 11). The presence of this atmospheric feedback is

270

supported by the positive shift in anomalous surface shortwave fluxes even for moderate

271

positive anomalous BBA emissions (see the dotted lines in Fig. 3d).

272
273

Prescribed soil moisture experiment274

The soil moisture feedback on summer surface climate is demonstrated by conducting

275

additional two sets of simulations with prescribed soil moisture calculated from CESM2-

276

LE. The purpose of this experiment is to quantify to what extent surface climate

277

differences in BBA_CMIP6 and BBA_smooth can be explained by soil moisture278

differences. The two simulation sets were all restarted from the year of 2000 and 2010279

from original 40 BBA_smooth members, and ran for a year under the same forcings

280

including BBA emissions. (The choice of the years is to utilize existing restart files.) The
10

281

only difference between the two sets is that for one set we prescribed daily mean soil

282

moisture fields (liquid water and ice) obtained from CESM2-LE climatology across the283

50 BBA_CMIP6 ensemble members and across 15 years (2000-2014), while for the other284

set from the 50 BBA_smooth members.

285
286

The surface temperature difference between the two sets in Fig. 5 demonstrates that287

reduced soil moisture induces substantial surface warming. As soil moisture is more

288

depleted at higher latitudes in BBA_CMIP6 compared to BBA_smooth (Fig. 4b), the land289

temperature response to soil moisture changes becomes larger closer to the Arctic (Fig.

290

5). The contribution of surface warming solely by soil drying (Fig. 5) relative to all291

contributions (Fig. 4a) is 22% over 50°N-60°N and 43% over 60°N-70°N. The soil

292

drying-induced warming is especially severe over Siberia (Fig. 5). This prescribed soil

293

moisture experiment confirms that the overall surface warming in BBA_CMIP6 relative294

to BBA_smooth in Fig. 4a is partly due to asymmetric response to reduced and enhanced

295

BBA fluctuations on interannual scales, and partly due to permafrost thawing-induced

296

soil moisture drying on decadal scales.

297
298

Discussion

299

Climate models are highly sensitive to a variety of aerosol forcings, including those from300

biomass burning. However, the sensitivity of the climate system to changes in the

301

interannual variance of aerosols has not previously been considered in experimental

302

designs for coordinated simulations such as CMIP6. Given that wildfires vary seasonally303

and from year-to-year, we have chosen to address this important question over the304

satellite observation period of fire emissions (1997-2014) using two groups of large

305

ensemble simulations with 50 members each. Our large ensemble approach using a single

306

Earth system model, comprised of 50 members forced by pronounced interannual BBA

307

variations during the satellite era following the CMIP6 protocol, and 50 members in308

which these interannual BBA variations are smoothed out 22, enables us to identify the
11

309

forced response signal that can otherwise be easily obscured by strong internal

310

variability. The 100 ensemble members for the historical period spanning 1850-1990

311

forced with identical BBA fluxes facilitate quantification of residuals from internal

312

variability that are not completely neutralized even after averaging 50 members

313(Methods), providing higher confidence in identifying a signal above natural variability

314

noise for the time interval of interest.

315
316

We found here that there is a strong decadal rectified response to interannual BBA

317

variability over the NH subarctic land regions, with a net summer warming in the318

presence of elevated interannual variability of biomass burning emissions. This mean

319

state response is mediated through interactions with subarctic soil moisture, resulting in

320

distinct summertime footprints over the NH high latitude land at the decadal timescale.

321

Although recent studies 7,8,23 have identified a similar warming response in the annual

322

mean to interannual fluctuations of BBA forcing, their proposed mechanistic pathways

323

for the climate response are less clear and distinct from the rectifier mechanisms

324

identified here. The earlier studies have hypothesized that asymptotic behaviors of325

aerosol effects on clouds with skewness in the emissions lead to a net reduction in low

326

cloud amount, which increases incoming shortwave radiation over northern high latitudes

327

7

328

as aerosol forcings along with the invoked atmospheric processes cannot sustain longer-329

term rectified responses (see the schematic timeseries of the nonlinearity-only response in

330

Fig. 6). Our study elucidates that indirect aerosol effects are only one of the triggers of331

decadal permafrost thawing, through which the decadal summer land warming is332

mediated (Fig. 6). With a simple conceptual framework for the combined effects of333

nonlinearity (by atmospheric processes and soil ice melting) and reddening (by soil

334

wetness serving as a climate memory), we can explain a decadal warming response

335

during the high BBA variability period of 1997-2014 and the recovery of the climate

336

system afterward (Fig. 6 and Supplementary Fig. 12).

. However, their proposed mechanism cannot explain the decadal scale climate response,

337
12

338

The rectified responses identified for surface variables such as temperature and soil

339

moisture are largely limited to the time interval of the perturbation (aerosol variance

340

modulations) itself, but longer-term responses do persist for the subsurface and deeper

341

soil layers in permafrost. The total loss of soil ice in BBA_CMIP6 is significant and its

342

impact lasts for decades (Supplementary Fig. 13). Depending on the season,

343

BBA_CMIP6 loses soil ice by a few percent up to ten percent relative to BBA_smooth344(Supplementary Fig. 14). In reality such a permafrost perturbation could potentially345

contribute to the release of methane, promoting additional warming. The permafrost

346

thawing-induced warming is also expected to modulate wildfire activity that can have347

feedback effects through ensuing carbon and aerosol releases. It is our hope that the348

rectifier mechanism identified in this study motivates improved representation of349

processes that impact BBA emissions, including through the inclusion of interactive fires

350

in models.

351
352

Our analyses have revealed that modulations of variance in aerosols and shortwave

353

radiation can be rectified by inherent nonlinearities and slow responses (reddening) in the354

soil system (Fig. 6). This new conceptual framework may have further implications for355

understanding the generation of long-term climate variability in response e.g. to year-to-356

year changes in fires, volcanic activity or even randomly-occurring interannual changes

357

in Arctic and subarctic cloudiness. The soil system serves as a memory for subarctic

358

summer climate in this study, but a similar conceptual model with Arctic sea ice and/or

359

ocean as a memory may also be applied to other components of the climate system. Our

360

results, which demonstrate mechanistically that spurious modulations of variability in

361

aerosol fluxes can lead to biases in the forced response, should be considered in the362

design of aerosol forcing protocols for the upcoming Coupled Model Intercomparison

363

Project Phase 7.

364
365
13

366

Methods367
368

CESM2-LE simulations

369

The 100 ensemble members comprising the CESM2-LE are divided into two subgroups

370

that are distinguished by their respective BBA emissions 22. The first group of 50

371

members follows the CMIP6 protocols for BBA and thereby includes strong interannual

372

variations during 1997-2014 due to the use of satellite-based Global Fire Emissions

373

Database (GFED) observations 3,6 (Supplementary Fig. 1b). The second group with 50

374

members has much smaller interannual variability, as the prescribed BBA fields are375

smoothed by applying 11-year running means at each grid point from the CMIP6 BBAs.

376

The 11-year running mean is performed for each month of the year separately to retain377

the annual cycle of BBA emissions. These two groups are referred to as BBA_CMIP6

378

and BBA_smooth, respectively, in this study. Net aerosol emissions are nearly conserved379

between the two groups of ensembles (difference < 0.35 % over 30°N-90°N and < 0.1%

380

globally for 1997-2014), and the aerosol difference only exists at the interannual scale.

381

Detailed descriptions of the dataset can be found in 22,27. To focus on responses over the382high latitude NH, the BBA emissions over 30°N-90°N is considered in this study. We

383

note that the strongest contribution is emitted over 50°N-70°N as shown in Fig. 1b, and384

results are not sensitive to the choice of a BBA emission latitude band between 30°N-

385

90°N and 50°N-70°N.

386
387

Estimating significance of composites for low – high emission summers

388

A bootstrap method is used for estimating uncertainty ranges for the composites of low

389

emission summers minus high emission summers in BBA_CMIP6 in Fig. 1 and390

Supplementary Fig. 4. A bootstrapped field is calculated by subtracting the mean of391

BBA_CMIP6 fields for 4 randomly selected years from the mean of BBA_CMIP6 fields

392

for 8 randomly selected years. Random years are selected between 1997-2014 without

393

replacement, and all calculations are for May to September (MJJAS). Values within one
14

394

standard deviation of 1000 bootstrapped cases are stippled in Figs. 1c-f and395

Supplementary Fig. 4 representing non-significant areas.

396
397

Estimating internal variability

398

Time series of differences in subarctic land temperatures between the first 50

399(BBA_CMIP6) and second 50 (BBA_smooth) members exhibit non-zero residuals even400

during the earlier historical period 1850-1990 (Supplementary Fig. 15). As identical401

aerosol fields are applied for all 100 members throughout this period, the non-zero402

differences between the two ensemble groups result from an incomplete cancellation of403

internal variability, even after averaging across 50 members. This suggests that the404

difference simulated during ~1990-2020 would not be entirely due to a climate response

405

to interannual variability of BBAs. The data availability of the same size of ensembles

406

over the longer historical period from 1850-1990 enables us to estimate a range of natural407

variability. As we are interested in identifying a climate response to BBA fluctuations by

408

comparing the full sets of 50 BBA_CMIP6 and 50 BBA_smooth members, we estimate a409

range of internal variability from differences between a random selection of 50 members

410

and the rest of the 50 members over 1880-1990. The earlier period 1850-1879 was

411

excluded to avoid potential impacts of initial conditions. Values within one standard

412

deviation of 1000 bootstrapped cases are stippled in Figs. 2b-d and 4, and Supplementary413

Figs. 1a, 2, 5, 8a-b, 9c-e, and 12 representing non-significant areas, and grey-shaded in

414

Supplementary Figs. 10 and 11 representing non-significant ranges.

415
416

Data availability

417

The CESM2-LE model output is available from418

https://www.cesm.ucar.edu/projects/community-projects/LENS2/data-sets.html.

419
420

Code availability
15

421

All codes used for analyses are developed by the corresponding author and available

422

upon request.

423
424

References

425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460

1

2
3

4
5
6
7
8
9
10
11
12

Fyfe, J. C., Kharin, V. V., Santer, B. D., Cole, J. N. S. & Gillett, N. P. Significant
impact of forcing uncertainty in a large ensemble of climate model simulations.
Proceedings of the National Academy of Sciences 118, e2016549118,
doi:10.1073/pnas.2016549118 (2021).
Brown, H. et al. Biomass burning aerosols in most climate models are too
absorbing. Nature Communications 12, 277, doi:10.1038/s41467-020-20482-9(2021).
van Marle, M. J. E. et al. Historic global biomass burning emissions for CMIP6(BB4CMIP) based on merging satellite observations with proxies and fire models(1750–2015). Geosci. Model Dev. 10, 3329-3357, doi:10.5194/gmd-10-33292017 (2017).
Lamarque, J. F. et al. Historical (1850–2000) gridded anthropogenic and biomassburning emissions of reactive gases and aerosols: methodology and application.
Atmos. Chem. Phys. 10, 7017-7039, doi:10.5194/acp-10-7017-2010 (2010).
Hoesly, R. M. et al. Historical (1750–2014) anthropogenic emissions of reactive
gases and aerosols from the Community Emissions Data System (CEDS). Geosci.
Model Dev. 11, 369-408, doi:10.5194/gmd-11-369-2018 (2018).
van der Werf, G. R. et al. Global fire emissions estimates during 1997–2016.
Earth Syst. Sci. Data 9, 697-720, doi:10.5194/essd-9-697-2017 (2017).
Fasullo, J. T. et al. Spurious Late Historical-Era Warming in CESM2 Driven by
Prescribed Biomass Burning Emissions. Geophysical Research Letters 49,
e2021GL097420, doi:https://doi.org/10.1029/2021GL097420 (2022).
DeRepentigny, P. et al. Enhanced simulated early 21st century Arctic sea ice loss
due to CMIP6 biomass burning emissions. Science Advances 8, 42,doi:10.1126/sciadv.abo2405 (2022).
Kim, J.-S., Kug, J.-S., Jeong, S.-J., Park, H. & Schaepman-Strub, G. Extensive
fires in southeastern Siberian permafrost linked to preceding Arctic Oscillation.
Science Advances 6, eaax3308, doi:doi:10.1126/sciadv.aax3308 (2020).
Shabbar, A., Skinner, W. & Flannigan, M. D. Prediction of Seasonal Forest Fire
Severity in Canada from Large-Scale Climate Patterns. Journal of Applied
Meteorology and Climatology 50, 785-799, doi:10.1175/2010jamc2547.1 (2011).
Balzter, H. et al. Impact of the Arctic Oscillation pattern on interannual forest firevariability in Central Siberia. Geophysical Research Letters 32,
doi:https://doi.org/10.1029/2005GL022526 (2005).
Macias Fauria, M. & Johnson, E. A. Large-scale climatic patterns control large
lightning fire occurrence in Canada and Alaska forest regions. Journal of
16

461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503

13
14

15

16

17

18
19
20
21
22
23

24
25

Geophysical Research: Biogeosciences 111,
doi:https://doi.org/10.1029/2006JG000181 (2006).
Bowman, D. M. J. S. et al. Fire in the Earth System. Science 324, 481-484,
doi:doi:10.1126/science.1163886 (2009).
Liu, Z., Yang, J., Chang, Y., Weisberg, P. J. & He, H. S. Spatial patterns anddrivers of fire occurrence and its future trend under climate change in a boreal
forest of Northeast China. Global Change Biology 18, 2041-2056,
doi:https://doi.org/10.1111/j.1365-2486.2012.02649.x (2012).
Flannigan, M. D., Krawchuk, M. A., de Groot, W. J., Wotton, B. M. & Gowman,
L. M. Implications of changing climate for global wildland fire. International
Journal of Wildland Fire 18, 483-507, doi:https://doi.org/10.1071/WF08187(2009).
Kasischke, E. S. & Turetsky, M. R. Recent changes in the fire regime across theNorth American boreal region—Spatial and temporal patterns of burning acrossCanada and Alaska. Geophysical Research Letters 33,
doi:https://doi.org/10.1029/2006GL025677 (2006).
Balshi, M. S. et al. Assessing the response of area burned to changing climate inwestern boreal North America using a Multivariate Adaptive Regression Splines(MARS) approach. Global Change Biology 15, 578-600,
doi:https://doi.org/10.1111/j.1365-2486.2008.01679.x (2009).
Podur, J., Martell, D. L. & Knight, K. Statistical quality control analysis of forest
fire activity in Canada. Canadian Journal of Forest Research 32, 195-205,
doi:10.1139/x01-183 (2002).
Stocks, B. J. et al. Large forest fires in Canada, 1959–1997. Journal ofGeophysical Research: Atmospheres 107, FFR 5-1-FFR 5-12,
doi:https://doi.org/10.1029/2001JD000484 (2002).
de Groot, W. J., Flannigan, M. D. & Cantin, A. S. Climate change impacts on
future boreal fire regimes. Forest Ecology and Management 294, 35-44,
doi:https://doi.org/10.1016/j.foreco.2012.09.027 (2013).
Witze, A. Why Arctic fires are bad news for climate change. Nature 585, 336-337(2020).
Rodgers, K. B. et al. Ubiquity of human-induced changes in climate variability.
Earth Syst. Dynam. 12, 1393–1411, doi:10.5194/esd-12-1393-2021 (2021).
Heyblom, K. B., Singh, H. A., Rasch, P. J. & DeRepentigny, P. Increased
Variability of Biomass Burning Emissions in CMIP6 Amplifies Hydrologic Cycle
in the CESM2 Large Ensemble. Geophysical Research Letters 49,
e2021GL096868, doi:https://doi.org/10.1029/2021GL096868 (2022).
Chung, E.-S. et al. Cold-Season Arctic Amplification Driven by Arctic OceanMediated Seasonal Energy Transfer. Earth's Future 9, e2020EF001898,
doi:https://doi.org/10.1029/2020EF001898 (2021).
Sejas, S. A. et al. Individual Feedback Contributions to the Seasonality of Surface
Warming. Journal of Climate 27, 5653-5669, doi:10.1175/jcli-d-13-00658.1(2014).

17

504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548

26
27
28
29
30
31
32
33
34
35
36
37
38
39
40

41

Deser, C., Tomas, R., Alexander, M. & Lawrence, D. The Seasonal Atmospheric
Response to Projected Arctic Sea Ice Loss in the Late Twenty-First Century.
Journal of Climate 23, 333-351, doi:10.1175/2009jcli3053.1 (2010).
Danabasoglu, G. et al. The Community Earth System Model Version 2 (CESM2).
Journal of Advances in Modeling Earth Systems 12, e2019MS001916,
doi:https://doi.org/10.1029/2019MS001916 (2020).
Serreze, M. C. & Francis, J. A. The Arctic Amplification Debate. Climatic
Change 76, 241-264, doi:10.1007/s10584-005-9017-y (2006).
Screen, J. A. & Simmonds, I. The central role of diminishing sea ice in recent
Arctic temperature amplification. Nature 464, 1334-1337,
doi:10.1038/nature09051 (2010).
Bond, T. C. et al. Bounding the role of black carbon in the climate system: A
scientific assessment. Journal of Geophysical Research: Atmospheres 118, 53805552, doi:https://doi.org/10.1002/jgrd.50171 (2013).
Feng, Y., Ramanathan, V. & Kotamarthi, V. R. Brown carbon: a significant
atmospheric absorber of solar radiation? Atmos. Chem. Phys. 13, 8607-8621,
doi:10.5194/acp-13-8607-2013 (2013).
Chang, D. Y. et al. Direct radiative forcing of biomass burning aerosols from theextensive Australian wildfires in 2019–2020. Environmental Research Letters 16,
044041, doi:10.1088/1748-9326/abecfe (2021).
Shi, S. et al. Biomass burning aerosol characteristics for different vegetation types
in different aging periods. Environment International 126, 504-511,
doi:https://doi.org/10.1016/j.envint.2019.02.073 (2019).
Penner, J. E., Dickinson, R. E. & O'Neill, C. A. Effects of Aerosol from Biomass
Burning on the Global Radiation Budget. Science 256, 1432-1434,
doi:doi:10.1126/science.256.5062.1432 (1992).
Ward, D. S. et al. The changing radiative forcing of fires: global model estimatesfor past, present and future. Atmos. Chem. Phys. 12, 10857-10886,
doi:10.5194/acp-12-10857-2012 (2012).
Johnson, B. T. et al. Evaluation of biomass burning aerosols in the HadGEM3
climate model with observations from the SAMBBA field campaign. Atmos.
Chem. Phys. 16, 14657-14685, doi:10.5194/acp-16-14657-2016 (2016).
Liljedahl, A. K. et al. Pan-Arctic ice-wedge degradation in warming permafrost
and its influence on tundra hydrology. Nature Geoscience 9, 312-318,
doi:10.1038/ngeo2674 (2016).
Carslaw, K. S. et al. Large contribution of natural aerosols to uncertainty inindirect forcing. Nature 503, 67-71, doi:10.1038/nature12674 (2013).
Alexander, M. A., Tomas, R., Deser, C. & Lawrence, D. M. The Atmospheric
Response to Projected Terrestrial Snow Changes in the Late Twenty-First
Century. Journal of Climate 23, 6430-6437, doi:10.1175/2010jcli3899.1 (2010).
Lawrence, D. M. et al. The Community Land Model Version 5: Description of
New Features, Benchmarking, and Impact of Forcing Uncertainty. Journal ofAdvances in Modeling Earth Systems 11, 4245-4287,
doi:https://doi.org/10.1029/2018MS001583 (2019).
deVries, D. A. Thermal Properties of Soils. (North-Holland, 1963).
18

549
550
551

42
43

Hillel, D. Environmental soil physics. (Academic press, 1998).
Bowen, I. S. The ratio of heat losses by conduction and by evaporation from any
water surface. Physical Review 27, 779-787 (1926).

552
553

Acknowledgments

554

The CESM2 Large Ensemble (CESM2-LE) project was conducted through a partnership

555

between the IBS Center for Climate Physics (ICCP) in South Korea and the Community

556

Earth System Model (CESM) group at the National Center for Atmospheric Research

557(NCAR) in the US. The authors would like to thank all of the scientists, software

558

engineers, and administrators at both NCAR and the ICCP that contributed to the project.

559

The simulations presented here were carried out on the IBS/ICCP supercomputer

560“Aleph”, a 1.43 petaflop high-performance Cray XC50-LE Skylake computing system

561

with 18720 processor cores, 9.59 PB of disc storage, and 43 PB of tape archive storage.

562

The work of J.-E.K., R.Y., K.B.R., A.T., S.-S.L, and K.S. was supported by the Institute563

for Basic Sciences (IBS), Republic of Korea, under IBS-R028-D1. The CESM project is564

supported primarily by the US National Science Foundation (NSF). NCAR is a major

565

facility sponsored by the US NSF under cooperative agreement 1852977. M.F.S. was

566

supported by NOAA’s Climate Program Office’s Modeling, Analysis, Predictions, and567

Projections (MAPP) program grant NA20OAR4310445. This is IPRC publication X and

568

SOEST contribution Y.

569
570

Author information

571

Authors and Affiliations

572

Center for Climate Physics, Institute for Basic Science, Busan, South Korea.

573

Pusan National University, Busan, South Korea

574

Ji-Eun Kim, Keith B. Rodgers, Axel Timmermann, Sun-Seon Lee & Karl Stein

575
19

576

Research Institute for Global Change, Japan Agency of Marine-Earth Science and577

Technology, Kanagawa, Japan

578

Ryohei Yamaguchi

579
580

National Center for Atmospheric Research, Boulder, CO, USA581

Gokhan Danabasoglu, Jean-Francois Lamarque, John T. Fasullo, Clara Deser, Nan

582

Rosenbloom & Jim Edwards

583
584

Department of Oceanography and International Pacific Research Center, School of585

Ocean and Earth Science and Technology, University of Hawaiʻi at Mānoa, Honolulu,

586

HI, USA

587

Malte F. Stuecker

588
589
590

Contributions

591

J.-E.K., K.B.R., A.T., R. Y., G.D., C.D. and J.-F.L. designed the smoothed biomass

592

burning emissions for the 50 ensemble members of CESM2-LE. J.-E.K., R.Y., K.B.R.,

593

A.T., S.-S.L., and K.S. initiated the work. J.-E.K. and R.Y. performed scientific analysis

594

and J.-E.K., K.B.R. and A.T. wrote the manuscript. The CESM2-LE simulations were595

conducted by S.-S.L with the support of N.R. and J.E. J.-E.K. performed the prescribed596

soil moisture experiment. All authors discussed the results and reviewed the manuscript.

597
20

598

Corresponding author

599

Correspondence to Ji-Eun Kim

600
601

Ethics declarations

602

Competing Interest Statement: The authors declare no competing interests.

603
604

Supplementary information

605

Supplementary information file

606
607
608
609
610
611
612
613
614
615
616
617
618
21

619

Figures

620

621
622

Fig. 1. Prescribed BBA difference and atmospheric responses to reduced BBA

623

emissions. (a) Monthly time series of prescribed BBA emissions anomalies in

624

BBA_CMIP6 relative to BBA_smooth, integrated over 30°N-90°N. Dot colors

625

correspond to month labeled across the top. (b) The standard deviation of annual BBA

626

differences between BBA_CMIP6 and BBA_smooth over 1997-2014. Composites of 8627

low BBA emission years (1997, 1999, 2000, 2001, 2005, 2007, 2009, 2011) minus

628

composites of 4 high BBA emission years (1998, 2002, 2003, 2012) in BBA_CMIP6 for629(c) aerosol optical depth (AOD), (d) cloud condensation nuclei (CCN), (e) the net

630

shortwave flux at the surface (positive downward), and (f) surface temperature. The 8631

low and 4 high emission years in BBA_CMIP6 are selected based on the BBA emission632

differences, shown in (a), for summer months for May to September (MJJAS). The

633

composite difference maps in (c-f) are for the same summer months (MJJAS). Stippling

634

denotes non-significant areas (see Methods).

635
22

636
637

Fig. 2. Evolution of soil ice and surface water and temperature in response to BBA

638

fluctuations over 50°N-70°N land. (a) Monthly time series of BBA flux difference over639

30°N-90°N (same as Fig. 1a). Monthly differences of (b) soil ice, (c) soil liquid, and (d)640

net soil water (ice + liquid) show that melting of ice propagates into deeper soil and that641

upper layer soil moisture is diminished in BBA_CMIP6. (e) Monthly time series of642

differences of surface temperature and (f) upper 10 cm soil moisture show decadal scale

643

changes. Dot colors in time series represent months (same as in Fig. 1a) with enlarged

644pink dots for July. Unit is kg m-2 for (b-d) and (f). Stippling in (b-d) denotes non-645

significant areas whose values are within one standard deviation of an internal variability

646

range (see Methods). All values correspond to differences between the two ensemble

647

groups (BBA_CMIP6 – BBA_smooth).
23

648
649

Fig. 3. Nonlinear sensitivity of BBA-cloud-radiation and BBA-soil ice. Relationship

650

between anomalies of (a) prescribed BBA flux and CCN, (b) CCN and cloud droplet

651

concentration, (c) cloud droplet concentration and net surface shortwave flux by clouds652(positive downward), (d) BBA and net surface shortwave flux, (e) net surface shortwave

653

flux and subsurface (0.9-1.9 m depth) soil ice tendency, and (f) BBA flux and total soil654

ice tendency. All values are taken from BBA_CMIP6 – BBA_smooth over 1997 – 2014.

655

CCN and cloud droplet concentrations are vertically integrated throughout the656

troposphere. The BBA flux is integrated over 30°N-90°N, and all other simulated values

657

are averaged over the land domain between 50°N-70°N. In (a-f), each scatter dot

658

corresponds to a 1-month value for June (blue), July (orange), and August (green). In (e),659

the ice tendency is calculated as a change in subsurface ice during 2 months following

660

surface shortwave flux for a given month. In (f), each scatter dot corresponds to June-

661

July-August (JJA) BBA flux versus a soil ice change from a previous year. Since the soil662

ice response to summer BBA emissions propagates to deeper soil over following months663

as depicted in Fig. 2b, mean soil ice for summer to the following winter (June to next
24

664

February) was taken for each given year (e.g., the dot for 1998 corresponds to the BBA

665

flux in JJA 1998 versus the mean soil ice for June 1998 – February 1999 minus the mean

666

soil ice for June 1997 – February 1998). Spring values are excluded in the soil ice

667

calculation given the ambiguity as to whether spring soil ice changes are influenced by668

the previous summer or the same year’s spring emissions (Note that spring BBA

669

emissions are irregular as shown in Fig. 1a). The black curve is a least-square second-670

degree polynomial fit from filled dots. Unfilled dots are considered as outliers whose671

squared error is greater than the three-standard deviation of all squared errors.

672
673
674
675
676
677

678
679

Fig. 4. July responses in BBA_CMIP6 – BBA_smooth. Difference maps of (a)680

temperature, (b) upper 10 cm soil moisture, (c) net surface shortwave flux (positive681

downward), (d) low-level (> 700 hPa) relative humidity, and (e) low-level cloud fraction682

in July over 2000-2014. Stippling denotes non-significant areas whose values are within

683

one standard deviation of an internal variability range (see Methods).

684
25

685
686

687
688

Fig. 5. Surface temperature response from the prescribed soil moisture experiment689

using CESM2. Color shows surface temperature difference in July between the two

690

prescribed soil moisture cases of BBA_CMIP6 and BBA_smooth. Except for prescribing

691

different daily soil liquid water and ice fields, all other simulation conditions have692

remained the same for the two simulation sets. See simulation details in Results.

693
694

26

695
696

Fig. 6. Schematic of atmosphere-land-permafrost coupling and generation of697

decadal variability from interannual forcing. For each given year, reduced (enhanced)

698

aerosol emissions from biomass burning induces net warming (cooling) of the699

atmosphere and surface through direct and indirect aerosol impacts. As a consequence of

700

nonlinear sensitivities in aerosol-cloud-radiation processes and soil ice melting,701

interannual changes in aerosols result in a decadal scale permafrost loss, promoting

702

drainage of upper soil water to the deeper soil layer. The decadal decrease in upper soil

703

moisture by the melt-water drainage induces decadal warming of the surface and lower704

troposphere that is intensified by cloud feedbacks through changes in relative humidity.

705

The bottom panels illustrate three conceptual response models when (upper left)

706

symmetric interannual fluctuation of aerosols is present, mimicking summertime aerosol
27

707

emission pulses. (lower left) A response varies at the same interannual timescale as the708

emission timescale when there is only nonlinearity in atmospheric processes. (upper709

right) A response fluctuates interannually with a gradual decay of a signal (reddening)

710

after a symmetric linear response to an aerosol pulse when a climate system has a711

memory effect. The decay timescale of 2 years is used in this example. (lower right) A

712

response is rectified, having longer timescales than the emission timescale, when both713

nonlinearity and reddening are present. Units are arbitrary. See Supplementary Fig. 12 for714

conceptual model responses to aerosol fluctuations of BBA_CMIP6 – BBA_smooth.

28