ARTICLE https://doi.org/10.1038/s41467-022-29013-0 OPEN Wildfires enhance phytoplankton production in tropical oceans Dongyan Liu 1 ✉, Chongran Zhou 1, John K. Keesing 2 ✉, Oscar Serrano 3,4, Axel Werner 3, Yin Fang 5, Yingjun Chen 6 ✉, Pere Masque 3,7,8, Janine Kinloch 9, Aleksey Sadekov10 & Yan Du 11 1234567890():,; Wildfire magnitude and frequency have greatly escalated on a global scale. Wildfire products rich in biogenic elements can enter the ocean through atmospheric and river inputs, but their contribution to marine phytoplankton production is poorly understood. Here, using geo- chemical paleo-reconstructions, a century-long relationship between wildfire magnitude and marine phytoplankton production is established in a fire-prone region of Kimberley coast, Australia. A positive correlation is identified between wildfire and phytoplankton production on a decadal scale. The importance of wildfire on marine phytoplankton production is sta- tistically higher than that of tropical cyclones and rainfall, when strong El Niño Southern Oscillation coincides with the positive phase of Indian Ocean Dipole. Interdecadal chlorophyll-a variation along the Kimberley coast validates the spatial connection of this phenomenon. Findings from this study suggest that the role of additional nutrients from wildfires has to be considered when projecting impacts of global warming on marine phy- toplankton production. 1 State Key Laboratory of Estuarine and Coastal Research, Institute of Eco-Chongming, East China Normal University, Shanghai 200062, China. 2 CSIRO Oceans and Atmosphere Research, and University of Western Australia Oceans Institute, Indian Ocean Marine Research Centre, Crawley, WA, Australia. 3 School of Science and Centre for Marine Ecosystems Research, Edith Cowan University, Joondalup, WA, Australia. 4 Centro de Estudios Avanzados de Blanes, Consejo Superior de Investigaciones Científicas, Blanes, Spain. 5 College of Marine Ecology and Environment, Shanghai Ocean University, Shanghai 201306, China. 6 Shanghai Key Laboratory of Atmospheric Particle Pollution and Prevention, Department of Environmental Science and Engineering, Fudan University, Shanghai 200438, China. 7 Departament de Física and Institut de Ciència i Tecnologia Ambientals, Universitat Autònoma de Barcelona, Bellaterra, Spain. 8 International Atomic Energy Agency, 4a Quai Antoine 1er, 98000 Principality of Monaco, Monaco. 9 Biodiversity and Conservation Science, Department of Biodiversity, Conservation and Attractions, Bentley Delivery Centre, Bentley, WA, Australia. 10 Ocean Graduate School, ARC Centre of Excellence for Coral Reef Studies, University of Western Australia, Crawley, WA 6009, Australia. 11 State Key Laboratory of Tropical Oceanography, South China Sea Institute of Oceanology, Chinese Academy of Sciences, Guangzhou, China. ✉email: dyliu@sklec.ecnu.edu.cn; john.keesing@csiro.au; yjchenfd@fudan.edu.cn NATURE COMMUNICATIONS | (2022)13:1348 | https://doi.org/10.1038/s41467-022-29013-0 | www.nature.com/naturecommunications 1 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-29013-0 M arine phytoplanktons are critical primary producers, approximately 26.5 million tons enter the ocean every year17. contributing nearly half of the biosphere’s net primary Despite their importance, our understanding of the effect of production1. The impact of global warming on marine wildfires on the ocean is far less than our understanding of their phytoplankton production (MPP) has become increasingly evi- role in terrestrial ecosystems. dent in the past few decades2, but how MPP is affected at dif- Northern Australia is one of the most fire-prone savanna ferent latitudes remains controversial3–5. For example, regions of the world (Fig. 1a). As with most of Australia, hot assimilating ocean colour satellite data into a marine biogeo- weather associated with the El Niño Southern Oscillation (ENSO) chemical model showed that global net MPP experienced a small is the main driver of wildfires18; the risk and severity of wildfires yet significant decline of −0.8 Pg C yr−1 or −2.1% decade−1 from can increase when the Indian Ocean Dipole (IOD) is also in its 1998 to 20154. Ocean stratification enhanced by warming has positive phase (pIOD). The drying effect of the easterly shift in been invoked as a major mechanism to explain the decline in equatorial trade winds induced by pIOD can promote fire con- MPP: warming modifies the mixed layer depth and reduces the ditions and lead to extra-strong wildfires19. For instance, during vertical mixing of the surface layer and underlying cooler the 2019–2020 bushfires, a globally unprecedented 21% of the nutrient-rich waters below the permanent pycnocline6. Changes Australian temperate and broadleaf mixed forest biome in these physical processes lead to a reduced supply of nutrients burned18,19. Hypothetical conjecture, when strong ENSO occurs to the upper ocean and, consequently, have a negative impact on during the pIOD phase, escalating wildfires can increase the flux MPP; for example, large-sized phytoplankton species, such as of biogenic elements into the ocean via the pathways of atmo- diatoms, can significantly decrease in abundance and cause MPP spheric deposition and riverine input and lead to higher MPP. To decline in the upper ocean7,8. However, although this might be prove their connection in the context of climate modes, it is true in temperate oceans, it is not always the case in the tropics necessary to examine a decadal relationship between wildfire and Polar Regions9,10. Generally, the rise in temperature and magnitude and MPP. changes in stratification in tropical oceans are smaller than in Palaeocological methods can provide an effective pathway to temperate oceans because of their larger heat capacity and ther- reconstruct the interdecadal correlations between environmental mal inertia9. Moreover, warming-induced secondary climate and biological information. Total organic carbon (TOC), total effects, such as increased tropical cyclones and upwellings, can nitrogen (TN), and biogenic silicate (BSi) preserved in sediment compensate for nutrient depletion in the upper ocean by accel- cores are important geochemical proxies to reconstruct MPP in erating turbulent mixing and promoting rainfall11. Phytoplankton the upper ocean. TOC and TN reflect the accumulation of organic metabolic capacity in the tropics is also higher than that at high matter in the seafloor, while BSi is a frustule component of dia- latitudes12. These factors can partially offset the negative effect of toms and is often used as a proxy for diatom biomass in the ocean greater stratification and even lead to a stage increase in MPP in because of its major contribution to MPP20,21. The effectiveness tropical oceans5,11. of BSi as an indicator of MPP in northern Australia has been In contrast to the impact of physical processes on MPP in validated using biomarkers22. Wildfire magnitude can be recon- tropical oceans, the role of wildfires has received minimal structed using black carbon (BC) contents preserved in the attention. The risk and severity of wildfires in the southern sediments. BC is an organic, molecularly diverse product result- hemisphere have greatly escalated on a global scale as a con- ing from the incomplete combustion of biomass and fossil fuels23 sequence of rising temperatures and more frequent heat and it decomposes slowly after burial in marine sediments as a waves13,14; for example, there was a distinct lengthening of the component of TOC17. Moreover, the source of BC from biomass fire weather season between 1980 and 2018 (Fig. 1a). The emis- or fossil fuels can be distinguished via the ratios char/soot. In sions and ash from wildfire are rich in biogenic elements15, such high-temperature fossil fuel combustion (e.g. vehicle emissions as nitrogen, phosphate, silicate, and iron, and exert a distinct and industrial coal combustion) the ratio of char/soot ratio is less impact on atmospheric and aquatic environments16,17. For than one, while for the relatively low-temperature biomass example, particles emitted by wildfires account for approximately burning (e.g. wildfire), the ratio of char/soot is much higher than 62% of the global annual emissions of organic matter from bio- one24. mass burning16. The global flux of soluble charcoal from biomass Therefore, we measured these geochemical parameters in three burning is estimated to be 40–250 million tons yr−1, and sediment cores (ID: 185, 200, and KGR) collected from the a b Kimberley Change in the length of the fire weather season from 1980 to 2018 (Days) Fig. 1 Maps showing changes in the length of the fire weather season on a global scale and sampling location. a Map showing the change in the length of the fire weather season between 1980 and 2018 (data from the ERA5 data set13,14; the green circle indicates the Kimberley coast, a fire-prone region in northern Australia). b Sampling sites of core 185, core 200, and core KGR (red dots) in the Kimberley coast. 2 NATURE COMMUNICATIONS | (2022)13:1348 | https://doi.org/10.1038/s41467-022-29013-0 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-29013-0 ARTICLE Kimberley coast of northern Australia (Fig. 1b), aiming to in Core 200 showed a distinct increase during the 1990s, reconstruct a 100-year time-series relationship between wildfire corresponding to relatively higher BC contents. However, these magnitude and MPP variation. Kimberley region is one of differences in concentrations did not affect the positive correla- the most fire-prone savanna regions of the world (Fig. 1a). The tion between BSi and BC during pIOD phase. These results are Kimberley coastline and its river catchments are sparsely popu- consistent with the findings in northern Australia, where there lated and remote, with relatively low levels of fuel consumption25. was a 1.5–3 times increase in primary production between The ratios char/soot in the three cores ranged 6.9 to 22.3 for core the 1960s and 2010s, mainly contributed by diatoms and high 185, from 9.8 to 26.2 for core 200, and from 12.3 to 31.7 for core primary production occurred during the ENSO years22. KGR, respectively (Supplementary Fig. 1), these numbers indicate that the BC in marine sediments of Kimberley coast pre- Observational data to verify palaeo-reconstruction. To verify dominantly originates from wildfires. Thus, it is reasonable to the accuracy of palaeo-reconstruction using geochemical meth- reconstruct a chronological relationship between wildfire mag- ods, the relationship between BC and fire data in the catchment nitude and MPP, using BC and BSi, respectively. The correlation from 1988 to 2018 was analysed. An archive of burnt areas, between BC and BSi and their characteristics over time are ana- mapped using NOAA imagery26, was used to examine fire extent lysed in the context of IOD and ENSO climate modes (Fig. 2a), between 1988 and 2018 in the catchments (Fig. 3a–c), and these referring to the fire record of the National Oceanic and Atmo- can be usefully divided into the less intense early dry season (EDS: spheric Administration (NOAA) during 1988–2018 and the January–June) fires and the more intense and potentially more variation of satellite chlorophyll-a (Chl-a) of Kimberley coast destructive late dry season (LDS: July–December) fires27. EDS from 2003 to 2018. The importance of wildfire effect on phyto- fires include managed burn regimes by government and local plankton production is statistically compared with that of tropical indigenous authorities to reduce ground fuel loads intended to cyclones and rainfall during pIOD and negative IOD (nIOD) prevent more severe fires later in the year. EDS-managed burns phases, respectively. were implemented in 2011 and are attributed with reducing the LDS wildfires by about half although the total area burnt on average remained the same28. This can be reflected by the relative Results changes of EDS and LDS areas before and after 2011 (Fig. 3a–c): The positive effect of wildfire on MPP. The period between 1920 increased EDS fires and decreased LDS fires. and 2017 can be divided into three climate modes, including EDS lower temperature managed burning results in less nIOD phase (1920–1960), a fluctuation between nIOD and pIOD complete combustion and thus produces more BC. Burnt areas phases (1960–2005), and pIOD phase (2005–2018) (Fig. 2a). The during the EDS showed an obvious increase during 2011–2017, three sedimentary cores (185, 200, and KGR) cover the time span corresponding to increased BC in the three cores after 2010 from the 1920s to 2010s, representing two phases of pIOD (Fig. 2c, d, g). Moreover, during the EDS, high rainfall during dominance (1991–2017) and nIOD dominance (1926–1990). January to March and southward, offshore winds from March to There are two common characteristics in the three cores: (1) BSi June would readily transport the emission and ash of wildfires and BC displayed a strongly positive correlation in the three cores from the land to the sea (Fig. 3d). In contrast, northerly, onshore during the phase of pIOD dominance (Table 1). (2) The fre- winds and dry weather dominate during the LDS, which can quency and magnitude of RSI increased evidently in the three facilitate fire products to remain inland and reduce their cores during the phase of pIOD dominance (Fig. 2b–g); a com- contribution to the sea. This helps to explain the low BC mon increase for TOC, TN, BSi, and BC appeared after 2010s contents in the sediments at the peak of the LDS fire around 2005 when strong ENSO conditions coincide with pIOD phase (Fig. 2). (Figs. 2 and 3). These findings demonstrate the effectiveness of These positive signals indicate an important linkage between BC as a proxy of wildfire frequency and/or intensity over this wildfire and MPP during pIOD phase. period and suggest that the impact of EDS fires on the ocean may During the phase of nIOD dominance (1926–1990), no positive be more important than that of LDS fires. correlation was found between TOC, TN, BSi, and BC in Core In order to assess the broader spatial effect of this 185 and KGR, but BC displayed positive correlations with TOC phenomenon, we further analysed the relationships between (r = 0.56, p < 0.01) and TN (r = 0.56, p < 0.05) in Core 200 satellite Chl-a concentrations from 2003 to 2017 in the Kimberley (Table 1). During the phase of pIOD dominance (1991–2017), coast (total sites: 94,190) and the contents of BC and BSi in Core BSi and BC displayed strongly positive correlations with TOC 200 (2003–2017) (Fig. 4). Chl-a was positively correlated with BC and TN in Cores 200, a positive correlation between BC and TOC at 34.4% of the sites (Fig. 4a) and with BSi at 54.0% of the sites (r = 0.67, p < 0.05) was found in Core KGR (Table 1). In (Fig. 4b). These results suggest that wildfires may have an comparison, the positive contribution of BSi and BC to organic extensive spatial effect on MPP in the Kimberley coast. matter greatly increased during the phase of pIOD dominance. After the 1990s, the contents of BSi and BC in the three cores changed in a V-shape (Fig. 2c, e, g). These two peaks are basically Discussion consistent with the time when strong ENSO (2003, 2015) The effects of increasing wildfires on biogeochemical cycling and overlapped with pIOD (Fig. 2a). The synchronous variations primary production in terrestrial ecosystems have attracted great between BSi and BC in the pIOD phase indicate the impact of attention29, but are barely taken into account in marine ecosys- wildfires. tems. The results from geochemical reconstruction and observa- The contents of TOC, TN, BSi, and BC in the three cores are tional data in the Kimberley coast indicated that wildfire plays a different. In comparison, they were higher in Core185 than in role in MPP enhancement during the pIOD phase. North-western Core 200 and KGR (Fig. 2b–g). This might be related to the Australia is one of the world’s tropical cyclone hot spots30. difference of grain sizes in the three cores and the variations of Research investigating MPP increases in northern Australia has freshwater discharge from the three rivers or differences in mainly focussed on the contribution of intensified tropical catchment size or characteristics such as size soil and vegetation cyclones and rainfall forced by rising temperatures during ENSO type (Supplementary Fig. 2), e.g. the median grain size (d50) in years22,31,32. Modelling results have suggested that tropical Core 185 showed a distinct decrease during 1998–2018, cyclone‐induced phytoplankton blooms in north-western Aus- corresponding to relatively lower BC (Fig. 2c); in contrast, d50 tralia could contribute to 20% of annual primary production11. NATURE COMMUNICATIONS | (2022)13:1348 | https://doi.org/10.1038/s41467-022-29013-0 | www.nature.com/naturecommunications 3 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-29013-0 IOD ENSO a 0.6 1.0 Nino 3.4 DMI 0.0 0.0 -0.6 -1.0 1920 1930 1940 1950 1960 1970 1980 1990 2000 2010 2020 b 2.4 Core 185 1.41 0.21 0.64 0.55 0.65 TOC (%) 2.1 TN (%) 0.14 1.8 0.86 0.27 0.78 1.5 0.07 1920 1930 1940 1950 1960 1970 1980 1990 2000 2010 2020 c 4.0 Core 185 0.55 -0.77 BSi (%) 0.45 BC (%) 3.0 0.35 -0.67 0.02 -0.40 2.0 0.25 1920 1930 1940 1950 1960 1970 1980 1990 2000 2010 2020 d 1.5 Core 200 1.20 0.25 0.23 TOC (%) 1.2 TN (%) 0.20 0.15 0.9 0.86 0.60 0.21 0.6 0.05 1920 1930 1940 1950 1960 1970 1980 1990 2000 2010 2020 e 2.0 Core 200 0.05 0.45 -0.43 0.17 BSi (%) 0.35 BC (%) 1.5 0.50 0.25 0.36 1.0 0.15 1920 1930 1940 1950 1960 1970 1980 1990 2000 2010 2020 f 1.0 Core KGR 0.11 TOC (%) 0.9 TN (%) 0.08 0.8 0.24 -1.18 0.7 0.05 1920 1930 1940 1950 1960 1970 1980 1990 2000 2010 2020 g 1.5 Core KGR 0.13 BSi (%) 1.3 BC (%) -0.56 0.50 0.07 1.1 0.02 0.9 0.01 1920 1930 1940 1950 1960 1970 1980 1990 2000 2010 2020 Year Fig. 2 The profiles to illustrate the variations of El Niño Southern Oscillation (ENSO), the Indian Ocean Dipole (IOD), and multiple geochemical proxies in the three sediment cores from the 1920s to 2010s. a Low-frequency signals of dipole mode index (DMI) and Niño 3.4 index showing the variations of IOD and ENSO, respectively. b Core 185: total organic carbon (TOC) and total nitrogen (TN). c Core 185: biogenic silicate (BSi) and black carbon (BC). d Core 200: TOC and TN. e Core 200: BSi and BC. f Core KGR: TOC and TN. g Core KGR: BSi and BC. Blue lines represent the shift changes assessed by sequential t test analysis of regime shift and numbers on lines were regime shift index (RSI) to show shifting magnitude. 4 NATURE COMMUNICATIONS | (2022)13:1348 | https://doi.org/10.1038/s41467-022-29013-0 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-29013-0 ARTICLE The climate data from the present study indicated that an nIOD (n = 7) >0.05/>0.05 >0.05/>0.05 0.00/0.04 Table 1 Pearson correlations (two sided) between each pair of variables in the three cores during the entire time series (ETS: 1926–2017) and the phases of positive Indian Ocean Dipole (pIOD) dominance (1991–2017) and negative Indian Ocean Dipole (nIOD) dominance (1926–1990), respectively (n: number for statistical analysis; r: correlation 0.56/0.00 coefficient; p: the level of significance; significant: p < 0.05; highly significant: p < 0.01; no significant: p > 0.05; bold numbers indicate a significant correlation between two increase in wildfire, sea surface temperature (SST), tropical —/−0.46 —/>0.05 cyclone frequency, and rainfall are either synchronised or lagged during the pIOD phase (Fig. 2 and Supplementary Fig. 3). Therefore, it is necessary to assess their relative importance in explaining the variance in BSi over time. Linear Fixed Effect Models and multiple linear regression were pIOD (n = 10) >0.05/<0.05 >0.05/>0.05 applied for the data of BC, SST, TCF, and rainfall in explaining 0.24/0.67 0.13/0.48 the annual variance of BSi during ETS, pIOD, and nIOD, —/<0.05 —/0.71 respectively (Table 2). Model results from cores 185 and 200 showed the positive significance of BC for BSi during ETS; in contrast, the significance of TCF and rainfall for BSi was lower Core KGR: BSi/BC than BC and only appear in Core 200 (Table 2). Multiple linear regression for cores 185 and 200 further showed that BC can >0.05/<0.05 >0.05/>0.05 ETS (n = 17) explain 36.8 and 8.2% of the variance in annual BSi in the nIOD 0.16/ 0.38 0.10/0.55 —/>0.05 phase, respectively, but its importance increased to 47.6 and 62% —/0.47 in the pIOD phase, respectively (Table 2). In contrast, SST, tro- pical cyclones, and rainfall can explain more about the variance in annual BSi in the nIOD phase (Table 2). The result is consistent with the model evaluations in the eastern Indian Ocean between nIOD (n = 26) 1997 and 2009, which showed that the atmospheric deposition of >0.05/<0.05 >0.05/<0.01 −0.05/0.56 −0.13/0.56 nutrients from wildfire accounted for a higher proportion of MPP —/−0.10 —/>0.05 increase than that of wind-driven nutrient upwelling in pIOD phases33. Although the importance of wildfires to MPP in the pIOD phase far exceeds our expectations, the phenomenon is reason- able to some extent. Atmospheric and river transport of macro- pIOD (n = 28) nutrients (e.g. nitrogen) and bio-essential metals (e.g. iron) <0.05/<0.01 <0.01/<0.01 produced by wildfire emission have been considered important 0.48/0.67 0.57/0.84 —/<0.01 —/0.68 contributors to marine primary productivity16,17. Albeit the contribution of wildfire elements to the ocean biogeochemical cycle remains unclear, limited observations have indicated their significance. For example, atmospheric nitrogen deposition dur- Core 200: BSi/BC ing the 2006 Indonesian wildfires was three to eight times higher ETS (n = 54) <0.01/<0.01 <0.01/<0.01 than during non-fire periods and supported the observation of 0.50/0.80 0.56/0.88 continuously increasing MPP in Sumatra34. More recently, —/<0.01 —/0.61 anomalously widespread phytoplankton blooms were observed from December 2019 to March 2020 in the Southern Ocean downwind of Australia, attributed to aerosol transport of 2019–2020 Australia wildfires, and high iron contents were observed in the aerosol samples35. Iron is an important trace nIOD (n = 31) >0.05/>0.05 >0.05/>0.05 −0.04/0.35 −0.30/0.28 element stimulating phytoplankton growth in oligotrophic —/>0.05 oceans36. Wildfire is also a major source of soluble iron into the —/0.19 ocean via atmospheric aerosols and river input37. In this study, iron and potassium (a typical element in ash) contents in core 200 were measured and they showed distinct increases during pIOD phase (Supplementary Fig. 5). Although the bioavailability of iron pIOD (n = 25) and potassium can reduce their significance correlated with BC in −0.08/−0.24 −0.32/−0.63 >0.05/>0.05 >0.05/<0.01 the ocean, their synchronous increase in the pIOD phase and —/<0.05 positive correlation (Iron: r = 0.55, p < 0.01; potassium: r = 0.72, —/0.41 p < 0.01) with BC (Supplementary Fig. 5) indicates a connection with wildfire occurrences. The results highlight the need, when strong ENSO conditions Core 185: BSi/BC coincide with pIOD phase, to consider the contribution of wild- fire to the functioning of oligotrophic tropical oceans, not just the −0.48/−0.39 −0.42/−0.20 ETS (n = 56) <0.01/>0.05 <0.01/<0.01 role of physical mixing mechanisms (e.g. upwelling, tropical —/<0.01 —/0.41 cyclones). It should be noted that the drying effect of the easterly shift in equatorial trade winds induced by pIOD can promote fire conditions during EDS (January–June) and have an important impact on the ocean. The estimation of atmospheric iron deposition shows that the contribution of iron from wildfire is p p p much higher than that of dust iron, particularly in the equatorial r r r variables). Pacific, which has a significant impact on MPP37. Thus, it is Variables necessary to strengthen observations of atmospheric deposition TOC and riverine input for estimating the flux of wildfire elements BSi TN transported to the ocean. Such understanding will provide further NATURE COMMUNICATIONS | (2022)13:1348 | https://doi.org/10.1038/s41467-022-29013-0 | www.nature.com/naturecommunications 5 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-29013-0 Fig. 3 Burnt areas of early dry season (EDS: January–June) and late dry season (LDS: July–December) in the catchment of King Edward, Drysdale, and King George Rivers from 1988 to 2017 and seasonal patterns of wind and rainfall in the catchment. a EDS and LDS fires in the catchment of King Edward River. b EDS and LDS fires in the catchment of Drysdale River. c EDS and LDS fires in the catchment of King George River. d Climatological monthly wind vector and rainfall during EDS and LDS periods. knowledge on the asymmetry of phytoplankton community excess 210Pb used to obtain the age models was determined as the difference responses to climate change in the northern and southern between the total 210Pb and 226Ra (210Pbsupported). The concentrations of 226Ra were determined for selected samples along the cores through gamma spectro- hemispheres. metry: calibrated geometries in HPGe detectors (CANBERRA, Mod. SAGe Well) were used to measure the decay product of 226Ra, 214Pb, at 295 and 352 keV. These Methods concentrations agreed with the total 210Pb concentrations at depths below the Core collection and chronological analyses. The three sediment cores (cores 185, excess 210Pb horizons (Supplementary Fig. 4). For Core 185 and 200, the Constant 200, and KGR) used in this study were collected from the northern Kimberley coast Flux: Constant Sedimentation (CF:CS) model39 was applied below the mixed layer (14°02′ S, 126°35′ E; 13°53′ S, 126°45′ E; 13°55′ S, 127°19′ E) (Fig. 1b), using a for the 14 to 28 cm sediment core sections, obtaining an average sedimentation rate vibrating head corer (Specialty Devices, Texas, US). Core 185 (=1.31 m long) was of 0.247 ± 0.018 g cm−2 yr−1 for core 185, and 0.338 ± 0.017 g cm−2 yr−1 for core extracted near to the King Edward River mouth in 20 m water depth, core 200 200. For Core KGR, the Constant Rate of Supply (CRS) model40 was applied, (=1.35 m long) was sampled near the Drysdale River mouth in 14 m water depth, obtaining sedimentation rates varying from 0.03 to 0.20 g cm−2 yr−1. and core KGR was sampled near the King George River mouth in 11.8 m water depth (Fig. 1b). The cores were sliced every 0.5 cm along the upper 10 cm, and the rest of the cores were sectioned at 1 cm-thick intervals. All samples were stored in a Geochemical analysis. Grain sizes were measured to understand their impact on freezer at −20 °C before dating and chemical analyses. organic matter preservation, using a Malvern Mastersizer 2000F Laser Particle Sizer. Samples were pre-treated using 10% H2O2 and 10% HCl to remove organic matter and carbonate, respectively. Then, samples were dispersed in a 0.05% Chronological analysis. The concentration profiles of 210Pb in the cores were (NaPO3)6 solution to separate particles for measurement. determined by measuring its granddaughter 210Po radioactive equilibrium using The TOC and TN contents in the sediment core were measured using an alpha spectrometry38 at Edith Cowan University, Australia. The concentration of elemental analyser (FlashSmart NC Soil, Thermo Scientific). Freeze-dried sediment 6 NATURE COMMUNICATIONS | (2022)13:1348 | https://doi.org/10.1038/s41467-022-29013-0 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-29013-0 ARTICLE a b Drysdale River Drysdale River National Park National Park Mt Trafalgar Mt Trafalgar Prince Regent Kimberley Prince Regent Kimberley Nature Reserve Nature Reserve Pearson correlation coefficient -0.5 -0.4 -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 0.4 0.5 Fig. 4 Spatial distribution of Pearson correlation coefficients between Chl-a, BC, and BSi during 2003–2017. a Chl-a and BC. b Chl-a and BSi. (two-sided analysis, n = 17; Inverted black triangles on the colour bar showing significance (p < 0.05) when the correlation coefficients were below −0.456 and over 0.456). Table 2 Relative importance in predicting marine phytoplankton production (MPP). Estimate Std. error t value p value Relative importance (%) Cores 185/200 Cores 185/200 Cores 185/200 Cores 185/200 Cores 185/200 ETS (1926–2017) Intercept 3.790/0.908 4.333/1.313 0.875/0.691 0.386/0.493 — Core 185: r2 = 0.217; p = 0.012; BC 0.495/0.164 0.169/0.033 2.924/4.994 0.005/<0.001 90.03/87.72 F4,52 = 3.594 SST −0.105/0.007 0.144/0.047 −0.730/0.140 0.468/0.889 4.76/2.61 Delta_AICc = 14.24 TCF 0.026/−0.024 0.026/0.011 0.990/−2.298 0.327/0.026 4.45/6.07 AICcWt = 0.00 Rainfall −0.0002/0.0001 0.0001/0.00004 −1.640/2.688 0.107/0.010 0.76/3.60 Core 200: r2 = 0.489; p < 0.001; F4,49 = 11.71 Delta_AICc = 0.00 AICcWt = 1.00 pIOD (1991–2017) Intercept 5.264/3.308 6.648/2.734 0.792/1.210 0.438/0.239 — Core 185: r2 = 0.289; p = 0.129; BC 0.280/0.201 0.207/0.048 1.351/4.155 0.192/<0.001 47.62/61.98 F4,20 = 2.031 SST −0.137/−0.082 0.220/0.097 −0.624/−0.841 0.539/0.409 4.90/1.94 Delta_AICc = 0.00 TCF 0.066/0.005 0.038/0.018 1.769/0.288 0.092/0.776 45.48/1.07 AICcWt = 0.98 Rainfall −8.30 × 10−5/8.06 × 10−5 1.22 × 10−4/5.52 × 10−5 −0.682/1.461 0.503/0.158 1.99/35.01 Core 200: r2 = 0.543; p < 0.001; F4,23 = 6.821 Delta_AICc = 41.83 AICcWt = 0.00 nIOD (1926–1990) Intercept −2.142/3.527 5.855/2.061 −0.366/1.711 0.717/0.101 — Core 185: r2 = 0.038; p = 0.890; BC 0.159/−0.016 0.298/0.108 0.532/−0.152 0.599/0.881 36.78/8.24 F4,28 = 0.277 SST 0.153/−0.072 0.205/0.069 0.749/−1.040 0.460/0.310 46.10/7.54 Delta_AICc = 11.02 TCF −0.002/−0.027 0.035/0.017 −0.065/−1.576 0.949/0.129 10.62/61.69 AICcWt = 0.00 Rainfall −3.26 × 10−5/1.51 × 10−4 1.64 × 10−4/7.65 × 10−5 −0.199/1.974 0.844/0.061 6.49/22.54 Core 200: r2 = 0.298; p = 0.087 F4,22 = 2.338 Delta_AICc = 40.67 AICcWt = 0.00 Linear Fixed Effect Models and multiple linear regression, with black carbon (BC), sea surface temperature (SST), tropical cyclones frequency (TCF), and rainfall in explaining the annual variance of biosilicate (BSi) during the entire time series (ETS: 1926–2017), the phases of positive Indian Ocean Dipole (pIOD) dominance (1991–2017) and negative Indian Ocean Dipole (nIOD) dominance (1926–1990), respectively, for Core 185 and Core 200. AIC (Akaike information criterion) model outputs Delta-AICc (difference in AIC score between the best model and the model being compared) and AICcWt (proportion of the total amount of predictive power provided) were used to identify the best model for each core and period (bold numbers indicate the variables contribute greatly to BSi or are significant in Linear Fixed Effect Models). samples were homogenised by grinding and then acidified with 2 M HCl to remove molybdenum blue was measured at 660 nm with pure water as a reference, and inorganic carbonate. The acidified samples were then dried in an oven at 60 °C for linear regression was performed for each sample through five absorbance values to 1 day, and then washed using Milli-Q water before measurement. Aliquots of estimate the absorbance at t = 0 h. The standard curve was drawn according to the approximately 50 mg of the pre-treated samples were used for analysis, and the known silicon concentration in the solution and the corresponding absorbance absolute error for the measurements was <0.3%. value42. The concentration of BSi in the sediment was calculated from the BSi measurements were performed using an INESA-L8 ultraviolet-visible concentration of silicon in the solution and the pre-weighted mass of the sediment spectrophotometer, according to the Silicon-Molybdenum Blue method41,42. sample. Solutions of 10% H2O2 and 10% HCl were added to freeze-dried and well-milled Sedimentary BC was quantified using the wet-chemical pre-treatment sediment samples (0.1 g) to remove any organic matter and carbonate. The washed integrated with the thermal optical reflectance (TOR) method43, and the TOR and dried samples were digested with 40 ml 2 M Na2CO3 at 85 °C for 5 h. At each method has been proven to effectively discriminate between char and soot44,45, the hour, 0.1 ml of the solution was extracted for absorbance measurements. The two subtypes of BC with different formation mechanisms and physicochemical volume of each sample was determined to be 10 ml with Milli-Q water, and then properties. Prior to the wet-chemical pre-treatment, the sediment samples were add 0.2 ml of HCl, 0.4 ml of ammonium molybdate solution (10.0 g/100 ml), 0.4 ml firstly thawed, freeze-dried, and homogenised. Approximately 0.20 g of each of ethanedioic acid solution (7.5 g/100 ml), and 0.4 ml of L-Ascorbic acid solution sample was then digested with HCl/HF to eliminate the inorganic carbonates, (5.0 g/100 ml) to make the reaction of silicon-molybdenum blue. The absorbance of metals and metal oxides, and silicates. The remaining residue was filtered through NATURE COMMUNICATIONS | (2022)13:1348 | https://doi.org/10.1038/s41467-022-29013-0 | www.nature.com/naturecommunications 7 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-29013-0 pre-fired (at 450 °C for 4 h) 47 mm-diameter quartz fibre filters (Whatman) and (0.01° × 0.01°) between 2003 and 2017, corresponding to the annual contents of BC dried in a constant temperature and humidity chamber following standard and BSi in the sediment core. methods45. BC was separated and quantified on a thermal optical carbon analyser manufactured by Desert Research Institute, Chinese Academy of Sciences. During Reporting summary. Further information on research design is available in the Nature the analysis, a 0.5 cm2 circular filter punch was placed in an oven. First, the oven Research Reporting Summary linked to this article. was heated in 100% He environment for pyrolysis of organic carbon (defined as OCPyro) and it can be monitored by the 635 nm diode laser. Then the analytical atmosphere was shifted to a mixture of 2% O2 and 98% He, three BC sub-fractions Data availability (defined as BC1, BC2, and BC3) will be generated respectively in three temperature The authors declare that the data supporting the findings of this study are available stages (580, 740, and 840 °C). All released carbon fractions were oxidised to CO2 within the paper and as a Supplementary Information file and may also be requested with MnO2 as the catalyst and estimated using a non-dispersive infrared detector from D.L. Source data are provided with this paper. (NDIR). BC was calculated as the sum of three BC sub-fractions minus the OCPyro (i.e. BC = BC1 + BC2 + BC3 − OCPyro). For quality assurance and control (QA/ QC), a random selection and analysis of 10% of the total filters showed that the Code availability relative standard deviation (RSD, %) of the measured BC concentration from There are no custom codes or custom mathematical algorithms used in this paper. different positions within a similar filter was less than 10%, demonstrating the even distribution of the residues onto the filters. In addition, blanks and replicate samples were analysed simultaneously at a frequency of one per ten samples. The Received: 15 June 2021; Accepted: 15 February 2022; blank samples yielded non-detectable BC, and the RSD of replicate analysis averaged 5%. The thermograms of a series of char and soot reference materials show that char evolves almost exclusively in the element carbon one (EC1) stage and soot in the EC2 and EC3 stages46,47. Char content is thus calculated as EC1 minus OCPyro and soot content as EC2 plus EC3. The high-temperature fossil fuel combustion, such References as motor vehicle emissions and industrial coal combustion, show a char/soot ratio 1. Field, C. B., Behrenfeld, M. J., Randerson, J. T. & Falkowski, P. Primary of typically less than 1.0, while the relatively low-temperature biomass burning production of the biosphere: integrating terrestrial and oceanic components. yielded a char/soot ratio significantly higher than 1.0, ranging from 1.2 to ~68 Science 281, 237–240 (1998). (depending on the fuel type, combustion temperature, air/fuel ratio, and so forth) 2. Doney, S. C. Plankton in a warmer world. Nature 444, 695–696 (2006). (Han et al., 2010). The ratios char/soot in the three cores ranges from 6.9 to 22.3 for 3. Moore, J. K. et al. Sustained climate warming drives declining marine core 185, from 9.8 to 26.2 for core 200, and from 12.3 to 31.7 for core KGR, with biological productivity. Science 359, 1139–1143 (2018). averages of 13.4 ± 3.3, 15.5 ± 3.5, and 21.0 ± 4.4 respectively (Supplementary Fig. 1). 4. Gregg, W. W. & Rousseaux, C. S. Global Ocean primary production trends in Given the Kimberley characteristics, where there is sparse population, low human the modern ocean color satellite record (1998-2015). Environ. Res. Lett. 14, activity, and massive distance to major cities, together with the relatively higher 124011 (2019). char/soot ratios, BC in the three cores is dominantly from biomass burning. 5. Sharma, P., Marinov, I., Cabre, A., Kostadinov, T. & Singh, A. Increasing Potassium and iron were measured following EPA method 305248. For each biomass in the warm oceans: unexpected new insights from SeaWiFS. sample, 0.25 g dried and homogenised bulk soil was digested at 240 °C using a Geophys. Res. Lett. 46, 3900–3910 (2019). mixture of concentrated HNO3, HF, and HCl. After total digestion, the samples 6. Behrenfeld, M. J. et al. Climate‐driven trends in contemporary ocean were diluted with 1% HCl prior to injection into the ICP-MS. Reagent blanks and productivity. Nature 444, 752–755 (2006). standard reference materials NIST 2702 (Inorganics in Marine Sediment) and 7. Laufkötter, C. et al. Drivers and uncertainties of future global marine primary MESS-3 (National Research Council of Canada) were run in parallel to the production in marine ecosystem models. Biogeosciences 12, 6955–6984 (2015). samples. 8. Tréguer, P. et al. Influence of diatom diversity on the ocean biological carbon pump. Nat. Geosci. 11, 27–37 (2018). 9. Marinov, I. et al. North-South asymmetry in the modeled phytoplankton Climate data and satellite Chl-a. The dipole mode index is used to indicate the community response to climate change over the 21st century. Glob. IOD phase, which was calculated by the anomalous SST gradient between the Biogeochem. Cycles 27, 1274–1290 (2013). western equatorial Indian Ocean (50°E–70°E and 10°S–10°N) and the south- eastern equatorial Indian Ocean (90°E–110°E and 10°S–0°N), based on the ERSST 10. Lewis, K. M., van Dijken, G. L. & Arrigo, K. R. Changes in phytoplankton V5 dataset (https://psl.noaa.gov/data/gridded/data.noaa.ersst.v5.html). The Niño concentration now drive increased Arctic Ocean primary production. Science 3.4 SST index indicates ENSO, which was calculated by the average SST over the 369, 192–202 (2020). area 5°S–5°N and 170°E–120°W. SST and tropical cyclones from 1920 to 2017 on 11. Menkes, C. E. et al. Global impact of tropical cyclones on primary production. the Kimberley coast were obtained from the HadISST1 dataset (https:// Glob. Biogeochem. Cycles 30, 767–786 (2016). www.metoffice.gov.uk/hadobs/hadisst/) and Australian Tropical Cyclone Database 12. Righetti, D., Vogt, M., Gruber, N., Psomas, A. & Zimmermann, N. E. Global (http://www.bom.gov.au/), respectively. Rainfall data were obtained from two pattern of phytoplankton diversity driven by temperature and environmental weather stations in Kalumburu (14.30°S, 126.65°E; http://www.bom.gov.au/). The variability. Sci. Adv. 5, eaau6253 (2019). monthly mean wind vector was obtained from NCEP Reanalysis data (https:// 13. Jolly, W. M. et al. Climate-induced variations in global wildfire danger from psl.noaa.gov/data/gridded/data.ncep.reanalysis.html). Fire records with a resolution 1979 to 2013. Nat. Commun. 6, 7537–7537 (2015). of 1.1 km pixels were obtained from NOAA Monthly Fire Burnt Areas 1988–2018. 14. Vitolo, C. et al. ERA5-based global meteorological wildfire danger maps. Sci. Satellite Chl-a data from 2003 to 2017 on the Kimberley coast were MODIS Aqua Data 7, 216 (2020). 4 km Chlorophyll-a level3 products derived from the Asia-Pacific Data Research 15. Etiegni, L. & Campbell, A. G. Physical and chemical characteristics of wood Center public dataset (http://apdrc.soest.hawaii.edu/data/data.php). ash. Bioresour. Technol. 37, 173–178 (1991). 16. Wiedinmyer, C. et al. The Fire Inventory from NCAR (FINN): a high resolution global model to estimate the emissions from open burning. Geosci. Statistical analysis. The correlations between each pair of variables in Table 1 Model Dev. 4, 625–641 (2011). were detected by Pearson correlation analysis, with r representing the correlation 17. Jaffé, R. et al. Global charcoal mobilization from soil via dissolution and coefficient and p representing the level of significance (significant: p < 0.05; highly riverine transport to the oceans. Science 340, 345–347 (2013). significant: p < 0.01; not significant: p > 0.05). Linear Fixed Effect Models, using 18. Harris, S. & Lucas, C. Understanding the variability of Australian fire weather Ime4 script in R49 and multiple regression analysis, were used to quantify the between 1973 and 2017. PLoS ONE 14, e0222328 (2019). relative importance of each variable in predicting BSi in Table 2. AIC (Akaike 19. Cai, W. et al. Opposite response of strong and moderate positive Indian Ocean information criterion) model selection was used to distinguish among a set of Dipole to global warming. Nat. Clim. Change 11, 27–32 (2021). possible models describing the relationships among the variables studied using 20. Khan, M. Z., Feng, Q., Zhang, K. & Guo, W. Biogenic silica and organic AICcmodavg script in R50. Normality and homoscedasticity of model residuals carbon fluxes provide evidence of enhanced marine productivity in the Upper were estimated visually. Ordovician-Lower Silurian of South China. Palaeogeogr. Palaeoclimatol. The sequential t test analysis of regime shifts (STARS)51 was used to examine Palaeoecol. 534, 109278 (2019). the shifting time and magnitude of geochemical parameters in Fig. 2. The STARS 21. Nelson, D. M., Treguer, P., Brzezinski, M. A., Leynaert, A. & Queguiner, B. algorithm converted to VBA for Excel is available from Production and dissolution of biogenic silica in the ocean: revised global www.BeringClimate.noaa.gov. Briefly, the method uses a t test to determine if estimates, comparison with regional data and relationship to biogenic sequential records in a time series depart significantly (p < 0.05) from mean values sedimentation. Glob. Biogeochem. Cycles 9, 359–372 (1995). in the preceding period (cut-off length set to 10 years). The regime shift index (RSI) 22. Yuan, Z. et al. Phytoplankton responses to climate‐induced warming and represents a cumulative sum of the normalised anomalies, indicating the interdecadal oscillation in North‐Western Australia. Paleoceanogr. magnitude of the shift. The spatial distribution of the Pearson correlation Paleoclimatol. 35, e2019PA003712 (2020). coefficients in Fig. 3 was determined based on the annual Chl-a in each grid 8 NATURE COMMUNICATIONS | (2022)13:1348 | https://doi.org/10.1038/s41467-022-29013-0 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-29013-0 ARTICLE 23. Goldberg, E. D. Black Carbon in the Environment: Properties and Distribution 50. Mazerolle, M. J. AICcmodavg: model selection and multimodel inference (Wiley, 1985). based on (Q)AIC(c). R package version 2.3-1 https://cran.r-project.org/web/ 24. Han, Y. M., Cao, J. J., Lee, S. C., Ho, K. F. & An, Z. S. Different characteristics packages/AICcmodavg/index.html (2020). of char and soot in the atmosphere and their ratio as an indicator for source 51. Rodionov, S. N. Use of prewhitening in climate regime shift detection. identification in Xi’an, China. Atmos. Chem. Phys. 10, 595–607 (2010). Geophys. Res. Lett. 33, L12707 (2006). 25. Halpern, B. S. et al. A global map of human impact on marine ecosystems. Science 319, 948–952 (2008). 26. Landgate. NOAA Monthly Fire Burnt Areas 1988–2018. Data Extracted Acknowledgements November 2019. (Landgate, 2019). This work was supported by the National Natural Science Foundation of China 27. Perry, J. J. et al. Regional seasonality of fire size and fire weather conditions (42030402, 41876127) and the Western Australian Marine Sciences Institution. O.S. was across Australia’s northern Savanna. Int. J. Wildl. Fire 29, 1–10 (2020). supported by I + D + i projects RYC2019-027073-I and PIE HOLOCENO 20213AT014 28. Corey, B. & Radford, I. J. North-Kimberley Landscape Conservation Initiative: funded by MCIN/AEI/10.13039/501100011033 and FEDER. Funding was provided to 2015–16 Monitoring, Evaluation, Research & Improvement Report (WA P.M. through an Australian Research Council LIEF Project (LE170100219). This work is Department of Parks and Wildlife, 2017). contributing to the ICTA “Unit of Excellence” (MinECo, MDM2015-0552). The IAEA is 29. Smith, A. J. P., Jones, M. W., Abatzoglou, J. T., Canadell, J. G. & Betts, R. A. grateful for the support provided to its Environment Laboratories by the Government of Climate Change Increases the Risk of Wildfires. ScienceBrief Review. https:// the Principality of Monaco. Fieldwork was carried out from the RV Solander (KGR) and sciencebrief.org/topics/climate-change-science/wildfires (2020). the PV Worndoom (Napier Broome Bay) and we thank Daryl Moncrieff, Jennifer 30. Lavender, S. L. & Walsh, K. J. E. Dynamically downscaled simulations of Munro, Matt de Candia, Michael Hourn, and Ryan Crossing for facilitating and/or Australian region tropical cyclones in current and future climates. Geophys. assisting with fieldwork. Georgina Pitt and Yueqi Wang assisted with data analysis. Res. Lett. 38, L10705 (2011). Zineng Yuan assisted with core sectioning; Peter Scott analysed the metals’ content in the 31. Thompson, P. A. et al. Climate variability drives plankton community sediment samples. Peter Thompson provided scientific comments for manuscript. All composition changes: The 2010–2011 El Niño to La Niña transition around necessary site access and collecting approvals were obtained from the Western Australian Australia. J. Plankton Res. 37, 1–19 (2015). Departments of Biodiversity Conservation and Attractions and Primary Industries and 32. Furnas, M. J. & Carpenter, E. J. Primary production in the tropical Regional Development and Balanggarra traditional owner representatives William continental shelf seas bordering northern Australia. Cont. Shelf Res. 129, Maraltadj Jnr and Bradley Carlton accompanied and assisted in the fieldwork conducted 33–48 (2016). in Napier Broome Bay. 33. Siswanto, E. Atmospheric deposition—another source of nutrients enhancing primary productivity in the eastern tropical Indian Ocean during positive Author contributions Indian Ocean Dipole phases. Geophys. Res. Lett. 42, 5378–5386 (2015). D.L. formulated concept and undertook the data analysis and writing; C.Z. measured 34. Sundarambal, P., Balasubramanian, R., Tkalich, P. & He, J. Impact of biomass geochemical parameters and carried out statistical analysis; J.K.K. and D.L. carried out burning on ocean water quality in Southeast Asia through atmospheric field programme and core collection; Y.F. and Y.C. undertook the black carbon analysis; deposition: field observations. Atmos. Chem. Phys. 10, 11323–11336 (2010). O.S., A.W. and A.S. undertook statistical analysis and interpretation of potassium and 35. Tang, W. et al. Widespread phytoplankton blooms triggered by 2019-2020 iron data; G.P. and J.K. sourced and analysed the burnt area data; P.M. carried out the Australian wildfires. Nature 370, 370–375 (2021). chronological analyses and Y.D. provided IOD and ENSO data. 36. Tagliabue, A. et al. The integral role of iron in ocean biogeochemistry. Nature 543, 51–59 (2017). 37. Hamilton, D. S. et al. Impact of changes to the atmospheric soluble iron Competing interests deposition flux on ocean biogeochemical cycles in the Anthropocene. Glob. The authors declare no competing interests. Biogeochem. Cycles 34, e2019GB006448 (2020). 38. Sánchez‐Cabeza, J., Masqué, P. & Ani‐Ragolta, I. 210Pb and 210Po analysis in sediments and soils by microwave acid digestion. J. Radioanal. Nucl. Chem. Additional information Supplementary information The online version contains supplementary material 227, 19–22 (1998). available at https://doi.org/10.1038/s41467-022-29013-0. 39. Krishnaswamy, S., Lal, D., Martin, J. M. & Meybeck, M. Geochronology of lake sediments. Earth Planet. Sci. Lett. 11, 407–414 (1971). Correspondence and requests for materials should be addressed to Dongyan Liu, John K. 40. Appleby, P. G. & Oldfield, F. The calculation of lead-210 dates assuming a Keesing or Yingjun Chen. constant rate of supply of unsupported 210Pb to the sediment. Catena 5, 1–8 (1978). Peer review information Nature Communications thanks Beizhan Yan, Maxine Mowe, 41. Mortlock, R. A. & Froelich, P. N. A simple method for the rapid determination Dan Hammarlund, and the other anonymous reviewer(s) for their contribution to the of biogenic opal in pelagic marine sediments. Deep Sea Res 36, 1415–1426 peer review of this work. Peer reviewer reports are available. (1989). 42. DeMaster, D. J. The supply and accumulation of silica in the marine Reprints and permission information is available at http://www.nature.com/reprints environment. Geochim. Cosmochim. Acta 45, 1715–1732 (1981). 43. Han, Y. et al. Evaluation of the thermal/optical reflectance method for Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in discrimination between char- and soot-EC. Chemosphere 69, 526–533 (2007). published maps and institutional affiliations. 44. Fang, Y. et al. Flux and budget of BC in the continental shelf seas adjacent to Chinese high BC emission source regions. Glob. Biogeochem. Cycles 29, 957–972 (2015). 45. Han, Y. M. et al. Evaluation of the thermal/optical reflectance method for Open Access This article is licensed under a Creative Commons discrimination between char- and soot-EC. Chemosphere 69, 569–574 (2007). Attribution 4.0 International License, which permits use, sharing, 46. Fang, Y. et al. Spatiotemporal trends of elemental carbon and char/soot ratios adaptation, distribution and reproduction in any medium or format, as long as you give in five sediment cores from Eastern China Marginal Seas: indicators of appropriate credit to the original author(s) and the source, provide a link to the Creative anthropogenic activities and transport patterns. Environ. Sci. Technol. 52, Commons license, and indicate if changes were made. The images or other third party 9704–9712 (2018). material in this article are included in the article’s Creative Commons license, unless 47. Han, Y. M., Cao, J. J., Lee, S. C., Ho, K. F. & An, Z. S. Different characteristics indicated otherwise in a credit line to the material. If material is not included in the of char and soot in the atmosphere and their ratio as an indicator for source article’s Creative Commons license and your intended use is not permitted by statutory identification in Xi’an, China. Atmos. Chem. Phys. 10, 595–607 (2010). regulation or exceeds the permitted use, you will need to obtain permission directly from 48. US EPA. Method 3052: microwave assisted acid digestion of siliceous and the copyright holder. To view a copy of this license, visit http://creativecommons.org/ organically based matrices. Revision 0 https://www.epa.gov/sites/default/files/ licenses/by/4.0/. 2015-12/documents/3052.pdf (1996). 49. Bates, D., Maechler, M., Bolker, B. & Walker, S. Fitting Linear Mixed-Effects Models using lme4. J. Stat. Softw. 67, 1–48 (2015). © The Author(s) 2022 NATURE COMMUNICATIONS | (2022)13:1348 | https://doi.org/10.1038/s41467-022-29013-0 | www.nature.com/naturecommunications 9