5 kyr of fire history in the High North Atlantic Region: natural variability and ancient human forcing

Biomass burning influences global atmospheric chemistry by releasing greenhouse gases and climate-forcing aerosols. There is controversy about the magnitude and timing of Holocene changes in biomass burning emissions from millennial to centennial time scales and, in particular, on the possible impact of ancient civilizations. Here we present a 5 kyr record of fire activity proxies levoglucosan, black carbon and ammonium measured in the RECAP ice core, drilled in the coastal East Greenland and therefore affected by processes occurring in the High North Atlantic Region. Levoglucosan and ammonium fluxes 5 show high levels from 5 to 4.5 kyr followed by an abrupt decline, possibly due to monotonic decline in Northern Hemisphere summer insolation. Levoglucosan and black carbon show an abrupt decline at 1.1 kyr BP (before 2000 AD), suggesting a decline in wildfire regime in the Icelandic territory due to the extensive land clearing caused by Viking colonizers. A minimum is reached at 0.5 kyr BP for all fire proxies, after which levoglucosan and ammonium fluxes increase again, in particular over the last 200 years. We find that the fire regime reconstructed from RECAP fluxes seems mainly related to climatic changes, 10 however over the last millennium human activities might have had a substantial influence controlling the occurrence of fire.

continuous biogenic emissions from vegetation and soils, while only peak values can be associated with biomass burning events (Fischer et al., 2015). K + , additionally, has been found to be highly sensitive to contamination, making it difficult to measure its species in ice (Legrand et al., 2016). Recently, however, most of the attention has been given to two specific fire proxies levoglucosan and black carbon (BC), the latter being a specific tracer of biomass burning in the pre-industrial times (Osmont et al., 2019). Despite that NH 4 + is influenced by biogenic emissions from Greenland ice-free areas, it has been included in this (MilliQ) and a buffer made from 35.8 g Na 2 HPO 4 ·12H 2 O, 1000 mL MilliQ, 600 µL NaOH (>32%), 100 µL HCHO (>37%) and 0.8 g Na 2 SO 3 , following 1 m of mixing coil at 80°C and 0.2 m mixing coil at room temperature the mixture was excited at 365 nm and detected at 400 nm by means of a photomultiplier based detector (PMT-FL, FIAlab instruments). Three standards 120 were used for calibration based on an 100 mg L -1 IC multielement standard (VII, Certipur, Merck) and diluted to 24, 50 and 200 ppb NH 4 + . The melt rate of the CFA was kept between 4 and 5 cm min -1 and with a response time of 12 seconds; the equivalent depth resolution of the NH 4 + dataset is less than 1 cm.
Based on the age scale released in Simonsen et al. (2019), the levoglucosan record covers the period from 0.089 to 5 kyr BP and a depth of 482 m, the BC record covers the period from 0.4 to 5 kyr BP and a depth of 482 m and the NH 4 + covers the 125 entire period back to 5 kyr.

Potential source regions of RECAP fluxes
To identify the potential forest fires regions able to influence the RECAP site, backward trajectories have been computed using the Hybrid Single-Particle Lagrangian Integrated Trajectory (HYSPLIT) model (Stein et al., 2015) with GDAS1 meteorological data. GDAS1 dataset is available at 1 x 1°resolution. The HYSPLIT model was run using the PySPLIT Python package 130 (Warner, 2018). Trajectories were calculated every 6 h in backward mode at 3 different altitudes (500, 1000 and 2000 m above Renland elevation, 2315 m a.s.l.) for the entire dataset available, that is the period 2006-2019. A 7-days run time was chosen based on estimations of BC maximum residence time (Ramanathan and Carmichael, 2008;Cape et al., 2012).
HYSPLIT modeling tool is useful to reconstruct the most likely transport areas over recent period and the information derived can be extended, over a longer timescales (Pre-industrial or Holocene), assuming that the main atmospheric circulation 135 mode has not changed significantly during those periods (Rubino et al., 2015). Stable isotope data (δ 18 O) suggest that the circulation pattern of air masses reaching the Greenland Ice Sheet did not significantly change over the last 10 000 years (Vinther et al., 2006) which supports extending to the past back trajectory analysis based on modern conditions.

140
Seven day long back trajectories calculated with HYSPLIT model over the period 2006-2019 suggest that the higher frequency of trajectories come from the nearby areas with respect to the Renland site (Figure 1), as also evidenced by Simonsen et al. (2019) and Maffezzoli et al. (2019) studying insoluble dust and sea ice respectively. In addition to the Renland ice cap surrounding area, the coastal area immediately on the North of the Renland site appears to be the most probable source region.
The South-East and South-West coastal areas of Greenland and Iceland are also ice-free source areas which provide the densest 145 and shortest trajectories considering the density of gridded points and travelling time. We define the region ranging from the longitudes 60°W -0°and the latitudes 90°N -55°N as the region with the most likely source areas of materials arriving at Renland (red box in Figure 1) and we hereafter call it High North Atlantic Region (HNAR).
Other contributors of air masses are North America, Northern Europe and Siberia. Boreal forests of North America and Siberia are also a possible source of impurities, as well as Northern Europe (Schüpbach et al., 2018), however they require a 150 longer travel path for fire emissions to reach the Renland site compared to coastal Greenland and Iceland and thus are expected to carry only a minor contribution.

Levoglucosan, BC and NH 4 + results
RECAP levoglucosan, BC and NH 4 + profiles display distinct variability on centennial timescales and do not exhibit a clear trend over the past 5 kyr ( Figure S1). The levoglucosan concentrations vary from 0.006 to 0.1 ng g -1 with a mean (considering 155 the whole 5 kyr record) of 0.033 ± 0.002 ng g -1 , BC concentration profile varies from 0.4 to 1.7 ng g -1 with a mean of 0.942 ± 0.006 ng g -1 and NH 4 + varies from 0.09 to 17 ng g -1 with a mean of 5.0 ± 0.2 ng g -1 .
Since Levoglucosan, BC and NH 4 + can be deposited in snow and glaciers by both wet and dry deposition (Stohl et al., 2007), fluxes were calculated as the product of concentration with annual accumulation and expressed in µg m -2 y -1 . Renland accumulation (Hughes et al., 2020;Corella et al., 2019;Simonsen et al., 2019) shows a stable profile, indicating that wet Levoglucosan flux shows high levels from 5 to 4.5 kyr BP with 31.5 µg m -2 y -1 on average. A rather stable deposition flux is determined between 4 to 2 kyr BP while higher variability is determined during the period 2 to 1 kyr BP, with an average of 165 17 µg m -2 y -1 . Lower fluxes are found in the period 1 to 0.5 kyr BP with an average of 10 µg m -2 y -1 . BC flux is stable with an average of 464 µg m -2 y -1 over the period 5 to 1 kyr BP ( Figure 2b). During this period some oscillations are detected with a higher depositional flux determined at ca. 3.6, 3, 2.5, 2 kyr BP. A decrease in BC flux is determined, similarly to levoglucosan, after 1.1 kyr BP with a value of 372 µg m -2 y -1 . NH 4 + flux shows higher values from 5 to 4.6 kyr BP with an average of 2889 µg m -2 y -1 and a rather stable profile from 4 to 1.25 kyr BP (2387 µg m -2 y -1 on average) and an abrupt decline at 2.7 kyr BP

Statistical analysis
To investigate possible similarities with available time series, we estimated pairwise Pearson correlations among available records (Table 1). Correlation coefficients were computed after interpolating all series to obtain a 20-years time resolution.

175
20-year bins have been chosen in order to maximize time resolution among all time series, limited by Northern Hemisphere temperature (Marcott et al., 2013). We investigated the correlation between our series with NEEM levoglucosan flux (Zennaro et al., 2014), Northern Hemisphere temperatures as tracer of main climate variability (Marcott et al., 2013), RECAP δ 18 O (Hughes et al., 2020) and sedimentary charcoal influx composites as regional fire reconstruction calculated over North Europe, North America, Siberia and more generally for the entire HLNH (High Latitude Northern Hemisphere, lat > 55°), available fluxes is low but significant (r = 0.21, p < 0.01), as well as between BC and NH 4 + (r = 0.27, p < 0.01). RECAP levoglucosan flux is significantly correlated with Northern Hemisphere reconstructed temperature (r = 0.46, p < 0.001) as well as with RECAP δ 18 O (r = 0.33, p < 0.01); RECAP BC flux is also significantly correlated with NH temperature (r = 0.32, p < 0.001) and with RECAP δ 18 O (r = 0.29, p < 0.01). No statistically significant correlation has been determined between RECAP biomass 185 burning fluxes and the regional sedimentary charcoal influx composites (Table 1).
To understand the mechanisms behind fire regime changes, we examine fire flux step changes and link them with climate and human history. To determine mean changes in time series we conduct the off-line change point analysis using the ruptures Python package (Truong et al., 2019). Three optimal breakpoints were found for levoglucosan flux using the Pelt method; Dynamic programming was then applied setting three breakpoints. The breakpoints of each time series are showed in Figure   190 3 by a change of color. As for levoglucosan flux, two main step changes are found at 4.45 and 1.15 kyr, where the profile significantly decreases, while at 0.2 kyr BP it increases; regarding the BC flux, one strong step decline is found at 1.03 kyr BP, while for the NH 4 + flux the main step changes are found at 4.57, 4.37 and 0.12 kyr BP (Figure 3a,c,e). It was not possible to detect the concomitant shift at 1.15 kyr BP on NH 4 + flux due to absence of data.
As previously discussed, climate and especially summer temperatures have a central role in explaining the trends observed 195 in RECAP fire reconstruction. However, other processes might be of relevance and in order to disentangle the contribution of climate in changing the fire regime, we calculate ratios between levoglucosan and BC fluxes with Northern Hemisphere reconstructed temperature from Marcott et al. (2013): a stable ratio suggests that climate is the main driver, significant variations in the ratio suggest additional processes to be in place. We perform change point analysis also on ratios between normalized levoglucosan and BC fluxes and Northern Hemisphere temperature (Levo/T and BC/T in Figure 3). The RECAP Levo/T ratio  Several studies suggest that climate is the main driver of global biomass burning (Marlon et al., 2008). Elevated summer temperatures and sustained droughts can affect fuel flammability and lead to increased global fire activity over seasonal to centennial timescales (Daniau et al., 2012). Climate conditions also determined forest expansion (and glacier retreat), with positive effects on fuel moisture and a dampening effect on biomass burning (Feurdean et al., 2020). Significant correlations of RECAP role in controlling forest fire activity in the HNAR. Fire reduced throughout the late Holocene in parallel with progressively decreasing summer insolation in the Northern Hemisphere through the late Holocene (Power et al., 2008). Greater-than-present summer insolation resulted in warmer and drier summers in the Northern Hemisphere with a consequent increasing of fire as showed by records from North America and Europe. Renland, being located at a high latitude, receives minimal or no insolation throughout the winter, meaning that summer insolation dominates (Hughes et al., 2020). Thus, summer solar input along 220 with temperatures may have the major influence in controlling fire activity in the HNAR.
High levels of levoglucosan and NH 4 + fluxes from 5 to 4.5 kyr BP could be linked to the NH warm temperatures of the earlyto-mid Holocene. In fact, Marcott et al. (2013) find warmer temperatures from 9.5 to 5.5 kyr B.P, followed by a cooling trend from 5.5 kyr BP onwards. Such cooling trend could explain levoglucosan and NH 4 + trends from 5 to 4.5 kyr BP. An additional source of NH 4 + fluxes, however, could also be soil emission from Greenland ice-free areas, especially during summer when decrease in biomass burning from the late nineteenth century to mid-to-late twentieth century (Marlon et al., 2008), however this is below the resolution of the RECAP fire tracer fluxes.
RECAP levoglucosan and BC fluxes do not correlate with NEEM levoglucosan flux (Table 1)  periods of Holocene, when the ice sheet retreats and the coasts are subject to the smoldering of peats.
Based on the absence of correlation between fire tracer fluxes and any charcoal record from the regions investigated as possible sources, the low correlation with NEEM levoglucosan flux and the results from back-trajectory analysis we suggest that the main source of fire markers that are deposited to Renland might the Icelandic region. Iceland has a high aeolian activity driven by climate, volcanic activity and glacial sediment supply. Atmospheric low-pressure systems are common, sometimes 260 referred as the "Icelandic low", which frequently result in relatively high wind speeds (Arnalds et al., 2016). It is thus likely that RECAP fire proxies respond more to a local influence rather than a hemispheric trend and are affected by fire emissions from Iceland. The fire history of Iceland has never been documented and the RECAP ice core might represent the first record preserving past Icelandic fire changes. The significant link between levoglucosan and BC with temperature reconstruction as well as RECAP δ 18 O suggests that the change of Icelandic fire regime might have been mostly driven by climate fluctuations.

265
The presence of birch forests covering the Icelandic territory is well documented for the last 2.5 kyr (Streeter et al., 2015) and the expansion of vegetation is strongly dependent on climate (Haraldsson and Ólafsdóttir, 2003).
Although the Icelandic territory is characterized by high humidity and frequent precipitation (Ólafsson et al., 2007)

Additional processes explaining fire proxy variability
Global climate alone may not be the only driver of fire activity. Additional processes are suggested to have an impact on the fire regime (Marlon et al., 2008), such as the anthropogenic activities including land clearance, ice sheet advance/retreat and peat fires accentuated by permafrost thawing, altering the temporal and spatial structure of fuel and the frequency of ignitions since the early Holocene (Pfeiffer et al., 2013;Vannière et al., 2016;Andela et al., 2017).

275
From change point analysis conducted on fire fluxes and on ratios with temperature we identify three main changes at approximately 4.5, 1.1 and 0.2 kyr BP, found to be common to more than one time series. In the following sections we formulate hypotheses about possible processes inducing fire regime variability.  population and agriculture in many areas of the Northern Hemisphere. Zennaro et al. (2014) found relatively high levoglucosan values in the NEEM ice core during the last two centuries, although a clear assessment was not possible due to lack of data.
Several charcoal records from North America and Europe evidence an increase since~1800 (as also reported in Figure 3f,g,i) and have been attributed both to climate and to the increasing use of fire for land clearance (Marlon et al., 2008(Marlon et al., , 2016. The abrupt increase of NH 4 + in the last 120 years BP may be attributed to anthropogenic emissions. Wendl et al. (2015) 320 argued that trends in NH 4 + concentrations in the Lomonosovfonna-09 ice core (Svalbard) after 1940 AD indicate a strong anthropogenic influence mainly from NH 3 emissions from agriculture and livestock in Eurasia. No BC data are available for the last 500 years, however, high values of BC flux are expected since several shallow cores from Svalbard and Greenland show broad concentration maxima starting from 1860 AD and are attributed to anthropogenic emissions (Osmont et al., 2018).

325
This paper provides 5 kyr records of fire proxies levoglucosan, BC and NH 4 + in the RECAP ice core in Greenland.
From back trajectory analysis and the comparison with regional fire reconstructions based on charcoal records, we find that the most likely source area of impurities arriving at Renland is the High North Atlantic Region (ranging from the longitudes 60°W -0°and the latitudes 90°N -55°N) and comprehends the Greenland Ice Sheet, the coasts of Greenland and Iceland.
Iceland, in particular, could have a great influence in driving fire regime changes in RECAP fire proxies due to its position with 330 respect to the Renland site and to the presence of large forested areas throughout the Holocene. In the Northern Hemisphere, Iceland is also the only region which is not up until now covered by fire reconstructions and RECAP fire proxies might represent the first reconstruction of past Icelandic fire regime.
We find that climate variability is the main control of changes in the High North Atlantic fire regime, and especially temperature driven by summer insolation. A downturn of levoglucosan and NH 4 + fluxes at 4.5 kyr BP may be associated with the 335 monotonic Holocene decline in Northern Hemisphere summer insolation that resulted in cooler conditions, cryosphere expansion and vegetation reduction in the HNAR. BC shows a different trend in this period, possibly explained by wildfires linked to smoldering of peats along the Greenland coasts due to warmer temperatures and melting of permafrost. While levoglucosan and NH 4 + can record peat fires, BC is not a good proxy because is emitted during flaming conditions when biomass burns.
Further studies are necessary to confirm our hypothesis. During the last millennium the fire regime may have been influenced 340 by the human impact in the Icelandic environment, which was uninhabited before the arrival of the Vikings. The massive land clearing and the active fire suppression probably led to a reduction of wildfires. The human impact on the fire regime probably accentuated over the last centuries due to population expansion and increases in land cover conversion rates for agricultural purposes.
Data availability. The dataset will be published in the Pangaea database.  Table 1. Correlation matrix of 5 kyr records with resolution of 20 years. *** indicates p < 0.01, ** indicates p < 0.05 and * indicates p < 0.1.
The colorbar represents the correlation from -1 to 1.