in

Cryogenic land surface processes shape vegetation biomass patterns in northern European tundra

[adace-ad id="91168"]

Study area

The study area (78 000 km2) is located between 68–71°N and 20–26°E, with strong climatic gradients, ranging from wet maritime to relatively dry continental, over tens of kilometers. The landscape of this climatically sensitive high-latitude region has been affected by multiple glaciations in the past. It includes the Scandes Mountains near the Arctic Ocean and low-relief areas to the south and east. The majority of the region (52%) is underlain by sporadic permafrost. Continuous and discontinuous permafrost are limited to the highest mountains of the study area (2% and 7%, respectively)17,26. This large proportion of sporadic, typically warm and shallow permafrost in the study area indicates that ground thermal response to climate warming can be rapid27. Our data do not cover low-relief plateaus of continuous permafrost (similar to northern Siberia and Alaska), where the generally high ice content of soil may lead to different and enhanced LSP responses under climate warming (e.g., ice wedge degradation and surface ponding) with altered AGB feedbacks43,44.

LSP observations

The data consist of 2917 study sites (each 25 m × 25 m) and includes previously combined observations (both in-situ [n = 581] and remote-sensing [n = 2336]) of the active surface features of three cryogenic LSP common in the area: cryoturbation, solifluction, and nivation. These LSP are mainly associated with seasonal freeze–thaw processes. Cryoturbation (i.e., frost churning) is a general term for soil movement caused by differential heave, and it creates typical surface features such as patterned ground, frost boils and hummocky terrain5. Solifluction is the slow mass wasting of surficial deposits through frost creep and permafrost flow, where gravitation causes frost-heaved soil to settle downwards during the summer thaw, creating features of lobes and terraces50. In addition, solifluction also includes gelifluction which is a mass wasting process caused by high porewater pressure in unconsolidated surface debris creating similar lobes and terraces5,50. We use the term nivation to collectively designate various weathering and fluvial processes which are intensified and depicted by the presence of snowbeds (which in general are melting in mid-July – late-August) and nivation hollows28,51. We expect the presence of such a snowbed to be an indication of active nivation processes, since in these environments the year-to-year spatial snow patters are fairly consistent31.

The rationale behind LSP sampling is described in previous geomorphic studies which served as a basis for the used protocol52,53. Due to the large study domain, study objectives (focus on distribution of active surface features, not on activity itself) and modeling resolution (50 m × 50 m), we used a visual method to estimate the presence/absence of the mapped LSP. We used high-resolution aerial photography (spatial resolution of 0.25 m−2) and targeted field surveys (GPS accuracy ~5 m; Garmin eTrex personal navigator) to construct the LSP dataset. A binary variable (1 = presence, 0 = absence) was assigned to each LSP to indicate their evident activity (or absence). The activity/absence of the LSP was visually estimated based on the evidence in ground surface, indicated by e.g., frost-heaving, cracking, microtopography (e.g., erosional and depositional forms), soil displacement indicative to a process form (e.g., solifluction lobes, patterned ground), changes in vegetation cover and late-lying snow. Such indicators average the LSP activity over several years. Even small areas with slight indication of activity were considered active processes. However, such a protocol based on a visual assessment is susceptible for incorrect activity classification; solifluction may be active despite having a complete vegetation cover19 and the presence of late-lying snowbed, although being a good indication28, does not necessarily mean that active nivation processes are present.

Remotely sensed vegetation index

For obtaining remotely sensed vegetation index for the study area, we employed a maximum-value compositing approach. We downloaded all available clear sky (less than 80% land cloud cover) Landsat OLI 8 images overlapping the study area from June to September between 2013 and 2017 (total of 1086 scenes) from the United States Geological Survey (USGS) database (http:earthexplorer.usgs.gov). Images were USGS surface reflectance products, which were preprocessed (georeferencing, projection, and atmospheric corrections) by USGS54. Landsat-8 satellite is the latest addition to the Landsat mission that has provided repeated land surface information globally since the 1970’s and is the most commonly used fine-scale satellite system for vegetation mapping. The native resolution of the Landsat OLI sensor is 30 m for the spectral bands used in the image processing steps of this study.

Normalized difference vegetation index (NDVI), a widely used spectral index to estimate the amount of green vegetation, was calculated as55:

$$({{{{{rm{rho }}}}}}{{{{{rm{NIR}}}}}}-{{{{{rm{rho }}}}}}{{{{{rm{red}}}}}})/({{{{{rm{rho }}}}}}{{{{{rm{NIR}}}}}}+{{{{{rm{rho }}}}}}{{{{{rm{red}}}}}})$$

(1)

where ρNIR and ρred are the surface reflectance for their respective Landsat bands, 0.851–0.879 (mu)m and 0.636–0.673 (mu)m.

USGS provides pixel-based quality assessment bands for all surface reflectance products. These bands were used to mask clouds, snow, water, and other low-quality pixels from the individual NDVI scenes. Additionally, if the NDVI images still had unphysical values over 1 or under -1, these pixels and their surroundings of 100 m radius were excluded. We determined maximum values for each 30 m resolution pixel of the study area individually. After masking cloud, snow, and water from the scenes, obvious scattered erroneous NDVI values remained in some scenes. Therefore, we excluded the values outside the pixel-based 95% percentile prior to maximum composite.

The CFmask cloud detection algorithm that is used to generate the quality assessment band has clear difficulties in distinguishing small snow patches from clouds. As such, a large portion of late-lying snow beds were repeatedly and incorrectly classified as clouds. Moreover, the CFmask algorithm creates buffers around the cloud pixels54, hence much information was lost around the snow patches that were incorrectly identified as clouds. After these processing steps, some pixels around the extreme late-lying snow beds had still too low number of NDVI records to provide reliable NDVI values for the maximum composite. To fill these small and scattered gaps in the initial maximum NDVI composite, we selected 74 mostly cloud-free scenes between August and September. For these 74 scenes, we manually digitized cloud masks to exclude cloud-contaminated pixels with high certainty. Moreover, every pixel must have passed the following quality checks to be included in the gap-filling composite: not classified as water in the USGS quality assessment band; normalized difference snow index (NDSI) value less than 0.4, and blue band reflectance less than 0.1 (to exclude snow); reflectance of red band between 0.03 and 0.4 (second check for water and snow, and deepest shadows); NDVI between 0 and 0.4 (lower threshold to exclude snow and water contamination; higher threshold to exclude erroneous values, as very late snowbed habitats always have very limited vegetation cover). Additionally, if the NDVI images had unphysical values over 1 or under -1, these pixels and their surroundings (200 m radius) were excluded. Pixels in the 74 selected images which passed these checks, were then used to create a secondary maximum NDVI composite that was used to fill the gaps in the initial maximum NDVI composite. The secondary composite comprised 0.4% of the pixels in the final composite. Among all 2917 LSP observation sites, 2.9% were located within the gaps in the initial maximum NDVI composite, and thus received their maximum NDVI values from the secondary NDVI composite.

In the used Landsat data, the nivation sites were not covered by snow, but instead were associated with generally lower AGB values as nival processes affect the vegetation’s structure and composition (Supplementary Table 1).

Above-ground biomass data

Above-ground biomass (AGB) reference data were collected from two regions, with a total of 433 sites that represent an area of > 4000 km2 (Supplementary Fig. 9). The first dataset (hereafter BM region 1; centering to ca. 69°N, 21°E) was collected between 2008 and 2011, and the second dataset (BM region 2; centering to ca. 70°N, 26.2°E) between 2015 and 2017. Both study regions are representative of an arctic and alpine treeline ecotone and include data from mountain birch forest to barren oroarctic tundra56,57.

The BM region 1 dataset consists of 309 field sites (each 10 m × 10 m), which are located around eight different massifs covering a wide range of environmental conditions (Supplementary Figs. 9–10). Sampling was performed in transects to cover various aspects of the slope (i.e., topoclimatic conditions), starting from the foothill of the mountain, and ending at the summit. A plot was systematically established at every 20 m increase in elevation and recorded with a GPS device. Four clip-harvest biomass samples (20 cm × 20 cm) were taken 5 m from the plot center in every cardinal direction. Two samples were used in bare mountaintops (north, south). The clip-harvest samples were dried for 48 h at +65 °C, and dry weight was recorded. The sample biomass values were converted to g m-2 and the average sample value was calculated for each site (Supplementary Fig. 9). The original BM region 1 dataset contains forest and treeline plots, but these were excluded from the final analyses due to an incomparable tree sampling strategy with BM region 2, which could introduce uncertainty into biomass estimates.

The BM region 2 data were collected from three different massifs having an elevation range from 120 m to 1064 m (Supplementary Fig. 9). The biomasses were sampled from 102 sites (each 24 m × 24 m in size) that were chosen using a stratified sampling to cover gradients of thermal radiation (potential incoming solar radiation), soil moisture (topographic wetness index, TWI) and vegetation zone (forest, treeline, and alpine zones). Radiation and TWI were calculated from a 10 m digital elevation model (DEM, provided by the National Land Surveys of Finland and Kartverket, the Norwegian mapping authority), and assigned to one of three classes based on observation percentiles (breaks at 20% and 80%) leading to total of 27 strata. Vegetation zones were digitized based on aerial imagery. After the first field survey, 22 sites were added to account for vegetation types that were not sufficiently represented by the GIS-based stratification. Thus, the total sample size of the BM region 2 dataset is 124 AGB sites.

The same clip-harvest sample protocol was used as in BM region 1; additional samples were also taken from 12 m in every cardinal direction, thus each site had eight AGB samples (Supplementary Fig. 9). Trees with diameter at breast height (DBH) greater than 20 mm were measured from a 900 m2 circular plot, which corresponds to the size of the NDVI product resolution. Large stems (DBH > 80 mm in the forest and 40 mm at the treeline) were measured from the whole plot, whereas smaller stems were measured from five subplots. Specifically, the center subplot was 100 m2, and the four subplots located at 8 m to every cardinal direction were each 12.5 m2. For the subplot observations, we used a plot expansion factor (900/150 = 6) to generalize the observations for the whole plot assuming a homogeneous forest structure i.e., each subplot stem represents six trees within the 900 m2 plot. A total of 98% of the measured stems were mountain birch (Betula pubescens ssp. czerepanovii), making it the most abundant species in the area. For predicting stem biomass, we used the average of three allometric equations58,59,60, in order to reduce the uncertainty related to the transferability of an individual allometric model. In addition, Populus tremula (1% of the observations) were found on low-altitude south-facing slopes, and Salix caprea (1%) in moist, nutrient-rich sites. Species-specific models61,62 were used to estimate their respective stem biomasses. Individual pines (Pinus sylvestris) were scattered in the area but were not present in any of the sampled plots.

The plots of above-ground tree biomass were converted to g m−2 and added to the mean clip-harvest AGB to obtain the total vascular plant AGB for each site. The BM region 1 and BM region 2 datasets were combined, and the NDVI value was extracted from the site center coordinates.

Spatial autocorrelation (SAC) is a common property of any spatial dataset and means that observations are related to one another by the geographical distance63. SAC in the model residuals violates the independence assumption commonly required by statistical models and can lead to inflated hypothesis testing and biased model estimates64. To investigate whether the plot-scale AGB data are spatially autocorrelated, we calculated semivariogram which describes the spatial dependency between the observations as a function of distance between the point pairs65. Semivariogram were calculated as:

$${{{{{rm{gamma }}}}}}(h)=frac{1}{2N(h)}mathop{sum }limits_{i=1}^{{N}_{h}}{left(Zleft({s}_{i}right)-Zleft({s}_{i}+hright)right)}^{2}$$

(2)

where N(h) denotes the number of data pairs within distance h, and (Zleft({s}_{i}right)) is an observation (or model residual) in location i. For the calculation, we used R package gstat66 (version 2.0-0). A visual inspection of the semivariogram indicated spatial autocorrelation at short distances (“AGB” in Supplementary Fig. 11). Therefore, for the NDVI-AGB conversion, we used a generalized least squares modeling (GLS, as implemented in R package nlme67 [version 3.1-137]) that can explicitly account for SAC in the data. For the modeling, the AGB values were log(x+0.1) transformed. The GLS, where AGB was modeled as a function of NDVI, were fitted assuming an exponential spatial correlation structure:

$$gamma left(hright)={c}_{0}+cleft(1-{e}^{-h/a}right)$$

(3)

where ({c}_{0}) is the difference between the intercept and origin (i.e., the “nugget” parameter in geostatistics), (c) is the amount of variance (i.e., the “sill”) and a represents the distance of spatial dependency (i.e., the “range”). The fitted GLS was as follows:

$${{log }}left({AGB}right)=-1.038629+9.725572times {NDVI}$$

(4)

The estimated spatial correlation parameters were c0 = 0.516, c = 0.484 and a = 260.605, indicating that the distance of spatial autocorrelation extends to ca. 261 m. The semivariogram for the model residuals indicated a notable reduction in the amount of spatial autocorrelation compared to the AGB data (Supplementary Fig. 11). The fitted model explained 70.6% of the deviance in the data. When the predicted values were converted back to the response scale, the model explained 60.5% of the deviance. Therefore, for the subsequent analyses we use the above-ground biomass estimated by the model.

Environmental predictors

In addition to LSP, we used climate, topography, and soil predictors to model AGB. Gridded monthly average temperatures and precipitation data (1981–2010; spatial resolution 50 m × 50 m) based on a large collection (> 950) of Fennoscandian meteorological stations were used in a spatial interpolation scheme17. Three climate predictors—growing degree days (GDD, °C, base temperature 5 °C), mean February air temperature (Tfeb, °C) and water balance (WAB, mm)—were calculated from the gridded climate data. WAB is the difference between total annual precipitation and potential evapotranspiration (PET), which was estimated from the monthly air temperature and precipitation data68:

$${PET}=58.93times {T}_{{above}0^circ C}/12$$

(5)

These climatic predictors were selected to represent different aspects of climate that are critical for tundra vegetation: heat requirements, cold tolerance and moisture availability. In addition, two local scale topographic predictors were calculated from a DEM (spatial resolution of 50 m × 50 m, provided by the National Land Survey Institutes of Finland, Norway, and Sweden): topographic wetness index69 (TWI, a proxy for soil moisture) and potential annual direct solar radiation70 (MJ cm-2 a-1). Slope angle was initially considered as a potential predictor for AGB but was later omitted due to the strong correlation with TWI (-0.93, P ≤ 0.001). We also calculated peat cover (%) from a digital land cover classification71. Here, the native resolution of 100 m was resampled at 50 m to match the resolution of the climatic and topographic predictors, using nearest-neighbor interpolation. The binary peat cover variable was transformed to a continuous scale using a spatial mean filter of 3 × 3 pixels52. Finally, the topmost soil layer of a global gridded soil database72 was used to obtain pH data. Again, the original resolution of 250 m was also resampled to 50 m resolution using bilinear interpolation.

Our fine-scale data revealed strong environmental gradients over the 78,000 km2 study area (Supplementary Table 1), most of which were only moderately inter-correlated (Spearman’s correlation coefficient < [0.7]; Supplementary Table 3). All predictors showed clear spatial autocorrelation structures (Supplementary Figs. 2–3).

Climate data for future warming scenarios

To investigate tundra AGB change along climate warming, we prepared climate data by modifying monthly air temperatures from +0.5 °C to +8 °C at 0.5 °C intervals. To account for the unequal warming among different months, we calculated monthly adjustment factors (Supplementary Table 4). The relative changes in current monthly air temperatures were compared to those projected by an ensemble of 23 global climate models derived from the archives of Coupled Model Intercomparison Project Phase 5 (CMIP5)73 for the period of 2070–2099, assuming a representative concentration pathway (RCP) 8.574. The monthly precipitation data needed to be adjusted to account for the air temperature increase. The amount of adjustment was determined by regressing the projected change in monthly precipitation ((triangle)Prec, %) by the CMIP5 models against the projected change in the corresponding monthly air temperature ((triangle)T, °C; Supplementary Fig. 12). The resulting relationship (triangle)Prec = −0.445 + 3.85 (times) (triangle)T was used to adjust the monthly precipitation data at each temperature scenario. Subsequently, all climate predictors (i.e., GDD, Tfeb and WAB) were recalculated for each corresponding warming scenario.

Future LSP data

We used the recently developed statistical modeling approach17 to predict future LSP distributions at the study area. In brief, this was done by relating each investigated LSP to climatic, topographic, and soil predictors using an ensemble of four statistical techniques ranging from regression (generalized linear models, generalized additive models) to machine learning (boosted regression trees, random forest). The LSP predictions of cryoturbation, solifluction and nivation over the study area (spatial resolution 50 m × 50 m) were produced for each corresponding warming scenario (i.e., from +0.5 °C to +8 °C at 0.5 °C intervals).

Statistical AGB modeling

Prior to modeling, the AGB values were log-transformed to approximate a normal distribution. All statistical analyses were conducted using the R statistical computing environment75. We used boosted regression trees (BRT, as implemented in the R package “gbm”76) to link AGB to the climatic, topo-edaphic and LSP predictors using all 2917 observations and assuming Gaussian error distribution. BRT is a sequential machine-learning method that combines a large number of iteratively fitted regression trees to a single model with improved prediction accuracy29. BRTs automatically incorporate interactions between predictors, are capable of modeling highly complex non-linear systems, and provide an estimate of the predictors’ relative influence in the models. We optimized BRTs using an “out-of-bag” estimator and setting the interaction depth to six, learning rate to 0.01 and bag fraction to 0.75. In order to control for potential differences caused by choice of modeling technique, we also ran parallel analyses using random forest30 (RF) with the maximum number of trees set to 500. To assess whether accounting for LSP improves the AGB predictions, the models were fitted using two configurations:

$${AGB}={GDD}+{WAB}+{T}_{{Feb}}+{TWI}+{Peat}+{RAD}+{pH}+{RA}{C}_{{base}}$$

(6)

[Base model]

$${AGB}= ,{GDD}+{WAB}+{T}_{{Feb}}+{TWI}+{Peat}+{RAD}+{pH}+{RA}{C}_{{LSP}} +{{{{{boldsymbol{Cryoturbation}}}}}}{{{{{boldsymbol{+}}}}}}{{{{{boldsymbol{Solifluction}}}}}}{{{{{boldsymbol{+}}}}}}{{{{{boldsymbol{Nivation}}}}}}$$

(7)

[LSP model]

where the RACBase and RACLSP refer to residual autocovariates (RAC). Following the ref. 77, these were calculated from the BRT residuals and then added back to the models as additional predictors to account for SAC. This was done despite a visual inspection of the BRT residuals without the RACs did not reveal evident signs of SAC (Supplementary Fig. 4). However, after calculating semiovariograms over very short distances (≤ 5 km) SAC patterns in the residuals could be detected (Supplementary Fig. 5). The two fitted exponential functions (gamma left(hright)=0.707(1-{e}^{-h/162.385})) for the Base model step and (gamma left(hright)=0.699(1-{e}^{-h/158.249})) for the LSP model step indicated that the SAC extends to ca. 162 m and 158 m, respectively. Consequently, the RACs were calculated using the function autocov_dist() as implemented in R-package spdep78 (version 1.1-3). For calculating the RACBase and RACLSP we used the above-mentioned neighborhood radius of 162 m and 158 m, respectively, with an inverse distance weighting scheme and symmetric weight matrix (style=“B” in the autocov_dist -function, following ref. 79).

The predictive capability of the models was assessed with a repeated cross-validation (CV) scheme, where the models were fitted 500 times at each round using a random sample of 99% of the data (out of 2917, with no replacement) and subsequently evaluated against the remaining 1%. A distance-block (hereafter h-block, i.e., the minimum distance between calibration and evaluation data) of 158 m (panel b in Supplementary Fig. 5) was specified to omit model calibration data located at the immediate vicinity of the evaluation data, which could lead to overly optimistic CV-statistics due to spatial autocorrelation80. The CV procedure yielded an average of 2873 observations available for model calibration, and 30 for evaluation per CV round. After each CV run the predicted and observed AGB were compared in terms of root mean squared error (RMSE) for both model steps.

To quantify the effects of LSP on AGB, the LSP models were used to predict AGB (using all 2917 observations and 500 bootstrap sampling-rounds) under two conditions. The first included that other environmental predictors (i.e., climate, topography, and soil) and RACLSP were fixed to their median values (constrained to an environmental envelope allowing each LSP occurrence), and each LSP is in turn set to absent (0). The second condition was that environmental predictors and RACLSP were fixed to their median values as mentioned above, but each LSP is in turn set to present (1). The total effects of LSP on AGB over the study area were determined by adding the mean effect sizes of each individual LSP as follows:

$${Total},{effect}= ,left({AG}{B}_{{cryot}0}-{AG}{B}_{{cryot}1}right)+left({AG}{B}_{{solif}0}-{AG}{B}_{{solif}1}right) +left({AG}{B}_{{nivat}0}-{AG}{B}_{{nivat}1}right)$$

(8)

Effect sizes for continuous predictors were calculated over the calibration data by comparing the predicted minimum and maximum AGB values while fixing other predictors to their median values.

To obtain spatial estimates, we predicted AGB over the study area at a spatial resolution of 50 m × 50 m, at each temperature scenario based on two model configurations: Base model and LSP model, using all observations. For the predictions, the RACs were fixed to their median values (i.e., zero). This was done in order to avoid additional interpolation over the study area that could increase uncertainty in the consequent analyses. Although the AGB predictions were produced over the entire study area, all comparisons of AGB values between the model steps were limited to areas suitable for the occurrence of LSP (i.e., the LSP realm)17. As the AGB predictions by BRT were in close agreement with those of RF (Spearman’s rank correlation coefficient = 0.96), we only present the BRT results.


Source: Ecology - nature.com

How marsh grass protects shorelines

Influence of historical changes in tropical reef habitat on the diversification of coral reef fishes