# Earlier snowmelt may lead to late season declines in plant productivity and carbon sequestration in Arctic tundra ecosystems

Donatella Zona, Peter Lafleur, Koen Hufkens, Barbara Bailey, Beniamino Gioli, George Burba, Jordan Goodrich, Anna Liljedahl, Eugénie Euskirchen, Jennifer Watts, et al.

<figure><img src="image_1.png" /></figure>

## To cite this version:

Donatella Zona, Peter Lafleur, Koen Hufkens, Barbara Bailey, Beniamino Gioli, et al. Earlier snowmelt may lead to late season declines in plant productivity and carbon sequestration in Arctic tundra ecosystems. Scientific Reports, 2022, 12 (1), pp.1-10. 10.1038/s41598-022-07561-1. hal-03642434

HAL Id: hal-03642434

https://hal.inrae.fr/hal-03642434

Submitted on 15 Apr 2022

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>
<figure><img src="image_3.png" /></figure>

Check for updates

OPEN

# Earlier snowmelt may lead to late season declines in plant productivity and carbon sequestration in Arctic tundra ecosystems

Donatella Zona<fnref n="1" />,²✉, Peter M. Lafleur³, Koen Hufkens⁴,⁵, Barbara Bailey<fnref n="1" />, Beniamino Gioli⁶, George Burba⁷,⁸, Jordan P. Goodrich⁹, Anna K. Liljedahl¹⁰,¹¹, Eugénie S. Euskirchen¹², Jennifer D. Watts¹¹,¹³, Mary Farina¹¹, John S. Kimball¹³, Martin Heimann¹⁴,¹⁵, Mathias Göckede¹⁴, Martijn Pallandt¹⁴, Torben R. Christensen¹⁶,¹⁷, Mikhail Mastepanov¹⁶,¹⁷, Efrén López-Blanco¹⁶,¹⁸, Marcin Jackowicz-Korczynski¹⁶,¹⁹, Albertus J. Dolman²⁰, Luca Belelli Marchesini²¹,²², Roisin Commane²³, Steven C. Wofsy²⁴, Charles E. Miller²⁵, David A. Lipson<fnref n="1" />, Josh Hashemi<fnref n="1" />, Kyle A. Arndt²⁶, Lars Kutzbach²⁷, David Holl²⁷, Julia Boike²⁸,²⁹, Christian Wille³⁰, Torsten Sachs³⁰, Aram Kalhori³⁰, Xia Song<fnref n="1" />, Xiaofeng Xu<fnref n="1" />, Elyn R. Humphreys³¹, Charles D. Koven³², Oliver Sonnentag³³, Gesa Meyer³³, Gabriel H. Gosselin³³, Philip Marsh³⁴ & Walter C. Oechel<fnref n="1" />,³⁵

Arctic warming is affecting snow cover and soil hydrology, with consequences for carbon sequestration in tundra ecosystems. The scarcity of observations in the Arctic has limited our understanding of the impact of covarying environmental drivers on the carbon balance of tundra ecosystems. In this study, we address some of these uncertainties through a novel record of 119 site-years of summer data from eddy covariance towers representing dominant tundra vegetation types located on continuous permafrost in the Arctic. Here we found that earlier snowmelt was associated with more tundra net CO₂ sequestration and higher gross primary productivity (GPP) only in June and July, but with lower net carbon sequestration and lower GPP in August. Although higher evapotranspiration (ET) can result in soil drying with the progression of the summer, we did not find significantly lower soil moisture with earlier snowmelt, nor evidence that water stress affected GPP in the late growing season. Our results suggest that the expected increased CO₂ sequestration arising from Arctic warming and the associated increase in growing season length may not materialize if tundra ecosystems are not able to continue sequestering CO₂ later in the season.

<footnote n="1">¹Department Biology, San Diego State University, San Diego, CA 92182, USA. ²School of Biosciences, University of Sheffield, Western Bank, Sheffield S10 2TN, UK. ³School of the Environment, Trent University, Peterborough, ON K9L 0G2, Canada. ⁴UMR 1391 ISPA, INRA, 71 Avenue Edouard Bourlaux, 33140 Villenave d'Ornon, France. ⁵Department of Applied Ecology and Environmental Biology, Ghent University, 653, 9000 Ghent, Belgium. ⁶Institute of BioEconomy, IBE, National Research Council (CNR), Via Giovanni Caproni 8, 50145 Firenze, Italy. ⁷LI-COR Biosciences, 4421 Superior St., Lincoln, NE 68504, USA. ⁸The Robert B. Daugherty Water for Food Global Institute, School of Natural Resources, University of Nebraska, Lincoln, NE 68583, USA. ⁹Department of Earth Sciences, University of Waikato, Hillcrest, Hamilton 3216, New Zealand. ¹⁰Water and Environmental Research Center, University of Alaska Fairbanks, Fairbanks, AK 99775-7340, USA. ¹¹Woodwell Climate Research Center, Falmouth, MA 02540, USA. ¹²Institute of Arctic Biology, University of Alaska Fairbanks, Fairbanks, AK 99775, USA. ¹³W.A. Franke College of Forestry & Conservation, The University of Montana, Missoula, MT 59812, USA. ¹⁴Max Planck Institute for Biogeochemistry, 07745 Jena, Germany. ¹⁵Institute for Atmospheric and Earth System Research (INAR)/Physics, Faculty of Science, University of Helsinki, Gustaf Hällströmin katu 2b, 00560, P.O. Box 64, 00014 Helsinki, Finland. ¹⁶Department of Ecoscience, Arctic Research Centre, Aarhus University, Frederiksborgvej 399, 4000 Roskilde, Denmark. ¹⁷Oulanka Research Station, Oulu University, Kuusamo, Finland. ¹⁸Department of Environment and Minerals, Greenland Institute of Natural Resources,</footnote>
3900 Nuuk, Greenland. <sup>19</sup>Department of Physical Geography and Ecosystem Science, Lund University, 22362 Lund, Sweden. <sup>20</sup>The Netherlands Institute for Sea Research, Den Burg, PO Box 59, 1790 AB Texel, The Netherlands. <sup>21</sup>Department of Sustainable Agro-Ecosystems and Bioresources, Research and Innovation Centre, Fondazione Edmund Mach, Via E. Mach 1, 38010 San Michele all'Adige, TN, Italy. <sup>22</sup>Agrarian-Technological Institute, RUDN University, 117198 Moscow, Russia. <sup>23</sup>Department of Earth and Environmental Sciences, Lamont-Doherty Earth Observatory, Columbia University, Palisades, NY 10964, USA. <sup>24</sup>School of Engineering and Applied Sciences, Harvard University, 20 Oxford St., Cambridge, MA 02138, USA. <sup>25</sup>Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena, CA 91109-8099, USA. <sup>26</sup>Earth Systems Research Center, Institute for the Study of Earth, Oceans, and Space, University of New Hampshire, 8 College Road, Durham, NH 03824, USA. <sup>27</sup>Institute of Soil Science, Center for Earth System Research and Sustainability (CEN), Universität Hamburg, 20146 Hamburg, Germany. <sup>28</sup>Permafrost Research, Alfred Wegener Institute, Helmholtz Centre for Polar and Marine Research, Potsdam, Germany. <sup>29</sup>Geography Department, Humboldt-Universität zu Berlin, 10099 Berlin, Germany. <sup>30</sup>GFZ German Research Centre for Geosciences, 14473 Potsdam, Germany. <sup>31</sup>Department of Geography and Environmental Studies, Carleton University, Ottawa, ON K1S 5B6, Canada. <sup>32</sup>Climate and Ecosystem Sciences Division, Lawrence Berkeley National Laboratory (LBNL), Berkeley, CA 94720, USA. <sup>33</sup>Département de Géographie, Université de Montréal, 1375 Avenue Thérèse-Lavoie-Roux, Montréal, QC H2V 0B3, Canada. <sup>34</sup>Department of Geography and Environmental Studies, Wilfrid Laurier University, 75 University Ave W., Waterloo, ON N2S 3C5, Canada. <sup>35</sup>Department of Geography, College of Life and Environmental Sciences, University of Exeter, Exeter EX4 4RJ, UK.<sup>36</sup>email: dzona@sdsu.edu

Climate change is affecting arctic ecosystems through temperature increase¹, hydrological changes², earlier snowmelt³,⁴, and the associated increase in growing season length⁵. Annual arctic air temperature has been increasing at more than double the magnitude of the global mean air temperature increase¹, and terrestrial snow cover in June has decreased by 15.2% per decade from 1981 to 2019⁴. Warming is the main driver of the earlier start of the growing season, and the greening of the Arctic⁶⁻⁸. Arctic greening is associated with enhanced vegetation height, biomass, cover, and abundance⁹. However, the complexity of arctic systems reveals an intricate patchwork of landscape greening and browning⁸,¹⁰,¹¹, with browning linked to a variety of stresses to vegetation⁸ including water stress¹²,¹³. The interconnected changes in temperature, soil moisture, snowmelt timing, etc. can have important effects on the carbon sequestered by arctic ecosystems¹⁴. The reservoir of carbon in arctic soil and vegetation depends on the interaction of two main processes: (1) changes in net CO₂ uptake by vegetation; and (2) increased net loss of CO₂ (from vegetation and soil) to the atmosphere via respiration. Therefore, defining the response of both plant productivity and ecosystem respiration to environmental changes is needed to predict the sensitivity of the net CO₂ fluxes of arctic systems to climate change.

An earlier snowmelt, and a longer growing season, do not necessarily translate into more carbon sequestered by high latitude ecosystems⁵. There is a large disagreement on the response of plant productivity and the net CO₂ uptake to early snowmelt in tundra ecosystems¹⁵⁻¹⁹. A warmer and longer growing season might not result in more net CO₂ uptake if CO₂ loss from respiration increases¹⁶, particularly later in the season, and surpasses the CO₂ sequestered by enhanced plant productivity in northern ecosystems¹⁶,²⁰. Moreover, snowmelt timing and the growing season length greatly affect hydrologic conditions of arctic soils²¹, as well as plant productivity²². Longer non-frozen periods earlier in the year²³ and earlier vegetation greening can increase evapotranspiration (ET), resulting in lower summer soil moisture²⁴⁻²⁶. The complexity in the hydrology of tundra systems arises from the tight link between the water drainage and the presence and depth of permafrost. The presence of permafrost reduces vertical water losses, preventing soil drainage in northern wetlands during most of the summer despite low precipitation input²⁷. Increasing rainfall²⁸ and increased permafrost degradation can increase soil wetness in continuous permafrost regions². Further permafrost degradation (e.g. ice-wedge melting) can increase hydrologic connectivity leading to increased lateral drainage of the landscape and subsequent soil drying²,²⁹.

Given the importance of soil moisture in affecting the carbon balance of arctic ecosystems, and its links with snowmelt timing, in this study, we investigated the correlation between summer fluxes of CO₂ (i.e., net ecosystem exchange (NEE), gross primary productivity (GPP) ecosystem respiration (ER)), ET, and environmental drivers such as soil moisture and snowmelt timing, while controlling for the other most important drivers of photosynthesis and respiration (i.e. solar radiation and air temperature). We expected earlier snowmelt to be correlated with larger ET and lower soil moisture, particularly during peak and late season, consistent with drying associated with a longer growing season. The lower soil moisture with earlier snowmelt should result in a negative correlation between snowmelt timing and GPP, particularly during the peak and late season (when we expect the most water stress), and in a positive correlation between snowmelt timing and ET during the entire growing season. This soil moisture limitation to plant productivity should result in lower net cumulative CO₂ sequestration during the entire summer, because of lower plant productivity if these ecosystems are water-limited due to lower soil moisture with earlier snowmelt.

**Testing the impact of snowmelt timing on the carbon dynamics and hydrology of tundra ecosystems.** The 11 sites were selected as among the longest-running tower sites in the circumpolar Arctic (including 6 to 19 years of fluxes per site and a total of 119 site-years of summer (June to August) eddy covariance CO₂ flux data, Table S1). All sites lie in the zone of continuous permafrost. The sites are representative of dominant tundra vegetation classes (wetland, graminoids, and shrub tundra), together accounting for 31% of all tundra vegetation types (Fig. 1³⁰ and Supplementary Information). Given the complex interactions among different variables (many covarying together), we used a variety of statistical analyses to identify the association between standardized anomalies of NEE, GPP, ER, and ET, and standardized anomalies of the main environ-
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Study sites. Locations of the 11 eddy covariance flux tower sites used in this study. Light blue regions delineate the total Circumpolar Arctic Vegetation Map (CAVM), green regions delineate the subset of CAVM vegetation types represented in this study (including all the vegetation types listed in Table S1). This map was created using QGIS.org, 2020, QGIS 3.10. Geographic Information System User Guide. QGIS Association: https://docs.qgis.org/3.10/en/docs/user_manual/index.html. The dataset used in the map was the CAVM map: CAVM Team. 2003. Circumpolar Arctic Vegetation Map. (1:7,500,000 scale), Conservation of Arctic Flora and Fauna (CAFF) Map No. 1. U.S. Fish and Wildlife Service, Anchorage, Alaska. ISBN: 0-9767525-0-6, ISBN-13: 978-0-9767525-0-9.</figcaption></figure>

mental controls during different periods of the summer corresponding to various stages in seasonal phenology (early season: June, peak season: July, and late season: August). We used a partial correlation analysis to identify if the timing of the snowmelt associates with anomalies of ET, soil moisture, NEE, GPP, ER, atmospheric vapor pressure deficit (VPD), or the Bowen ratio (the ratio between Sensible Heat (H) and Latent Heat (LE)) while statistically controlling for the main meteorological forcing such as air temperature and solar radiation (Methods). Identifying the correlation between ET (and the Bowen ratio) and snowmelt timing is a way to assess water limitation to ecosystems (in addition to testing their response to soil moisture changes), as H, and therefore, the Bowen ratio, are expected to increase with surface drying<sup>31,32</sup>. To identify the association between snowmelt timing, the main environmental variables (i.e., air temperature and solar radiation), and NEE, GPP, ER, and ET over time, we performed a maximum covariance analysis (MCA) on the monthly median standardized anomalies from 2004 to 2019 (a time period when data for most of the sites were available). MCA allowed us to find patterns in two space-time datasets that are highly correlated using a cross-covariance matrix<sup>26</sup>. We retained sites as the unit of variation (i.e., by estimating the standardized anomalies by site for each of the indicated variables, see "Methods"), hence the results of the MCA integrated the site level relationships between each of the variables over time). The goal of this analysis was to identify the most important environmental drivers associated with NEE, GPP, and ER across all the sites over time. MCA is particularly appropriate for this study as it can handle data with gaps and unequal lengths in the datasets. We also tested the relative importance of the abovementioned environmental drivers on the monthly median GPP, ER, and NEE using a linear mixed effect model, including site as a random effect to account for the site-to-site variability. The MCA and the mixed model analyses were conducted to test the relative importance of snowmelt and other variables at different times of the season. Finally, to evaluate the water balance at different times of the season, we estimated the difference between Potential Evapotranspiration (PET) and the actual ET, and the difference between precipitation (PPT) and ET for each of the sites, years, and months (e.g. June, July, and August). This study did not attempt to describe the long-term temporal changes in the anomalies of snowmelt and carbon fluxes, given the short data record available for some of the sites (i.e. less than 10 years, Table S1), but instead focused on understanding the association between environmental variables and the carbon balance at different times of the season. More details of these analyses are included in the Methods.

**Influence of snowmelt timing on NEE, GPP, ER, and hydrological status of tundra ecosystems.** Once statistically controlling for solar radiation and air temperature (in the partial correlation analysis, see "Methods"), we observed a significant positive relationship between the snowmelt timing anomalies and NEE anomalies (i.e. earlier snowmelt was associated with a higher net CO₂ sequestration) in June and July, but a negative correlation in August (Fig. 2a, Table 1). A significant relationship was also found between snowmelt date anomalies and GPP anomalies, with more positive GPP anomalies (i.e. higher plant productivity) with earlier snowmelt in June and July, and more negative GPP anomalies with earlier snowmelt in August (Fig. 2b, Table 1). Earlier snowmelt was associated with significantly higher ER in both June and July, but there was no significant relationship in August (Fig. 2c, Table 1), suggesting that the late-season correlation between NEE and snowmelt timing was mostly driven by the lower GPP and with earlier snowmelt in August. The MCA analysis showed that
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Relationships between the indicated median monthly anomalies using partial correlation analysis accounting for solar radiation and air temperature anomalies (retaining site as the unit of variation). Given that the interaction term between “month” and snowmelt timing was significant, we included the correlation coefficients and P of the regressions for each of the indicated months separately in each panel (also included in Table 1). Negative values indicate CO₂ uptake and positive values CO₂ release into the atmosphere, and shaded areas are 95% confidence intervals.</figcaption></figure>

<table><caption>**Table 1.** Significance (P) and Pearson's correlation coefficient (r) of the relationships between the indicated monthly median standardized anomalies for June, July, and August retaining site as a unit of variation. Site was retained as the unit of variation by estimating the standardized anomalies by month and site for each of the indicated variables. The anomalies of the indicated variables were regressed with snow depth anomalies using a partial correlation accounting for the anomalies of solar radiation and air temperature, as shown in Fig. 2. The P and r were only included when the P < 0.1 (given that for P > 0.1 we assumed that r is not significantly different from zero), N = 284.</caption><thead><tr><th>Regression model</th><th>Month</th><th>P</th><th>r</th></tr></thead><tbody><tr><td rowspan="3">NEE ~ snow melt | Rg & air T</td><td>June</td><td>&lt;0.001</td><td>0.42</td></tr><tr><td>July</td><td>0.040</td><td>0.21</td></tr><tr><td>August</td><td>&lt;0.001</td><td>-0.48</td></tr><tr><td rowspan="3">GPP ~ snow melt | Rg & air T</td><td>June</td><td>&lt;0.001</td><td>-0.52</td></tr><tr><td>July</td><td>0.001</td><td>-0.33</td></tr><tr><td>August</td><td>0.0074</td><td>0.27</td></tr><tr><td rowspan="3">ER ~ snow melt | Rg & air T</td><td>June</td><td>&lt;0.001</td><td>-0.38</td></tr><tr><td>July</td><td>&lt;0.001</td><td>-0.34</td></tr><tr><td>August</td><td>-</td><td>-</td></tr></tbody></table>

the anomalies in snowmelt timing had the highest squared covariance fraction (SCF) with the monthly median anomalies of GPP, NEE, and ER in June and July, and the lowest in August over the 2004–2019 period (Fig. 3). A similar result was observed in the linear mixed effect model, which showed a significant relationship between snowmelt date and GPP, and NEE, in all summer months, higher $R_m^2$ between the snowmelt date and GPP in June and July, and no significant relationship between snowmelt date and ER in August (Table S3). In late season, other environmental variables had a higher covariance with the GPP, NEE, and ER anomalies than the snowmelt timing (Fig. 3, Table S3).

Our results are consistent with the discrepancy between the observed increase in the maxNDVI over the last four decades and the time-integrated (TI) NDVI which instead has plateaued in the last two decades and even decreased over the last 10 years in several northern arctic ecosystems³³. TI-NDVI considers the length of the growing season and phenological variations³⁴ and, therefore, better integrates vegetation development during the entire growing season. Moisture was shown to be important for the NDVI trends³³,³⁵. Given the potential water limitation to summer carbon uptake in northern ecosystems¹²,²³⁻²⁵, we tested if an earlier snowmelt was associated with a decrease in soil moisture, which would affect GPP and NEE. We only observed a significant correlation between soil moisture anomalies and snowmelt date anomalies in June (i.e. higher soil moisture with earlier snowmelt, Fig. S1a, Table S2), but no significant correlation in July and August (Fig. S1a, Table S2). The higher soil moisture with earlier snowmelt in June is consistent with surface inundation after snowmelt³⁶,³⁷ and earlier soil thawing resulting in higher soil moisture (i.e., soil moisture is low while soils are frozen). A similar result was observed for the ET anomalies. Higher ET with earlier snowmelt in June (Fig. S1b) could be the result of surface inundation after snowmelt³². The standardized NEE anomalies were significantly correlated with the soil moisture anomalies in each of the summer months (Fig. S1d, Table S2). However, the relationship between
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Squared covariance fraction (SCF) of each couple of the indicated variables for the maximum covariance analysis (MCA) of the monthly median anomalies of GPP, ER, and NEE in June, July, and August. The first pair of singular vectors are the phase-space directions when projected that have the largest possible cross-covariance. The singular vectors describe the patterns in the anomalies that are linearly correlated. A higher SCF indicates a stronger association over time between the indicated variables.</figcaption></figure>

the GPP (and ER anomalies) and soil moisture anomalies was only significant in June (Fig. S1e,f, Table S2) suggesting an earlier activation of the vegetation with earlier soil thaw (and the associated higher soil moisture). A higher water loss from ET in early season (Fig. S1b) could have resulted in the drying of the surface moss layer with the progression of the summer, which would have been consistent with the observed lower GPP and the lower net CO₂ sequestration with earlier snowmelt observed in August (Fig. 2a,b, Table 1). A potential moisture limitation to plant productivity might have been consistent also with the higher SCF of NEE, or GPP and VPD anomalies in August than in June and July (Fig. 3). However, no significant relationship between ET (or soil moisture) and snowmelt date anomalies was observed in July and August (Fig. S1a,b) contrary to what would be expected if drying occurred following earlier snowmelt. No significant relationship was found between VPD anomalies and snowmelt date anomalies in any of the summer months (P = 0.14 in a partial correlation considering air temperature and solar radiation anomalies). Finally, surface drying should result in an increase in the Bowen ratio anomalies with the progression of the summer, given that H increases with a decrease in water table and surface drying³²,³⁸. However, the Bowen ratio showed no correlation with the standardized snowmelt date anomalies in any of the summer months (Fig. S1c, Table S2), and presented similar values in all the summer months (Fig. S2a). The lack of correlation between the soil moisture, VPD, Bowen ratio, and snowmelt date anomalies suggests that an earlier snowmelt did not result in significant surface drying in the sites of this study. The median PET-ET and PPT-ET for all years and sites included in this analysis (Fig. S2 b,c) was slightly higher in August, similar to reports by others for the Russian arctic tundra³⁸,³⁹, further supporting a lack of soil moisture limitation in late season. Although these analyses do not consider runoff, which can be significant²¹,²⁶,
overall our results do not suggest that an earlier snowmelt resulted in a water stress that significantly limited plant productivity in these arctic ecosystems over continuous permafrost.

The correlation between the anomalies in the August GPP and snowmelt timing is consistent with earlier senescence in northern plant species (e.g. *Eriophorum vaginatum*, a dominant species across these tundra types) compared to southern species growing in the same location in a common garden experiment⁴⁰. The phenotypic variation was shown to be persisting for decades⁴¹, and ecotypes may be unable to extend their effective growth period or take advantage of a longer growing season⁴⁰. Several studies across different plant functional types have shown that once plant growth is initiated after the snowmelt in northern ecosystems, it continues only for a fixed number of days until the occurrence of senescence⁴²⁻⁴⁴. Therefore, the lower GPP in August with earlier snowmelt might not be linked to water limitation on photosynthesis later in the season, but rather to an earlier senescence arising from the endogenous rhythms of growth and senescence, that plant functional types living in these extreme conditions have developed over decades. On a broader scale, earlier senescence with an earlier start of the growing season after snowmelt in northern ecosystems is also consistent with an earlier spring zero-crossing date and an earlier autumn zero-crossing date of the mean detrended seasonal CO₂ variations at Barrow, AK, USA (NOAA ESRL: https://www.esrl.noaa.gov/gmd/ccgg/obspack/) during 2013–2017 compared to 1980–1984⁵. The spring and autumn zero-crossing date is the time when the detrended seasonal CO₂ variations intersect the zero line in spring and autumn respectively and can be used as an indicator for the start and end of the net CO₂ uptake by vegetation⁴⁵,⁴⁶. On the other hand, NDVI measurements show both an earlier start of the season and a later end of season for 2008–2012 compared to 1982–1986⁵. The disagreement between the detrended seasonal atmospheric CO₂ concentration showing an earlier autumn zero-crossing date and the NDVI measurements showing a later end of the season has been explained by the increase in respiration in the fall²⁰. The disagreement between atmospheric CO₂ concentration trends (showing an earlier autumn zero-crossing date), and NDVI (showing a later end of the season,⁵ may also be explained by the challenges in using NDVI as a proxy for plant productivity in these arctic systems. The relationship between NDVI and CO₂ flux and plant productivity is highly variable and non-linear in arctic ecosystems⁴⁷. While some arctic ecosystems have shown that NDVI was strongly correlated with GPP (explaining 75% of the variation in GPP⁴⁸, other studies showed that NDVI was either not significantly correlated with GPP and NEE⁴⁹ or was only able to explain a minor fraction (maximum of 25%) of the variation in NEE and GPP in some arctic tundra ecosystems after accounting for the seasonal variation⁵⁰,⁵¹.

In conclusion, earlier snowmelt was associated with greater net CO₂ uptake and higher GPP in early and peak seasons, but with less net CO₂ uptake and lower GPP later in the summer, in the studied arctic tundra ecosystems. We did not find evidence of a late-season water limitation to GPP with earlier snowmelt. Although several hypotheses can be forwarded to explain the link between snowmelt and late season declines in plant productivity and carbon uptake, the current literature does not provide a definitive explanation (schematic Fig. 4). Future studies should investigate the potential interaction of different processes explaining the response of the carbon dynamics in the Arctic to earlier snowmelt and reconstruct the temporal changes in the carbon balance from these systems. The link between the long-term changes in the CO₂ fluxes and NDVI in tundra ecosystems needs closer examination. Studies should investigate if higher NDVI is definitively associated with higher net CO₂ uptake. Greening of the Arctic might not necessarily translate into more net CO₂ uptake, as early and peak season carbon gains might be offset by a late-season CO₂ loss, and respiration might counterbalance the increase in plant productivity. A better understanding of the processes driving these temporal changes is a fundamental step in advancing our prediction of the response of the arctic CO₂ balance to changing climate.

## Materials and methods

**Site description.** A total of 11 eddy covariance flux tower sites across the Arctic were used in this study, where each site had at least six summers of flux data available (SI, Table S1). Ecosystem-scale CO₂ fluxes were estimated using the eddy covariance method⁵²⁻⁵⁴. Details pertaining to the sites, data processing, and gap-filling are provided in the SI Appendix, Table S1. All sites are located in continuous permafrost tundra regions. Vegetation at the tower sites, the instruments used to measure fluxes, the average environmental conditions at each site, the datasets used in this study for each site, and the references describing the sites are indicated in SI Appendix, Table S1. As shown in Fig. S2, the Bowen ratio reported for this study showed similar values to what previously reported during the growing season months in the Arctic (from 3.9 in a dry heath to 1.6 in a wet fen in Greenland³¹; 0.83³⁸ to 0.20–0.25 in two Siberian Arctic sites³⁹; and 0.51–1.69 in a moist-tussock tundra in Alaska³².

To estimate the standardized anomalies in the soil moisture we selected the most consistent depths and sensors (i.e., the same sensor available for the entire time period in each site, or sensors at the most similar depths in each site and across sites when data from the same sensor was not available due to instrument failure). The number of sensors and the soil depths in each of the sites used for all the analyses were: (CA-DL1: N=2, one in a wet location and one in a dry location (both at -10 cm depth); US-Atq: N = 1 (2010–2019) (-10 cm depth); US-Ivo: N=4 at -5 cm depth; US-Bes: N = 2 (diagonally inserted at 0–10 cm); RU-Che: N = 2 (-8 cm and -16 cm depth); RU-Sam: N = 5 (at -5, -14, in rims at -5, -12, -15 cm depths in the center of ice-wedge polygons); US-ICs: N = 2 (at -2.5 cm depth); GL-ZaH: N = 2 (2000–2004 vertical 0–6 cm and from 2005 onward are at two horizontal depths: -5 cm, -10 cm); CA-TVC: included one sensor inserted horizontally at -20 cm depth. More details on the temporal coverage of the soil moisture data from each site are included in Supplementary Table S1.

The R package 'Evapotranspiration' (Version 1.15⁵⁵) was used to estimate the daily aggregated Priestley-Taylor potential evaporation⁵⁶,⁵⁷ at each of the study sites, then summed into monthly totals. We subtracted monthly total actual ET measured with eddy covariance at the respective sites to estimate the PET-ET shown in Fig. S2b. Raster files of monthly precipitation accumulation were acquired for the months of June–August from
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Schematic of the effect of earlier snowmelt on NEE, GPP, and ER at different times of the season. Earlier snowmelt results in an earlier activation of the vegetation, higher plant productivity, and higher net carbon uptake in June and July. This earlier activation could result in more carbon loss and lower plant productivity with earlier snowmelt in August, potentially related to either environmental stress, or to earlier senescence. Photo credit: Donatella Zona.</figcaption></figure>

TerraClimate⁵⁸ over the years 1959–2019. Precipitation data was then extracted for the Eddy Covariance tower coordinates using the terra package⁵⁹ in R⁶⁰ to estimate monthly total precipitation for June, July, and August for each of the sites (Tables S1). We did not use the precipitation collected by the meteorological sensors installed in the tower sites given large gaps in the site-level dataset. The other environmental variables used in this study were collected at the tower sites. The median difference between the total precipitation and the total ET in each site was estimated to evaluate the PPT-ET at each study site during each time of the season (June, July, and August as shown in Fig. S2c). The median was used as it is less affected by outliers.

**Statistical analysis.** Site-level data. For the analyses performed in this study, we separated the data into different times of the season (early season: June, peak season: July, late season: August), given that some of the environmental controls could be very different given the distinct stages of vegetation development. A partial correlation analysis was carried out to identify the correlation between the monthly median standardized anomalies (the ratio between the anomalies and the climatological standard deviation) of NEE, GPP, ER, ET, snowmelt date, and other environmental variables (most of which covary). The standardized monthly median anomalies were estimated by calculating the monthly median of the standardized weekly anomalies. The standardized weekly anomalies were estimated by subtracting the observed values (for each week, year, and site) from the respective average values (for each week, and site and for the entire period available in each site), and then standardized by dividing the anomalies by the standard deviation for each week, and site for the entire period available in each site. The standardized median monthly anomalies were used instead of weekly values, given that the median is less affected by outliers. The monthly scale was used for all the variables included in these analyses as the most appropriate temporal scale to identify the importance of the variability in soil moisture on CO₂ fluxes (given that soil moisture does not change much at the hourly and weekly scale at these tundra sites).

The NEE, GPP, ER data used in these analyses were gap-filled using standard methodologies as described in the Supplementary Information. The partial correlations tested the relationships between the standardized snowmelt date anomalies and the monthly median standardized anomalies of NEE, GPP, ER, ET, VPD, and soil moisture, retaining site as the unit of variation (while controlling for solar radiation and air temperature anomalies as the main controls on carbon fluxes). We also tested if the inclusion of a given site within a linear mixed model changed the results of the correlation analysis between the anomalies. To this purpose, linear mixed effects models (nlme package in R⁶⁰) were used to test the significance of the correlation between the above-mentioned anomalies, by including “site” as categorical random effects to account for pseudo-replication due to the different sites measured in different years. Model performance was evaluated based on the Akaike information criterion (AIC) values, on the marginal coefficient of determination ($R_m^2$ similar to the explanatory power of the linear
models) for generalized mixed-effects models as output by the “r.squaredGLMM” function within the “MuMln” package in R<sup>61,62</sup>. These analyses are included in Table 1, Table S2, and Table S3.

To maximize the dataset for each analysis we included all available time periods for the variables regressed in Fig. 2, but only selected the 2004–2019 period for the Maximum Covariance Analysis (MCA) to include a time period where most sites had data available. The MCA was performed on two fields (i.e., anomalies in NEE, GPP, and ER and anomalies in the environmental drivers, see Fig. 3); the columns of the two fields are spatial locations (each site was retained as a unit of variation of this analysis) and rows are temporal measurements. Site was retained as the unit of variation by estimating the median standardized anomalies by year, week, and site for each of the indicated variables, and then by estimating a monthly median from these weekly anomalies. The first pair of singular vectors are the phase-space directions when projected that have the largest possible cross-covariance. The singular vectors describe the patterns in the anomalies that are linearly correlated. We used the time series of the first singular value decomposition (SVD) mode to visualize the parts of the datasets that vary together and reported the squared covariance fraction (SCF) of the MCA in Fig. 3. Given the limited length of the dataset we did not discuss the long-term changes in the reported anomalies. However, the MCA allowed us to evaluate the influence of snowmelt timing on the carbon balance over time at different times of the growing season, and compare its relative importance to other variables. All analyses were carried out in R<sup>60</sup>.

## Data availability

The eddy covariance data from RU-Che, RU-Cok, and GL-ZaH (previously named DK-ZaH), CA-DL1, were obtained from the European Fluxes Database (http://www.europe-fluxdata.eu/home), from the Ameriflux Database (http://ameriflux.lbl.gov/), with some updated versions provided directly by the principal investigators of each site (e.g. the data from GL-ZaH are also available on: https://data.g-e-m.dk). The data from US-IC and US-ICs are stored in the http://aon.iab.uaf.edu/data_access. US-Bes, US-Atq, US-Ivo are stored in the Arctic Data Center (Donatella Zona. 2021. Greenhouse gas flux measurements at the zero curtain, North Slope, Alaska, 2012–2021. Arctic Data Center. doi:10.18739/A2ZG6G80B.). The R code for the MCA analysis is available in the webpage: https://edoras.sdsu.edu/~babailey/ (https://edoras.sdsu.edu/~babailey/mcaplot/).

Received: 7 October 2021; Accepted: 14 February 2022

Published online: 21 March 2022

## References

1. Overland, J. E., et al. The NOAA Arctic Report Card, Surface Air Temperature. https://arctic.noaa.gov/Report-Card/Report-Card2019/ArtMID/7916/ArticleID/835/Surface-Air-Temperature. (2019).

2. Liljedahl, A. K. et al. Pan-Arctic ice-wedge degradation in warming permafrost and its influence on tundra hydrology. *Nat. Geosci.* **9**, 312. https://doi.org/10.1038/ngeo2674 (2016).

3. Mudryk, L. R., Kushner, P. J., Derksen, C. & Thackeray, C. Snow cover response to temperature in observational and climate model ensembles. *Geophys. Res. Lett.* **44**, 919–926. https://doi.org/10.1002/2016GL071789 (2017).

4. Mudryk, L., Brown, R., Derksen C., Luojuks, K., Decharme B., & Helfrich S. Terrestrial snow cover. In *Arctic Report Card 2019*. (Richter-Menge, J., Druckenmiller, M. L., Jeffries, M. Eds.). https://www.arctic.noaa.gov/Report-Card. (2019).

5. Piao, S. et al. Characteristics, drivers and feedbacks of global greening. *Nat. Rev. Earth Environ.* **1**, 14–27 (2020).

6. Lucht, W. et al. Climatic control of the high- latitude vegetation greening trend and Pinatubo effect. *Science* **296**, 1687–1689 (2002).

7. Berner, L. T. et al. Summer warming explains widespread but not uniform greening in the Arctic tundra biome. *Nat. Commun.* **11**, 4621. https://doi.org/10.1038/s41467-020-18479-5 (2020).

8. Myers-Smith, I. H. et al. Complexity revealed in the greening of the Arctic. *Nat. Clim. Change* **10**, 106–117. https://doi.org/10.1038/s41558-019-0688-1 (2020).

9. Forbes, B. C., Fauria, M. M. & Zetterberg, P. Russian Arctic warming and 'greening' are closely tracked to tundra shrub willows. *Glob. Change Biol.* **16**, 1542–1554. https://doi.org/10.1111/j.1365-2486.2009.02047.x (2010).

10. Lara, M. J., Nitze, I., Grosse, G., Martin, P. & McGuire, A. D. Reduced arctic tundra productivity linked with landform and climate change interactions. *Sci. Rep.* **8**, 2345 (2018).

11. Miles, V. V. & Esau, I. Spatial heterogeneity of greening and browning between and within bioclimatic zones in northern West Siberia. *Environ. Res. Lett.* **11**, 115002 (2016).

12. Gonsamo, A., Ter-Mikaelian, M. T., Chen, J. M. & Chen, J. Does earlier and increased spring plant growth lead to reduced summer soil moisture and plant growth on landscapes typical of tundra–taiga interface? *Remote Sens.* **11**, 1989 (2019).

13. Gamm, C. M. et al. Declining growth of deciduous shrubs in the warming climate of continental western Greenland. *J. Ecol.* **106**, 640–654. https://doi.org/10.1111/1365-2745.12882 (2018).

14. Bruhwiler, L., Parmentier, F.-J.W., Crill, P., Leonard, M. & Palmer, P. I. The Arctic carbon cycle and its response to changing climate. *Curr. Clim. Change Rep.* **7**, 14–34. https://doi.org/10.1007/s40641-020-00169-5 (2021).

15. Humphreys, E. R. & Lafleur, P. M. Does earlier snowmelt lead to greater CO₂ sequestration in two low Arctic tundra ecosystems? *Geophys. Res. Lett.* https://doi.org/10.1029/2011GL047339 (2011).

16. Parmentier, F. J. W. et al. Spatial and temporal dynamics in eddy covariance observations of methane fluxes at a tundra site in northeastern Siberia. *J. Geophys. Res. Biogeosci.* **116**, G03016 (2011).

17. Lund, M. et al. Trends in CO₂ exchange in a high Arctic tundra heath, 2000–2010. *J. Geophys. Res. Biogeosci.* https://doi.org/10.1029/2011jg001901 (2012).

18. Ueyama, M. et al. Growing season and spatial variations of carbon fluxes of Arctic and boreal ecosystems in Alaska (USA). *Ecol. Appl.* **23**(8), 1798–1816. https://doi.org/10.1890/11-0875.1 (2013).

19. López-Blanco, E. et al. Multi-year data-model evaluation reveals the importance of nutrient availability over climate in arctic ecosystem C dynamics. *Environ. Res. Lett.* **15**, 094007 (2020).

20. Piao, S. et al. Net carbon dioxide losses of northern ecosystems in response to autumn warming. *Nature* **451**, 49–52 (2008).

21. Liljedahl, A. K. et al. Tundra water budget and implications of precipitation underestimation. *Water Resour. Res.* **53**, 6472–6486 (2017).

22. Park, T. et al. Changes in growing season duration and productivity of northern vegetation inferred from long-term remote sensing data. *Environ. Res. Lett.* **11**, 084001 (2016).

23. Parida, B. R. & Buermann, W. Increasing summer drying in North American ecosystems in response to longer nonfrozen periods. *Geophys. Res. Lett.* **41**, 5476–5483 (2014).
24. Angert, A. et al. Drier summers cancel out the CO₂ uptake enhancement induced by warmer springs. *Proc. Natl. Acad. Sci. U.S.A.* **102**, 10823–10827 (2005).

25. Buermann, W. et al. Widespread seasonal compensation effects of spring warming on northern plant productivity. *Nature* **562**, 110–114 (2018).

26. Lian, X. et al. Summer soil drying exacerbated by earlier spring greening of northern vegetation. *Sci. Adv.* **6**, eaax0255 (2020).

27. Rouse, W. R. The energy and water balance of high-latitude wetlands: Controls and extrapolation. *Glob. Change Biol.* **6**, 59–68 (2000).

28. Zhang, X. et al. On the variation of regional CO₂ exchange over temperate and boreal North America. *Global Biogeochem. Cycles* **27**, 2012GB004383 (2013).

29. Christensen, T. R. et al. Multiple ecosystem effects of extreme weather events in the Arctic. *Ecosystems* **24**, 122–136 (2021).

30. Walker, D. A. et al. The circumpolar Arctic vegetation map. *J. Veg. Sci.* **16**(3), 267–282. https://doi.org/10.1111/j.1654-1103.2005.tb02365.x (2005).

31. Stiegler, C., Lund, M., Christensen, T. R., Mastepanov, M. & Lindroth, A. Two years with extreme and little snowfall: Effects on energy partitioning and surface energy exchange in a high-Arctic tundra ecosystem. *Cryosphere* **10**, 1395–1413. https://doi.org/10.5194/tc-10-1395-2016 (2016).

32. Vourlitis, G. L. & Oechel, W. C. Landscape-scale CO₂ H₂O vapour and energy flux of moist-wet coastal tundra ecosystems over two growing seasons. *J. Ecol.* **85**, 575–590 (1997).

33. Bhatt, U. S. et al. Climate drivers of Arctic tundra variability and change using an indicators framework. *Environ. Res. Lett.* **16**, 055019 (2021).

34. Tucker, C. J. & Sellers, P. J. Satellite remote sensing of primary production. *Int. J. Remote Sens.* **7**, 1395–1416 (1986).

35. Arndt, K. A. et al. Arctic greening associated with lengthening growing seasons in Northern Alaska. *Environ. Res. Lett.* **14**, 125018 (2019).

36. Bowling, L. C., Kane, D. L., Gieck, R. E., Hinzman, L. D. & Lettenmaier, D. P. The role of surface storage in a low-gradient Arctic watershed. *Water Resour. Res.* **39**, 4 (2003).

37. Woo, M.-K., Young, K. L. & Brown, L. High Arctic patchy wetlands: Hydrologic variability and their sustainability. *Phys. Geogr.* **27**, 297–307 (2006).

38. Goekede, M. et al. Shifted energy fluxes, increased Bowen ratios, and reduced thaw depths linked with drainage-induced changes in permafrost ecosystem structure. *Cryosphere* **11**, 2975–2996 (2017).

39. Runkle, B. R. K., Wille, C., Gažović, M., Wilmking, M. & Kutzbach, L. The surface energy balance and its drivers in a boreal peatland fen of northwestern Russia. *J. Hydrol.* **511**, 359–373 (2014).

40. Parker, T. C., Tang, J., Clark, M. B., Moody, M. M. & Fetcher, N. Ecotypic differences in the phenology of the tundra species *Eriophorum vaginatum* reflect sites of origin. *Ecol. Evol.* **7**, 9775–9786 (2017).

41. Souther, S., Fetcher, N., Fowler, Z., Shaver, G. R. & McGraw, J. B. Ecotypic differentiation in photosynthesis and growth of *Eriophorum vaginatum* along a latitudinal gradient in the Arctic tundra. *Botany* **92**, 551–561 (2014).

42. Bjorkman, A. D., Elmendorf, S. C., Beamish, A. L., Vellend, M. & Henry, G. H. R. Contrasting effects of warming and increased snowfall on Arctic tundra plant phenology over the past two decades. *Global Change Biol.* **21**, 4651–4661 (2015).

43. Rosa, R. K. et al. Plant phenological responses to a long-term experimental extension of growing season and soil warming in the tussock tundra of Alaska. *Glob. Change Biol.* **21**, 4520–4532 (2015).

44. Semenchuk, P. R. et al. High Arctic plant phenology is determined by snowmelt patterns but duration of phenological periods is fixed: an example of periodicity. *Environ. Res. Lett.* **11**, 125006 (2016).

45. Keeling, C. D., Chin, J. F. S. & Whorf, T. P. Increased activity of northern vegetation inferred from atmospheric CO₂ measurements. *Nature* **382**, 146–149 (1996).

46. Piao, S. et al. Weakening temperature control on the interannual variations of spring carbon uptake across northern lands. *Nat. Clim. Change* **7**, 359–363 (2017).

47. Beamish, A. et al. Recent trends and remaining challenges for optical remote sensing of Arctic tundra vegetation: A review and outlook. *Remote Sens. Environ.* **246**, 111872 (2020).

48. Street, L. E., Shaver, G. R., Williams, M. & Van Wijk, M. T. What is the relationship between changes in canopy leaf area and changes in photosynthetic CO₂ flux in arctic ecosystems? *J. Ecol.* **95**, 139–150 (2007).

49. Zona, D. et al. Characterization of the carbon fluxes of a vegetated drained lake basin chronosequence on the Alaskan Arctic Coastal Plain. *Glob. Change Biol.* **16**, 1870–1882 (2010).

50. La Puma, I. P., Philippi, T. E. & Oberbauer, S. F. Relating NDVI to ecosystem CO₂ exchange patterns in response to season length and soil warming manipulations in arctic Alaska. *Remote Sens. Environ.* **109**, 225–236. https://doi.org/10.1016/j.rse.2007.01.001 (2007).

51. Olivas, P. C. et al. Effects of fine-scale topography on CO₂ flux components of Alaskan coastal plain tundra: Response to contrasting growing seasons. *Arct. Antarct. Alp. Res.* **43**, 256–266 (2011).

52. Burba, G. G., McDermitt, D. K., Grelle, A., Anderson, D. J. & Xu, L. Addressing the influence of instrument surface heat exchange on the measurements of CO₂ flux from open-path gas analyzers. *Glob. Change Biol.* **14**(8), 1854–1876 (2008).

53. Burba, G. et al. Calculating CO₂ and H₂O eddy covariance fluxes from an enclosed gas analyzer using an instantaneous mixing ratio. *Glob. Change Biol.* **18**(1), 385–399 (2012).

54. Burba, G. Eddy covariance method for scientific, industrial, agricultural and regulatory applications: A field book on measuring ecosystem gas exchange and areal emission rates. in *LI-COR Biosciences*. 331. ISBN: 978-0-61576827-4. (2013).

55. Guo, D., Westra, S. & Maier, H. R. An R package for modelling actual, potential and reference evapotranspiration. *Environ. Modell. Softw.* **78**, 216–224 (2016).

56. McMahon, T. A., Peel, M. C., Lowe, L., Srikanthan, R. & McVicar, T. R. Estimating actual, potential, reference crop and pan evaporation using standard meteorological data: A pragmatic synthesis. *Hydrol. Earth Syst.* **17**, 1331–1363 (2013).

57. Priestley, C. & Taylor, R. On the assessment of surface heat flux and evaporation using largescale parameters. *Mon. Weather Rev.* **100**(2), 81–92 (1972).

58. Abatzoglou, J. T., Dobrowski, S. Z., Parks, S. A. & Hegewisch, K. C. TerraClimate, a high-resolution global dataset of monthly climate and climatic water balance from 1958–2015. *Sci. Data* **5**, 170191 (2018).

59. Hijmans, R. J. terra: Spatial Data Analysis. R Package Version 1.1-4. https://CRAN.R-project.org/package=terra (2021).

60. R Core Development Team. R: A Language and Environment for Statistical Computing v. 4.0.2. (R Foundation for Statistical Computing, 2020).

61. Johnson, P. C. D. Extension of Nakagawa & Schielzeth's R2GLMM to random slopes models. *Methods Ecol. Evol.* **5**, 944–946 (2014).

62. Nakagawa, S. & Schielzeth, H. A general and simple method for obtaining R² from generalized linear mixed-effects models. *Methods Ecol. Evol.* **4**, 133–142 (2013).

## Acknowledgements

The complete list of funding bodies that supported this study is included in the SI Appendix. DZ, WCO, XX, and DAL acknowledge support from the Office of Polar Programs of the National Science Foundation (NSF) (award number 1204263, and 1702797) with additional logistical support funded by the NSF Office of Polar Programs,
from the NASA Carbon in Arctic Reservoirs Vulnerability Experiment (CARVE), an Earth Ventures (EV-1) investigation, under contract with the National Aeronautics and Space Administration, and from the NASA ABoVE (NNX15AT74A; NNX16AF94A) Program. JDW acknowledges support from NASA NNH17ZDA001N-NIP. The Alaskan sites are located on land owned by the Ukpeagvik Inupiat Corporation (UIC). This project has received funding from the European Union's Horizon 2020 research and innovation program under grant agreement No. 727890, and by the Natural Environment Research Council (NERC) UAMS Grant (NE/P002552/1), and from the NOAA Cooperative Science Center for Earth System Sciences and Remote Sensing Technologies (NOAA-CESSRST) under the Cooperative Agreement Grant # NA16SEC4810008. Part of the analysis was carried out in part at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration (80NM0018D0004).

**Author contributions**

D.Z. and P.L. wrote the main manuscript text and D.Z. performed most of the analyses and prepared most of the figures, K.A. helped to prepare Fig. 2, E.L.-B. helped to prepare Fig. 3. B.B. performed the MCA analysis and wrote the MCA R code. All authors reviewed the manuscript and provided feedback to the analysis and data interpretation.

**Competing interests**

The authors declare no competing interests.

**Additional information**

**Supplementary Information** The online version contains supplementary material available at https://doi.org/10.1038/s41598-022-07561-1.

**Correspondence** and requests for materials should be addressed to D.Z.

**Reprints and permissions information** is available at www.nature.com/reprints.

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_8.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2022