Coverage and Rainfall Response of Biological Soil Crusts Using Multi-Temporal Sentinel-2 Data in a Central European Temperate Dry Acid Grassland

Abstract

Biological soil crusts (BSCs) are thin microbiological vegetation layers that naturally develop in unfavorable higher plant conditions (i.e., low precipitation rates and high temperatures) in global drylands. They consist of poikilohydric organisms capable of adjusting their metabolic activities depending on the water availability. However, they, and with them, their ecosystem functions, are endangered by climate change and land-use intensification. Remote sensing (RS)-based studies estimated the BSC cover in global drylands through various multispectral indices, and few of them correlated the BSCs’ activity response to rainfall. However, the allocation of BSCs is not limited to drylands only as there are areas beyond where smaller patches have developed under intense human impact and frequent disturbance. Yet, those areas were not addressed in RS-based studies, raising the question of whether the methods developed in extensive drylands can be transferred easily. Our temperate climate study area, the ‘Lieberoser Heide’ in northeastern Germany, is home to the country’s largest BSC-covered area. We applied a Random Forest (RF) classification model incorporating multispectral Sentinel-2 (S2) data, indices derived from them, and topographic information to spatiotemporally map the BSC cover for the first time in Central Europe. We further monitored the BSC response to rainfall events over a period of around five years (June 2015 to end of December 2020). Therefore, we combined datasets of gridded NDVI as a measure of photosynthetic activity with daily precipitation data and conducted a change detection analysis. With an overall accuracy of 98.9%, our classification proved satisfactory. Detected changes in BSC activity between dry and wet conditions were found to be significant. Our study emphasizes a high transferability of established methods from extensive drylands to BSC-covered areas in the temperate climate. Therefore, we consider our study to provide essential impulses so that RS-based biocrust mapping in the future will be applied beyond the global drylands.

1. Introduction

Biological soil crusts (BSCs), or so-called biocrusts, are thin microbiological layers that “result from an intimate association between soil particles, cyanobacteria, algae, microfungi, lichens, and bryophytes, which live within, or immediately on top of the uppermost millimeters of soil” [1]. The presence or absence of water massively influences the physiognomy of BSCs as they consist of poikilohydric organisms, capable of temporally suspending their metabolic activities, such as photosynthesis, when no water is available [2]. These can be recommenced immediately by precipitation events and kept active for up to several days, depending on water availability [2].

BSCs cover around 12% of the Earth’s land surface [3]. While predominantly occurring within arid and semi-arid drylands, BSCs can be found globally throughout different climate zones [4]. They influence ecosystems by limiting the water availability for soil biota, nutrient cycles, and vascular vegetation [5] and highly contribute to the global carbon and nitrogen budgets [5]. Human society benefits from their ecological functions, such as mitigating water and wind erosion [6] and regulating air quality [7]. BSCs are reported as very sensitive to environmental changes, particularly to rising temperature and decreasing humidity [8]. Therefore, they are highly endangered by climate change and land-use intensification, whereby the global BSC cover is expected to decrease by 25 to 40% from 2005 to 2070 [3]. This also implicates the loss of their ecosystem functions. Acquiring more accurate information on their spatiotemporal distribution and dynamic is considered essential for developing sustainable management and protection strategies [1,9,10].

Beyond the arid and semi-arid drylands, a considerable number of drylands featuring BSCs allocates in temperate climates. In Central Europe, they are mainly associated with natural open ecosystems [11,12,13] or disturbed landscapes [14,15]. Dry grasslands with BSCs and soil lichen communities are widespread on quaternary inland dunes deposited during the last glacial time; for example, they are distributed along a belt stretching from the Netherlands, across northern Germany, including Brandenburg, to Poland, and along rivers in southern Germany [14,15,16,17,18]. These habitats are characterized by nutrient-poor soils, dryness, and high solar irradiation [11], allowing colonization by drought-adapted cryptogams without being outcompeted by vascular plants [13]. For instance, in open xeric grasslands, soil lichen and moss communities accompanied by microalgae are especially frequent and influence the ecosystem development [19]. As photoautotroph microbiotic communities, they contribute to carbon accumulation and soil formation of initial sandy soils [15,20]. Furthermore, they change the hydrological properties of the soil surface and alternate sub-critical water hydrophobicity [17,21].

The first use of remote sensing (RS) for BSC-related research dates back to 1986, when optical imagery from the Landsat TM sensor was analyzed [22]. Ever since, studies employed optical and radar RS data from space-borne, airborne, and Unmanned Airborne Vehicle (UAV)-based sensors in different spatial, radiometric, and temporal resolutions [23,24] (Table 1) for classifying and mapping the distribution of BSCs [3,24,25,26,27,28,29,30,31,32,33]. The spectral features of BSCs were identified and characterized using multi-and hyperspectral data to differentiate them from co-occurring land cover types (e.g., bare soil and vascular plants), to distinguish between different crust types, and to assess differences between dry and wet BSCs [33,34,35,36,37,38,39,40,41]. Additionally, RS data was used to model various processes and activities of BSCs related to their ecosystem functions [2,19,42,43,44,45,46,47,48].

The RS-based assessment of BSCs is predominantly based on the optical part of the electromagnetic spectra (Table 1). BSCs show a unique spectral signature that can be utilized via the calculation of spectral indices [31,36,37,38]. The most common one is the Normalized Difference Vegetation Index (NDVI) [49], capable of delineating different vegetation types and their condition. It estimates the amount of green biomass and chlorophyll using the difference in reflectance in the chlorophyll absorption feature at the red and the high reflection in the near-infrared (NIR) spectral domain [49]. With respect to the BSCs, this means that dry and wet conditions must be considered differently, as the spectral reflectance changes with the change in BSC activity and appearance, based on the presence and absence of water. Dry BSCs slightly differ from bare soil, whereas wet BSCs show stronger chlorophyll absorption and higher NIR reflectance resembling vascular plants [33,34,35]. Therefore, the NDVI was commonly applied for assessing the change in metabolic activity of BSCs after precipitation events, as wet biocrusts are reported to show a significantly higher NDVI than dry biocrusts [37,39,40].

Table 1. Remote sensing sensors used for characterizing and mapping biological soil crusts with the corresponding studies’ timeframes and focus areas. MS = multispectral, HS = hyperspectral, UAV = unmanned airborne vehicle, SAR = synthetic aperture radar. * not specified.

Sensor				Timeframes/Temporal Resolution	Focus Areas	References
Platform	Type	Name	Spatial Resolution
Spaceborne	MS	Kompsat-2	1 m	2 images	Niger	[50]
	MS	QuickBird, WorldView-2	2.5 m	2 images	Negev	[51]
	MS	SPOT	20 m	28 images during 2 years	Negev	[45]
	MS	Landsat 4–8	30–75 m	single images, 31 annual composites	Negev, China, Namibia, South Africa	[22,26,27,28,36,48,52]
	MS	Sentinel-2 MSI	10–20 m	single image, 20 images during 2 years, full time series during 2016	China, Negev, Idaho	[28,30,53]
	MS	Meteosat-9 SEVIRI	3 km	Every 15 min during 2 months	Negev	[46]
	MS	ASTER	90 m	“several” images *	Negev	[54]
	MS	MODIS	250 m–1 km	Every 16 days during growing season, every 8 days from 2000 to 2014	Patagonia, China	[55,56]
	SAR	SIR-C	13–30 m	single image	Negev and Sinai	[57]
	SAR	ASAR	- *	15 images during 1.5 years	Negev	[57]
UAV	MS	Ricoh GR II	0.5 cm	single image	Utah	[58]
	MS	DJI X5s	1–3 cm	single images	Hawaii	[59]
	MS	- *	- *	single images	China	[28]
Airborne	HS	AMS	5 m	single images	Australia	[25]
	HS	CASI 1&2	1–1.5 m	single images	Spain, South Africa	[9,10,29,31]
	HS	DAIS	8 m	single images	Negev	[26,36,54]

Table 1. Remote sensing sensors used for characterizing and mapping biological soil crusts with the corresponding studies’ timeframes and focus areas. MS = multispectral, HS = hyperspectral, UAV = unmanned airborne vehicle, SAR = synthetic aperture radar. * not specified.

Sensor				Timeframes/Temporal Resolution	Focus Areas	References
Platform	Type	Name	Spatial Resolution
Spaceborne	MS	Kompsat-2	1 m	2 images	Niger	[50]
	MS	QuickBird, WorldView-2	2.5 m	2 images	Negev	[51]
	MS	SPOT	20 m	28 images during 2 years	Negev	[45]
	MS	Landsat 4–8	30–75 m	single images, 31 annual composites	Negev, China, Namibia, South Africa	[22,26,27,28,36,48,52]
	MS	Sentinel-2 MSI	10–20 m	single image, 20 images during 2 years, full time series during 2016	China, Negev, Idaho	[28,30,53]
	MS	Meteosat-9 SEVIRI	3 km	Every 15 min during 2 months	Negev	[46]
	MS	ASTER	90 m	“several” images *	Negev	[54]
	MS	MODIS	250 m–1 km	Every 16 days during growing season, every 8 days from 2000 to 2014	Patagonia, China	[55,56]
	SAR	SIR-C	13–30 m	single image	Negev and Sinai	[57]
	SAR	ASAR	- *	15 images during 1.5 years	Negev	[57]
UAV	MS	Ricoh GR II	0.5 cm	single image	Utah	[58]
	MS	DJI X5s	1–3 cm	single images	Hawaii	[59]
	MS	- *	- *	single images	China	[28]
Airborne	HS	AMS	5 m	single images	Australia	[25]
	HS	CASI 1&2	1–1.5 m	single images	Spain, South Africa	[9,10,29,31]
	HS	DAIS	8 m	single images	Negev	[26,36,54]

However, in BSC cover estimation and mapping, the NDVI values can be easily confused with other land cover types, such as soils and sparsely woody vegetation, leading to imprecise results [33]. Therefore, other biocrust-specific indices were developed during the last decades. For example, the Crust Index (CI) [26] uses the red and blue spectral bands to detect BSCs based on an increase in reflectivity within the blue spectral domain that is caused by the presence of cyanobacteria. Later, the green and red band combination was found to be very useful, particularly for distinguishing BSCs from vascular plants [28]. The Biological Soil Crust Index (BSCI) [27] utilizes the green, red, and NIR reflectance. It was specifically designed to detect lichen-dominated BSCs [33]. Additionally, the thermal spectra were studied, building on the principle of variations in land surface temperature originating from differences in albedo and absorbance caused by BSCs [60]. The Thermal Crust Index (TCI) [54] was developed for discriminating BSCs from bare sand in deserts and was reported to be particularly powerful when combined with the NDVI and CI [24].

Besides the RS data, topographic information holds enormous potential for supporting soil and vegetation mapping approaches [61,62,63,64,65,66]. Yet, for classifying BSCs, topographic information is not common and has only been discussed recently. For example, topographic attributes (TAs) such as the slope angle, the wetness index, and the solar irradiation [32] were found to control patterns of BSCs and have been successfully employed as environmental covariates (i.e., predictors) in estimating and mapping BSCs [50].

Land cover classification, in general, can be performed using different approaches. Hereby, the ensemble machine learning (ML) algorithm Random Forest (RF) is known for high classification accuracy and providing meaningful results in the context of RS data [67,68]. It is a supervised technique based on reference areas for training and testing the algorithm by defining the characteristics of the classes. It is widely used in RS-based land cover and vegetation mapping applications and was already successfully used to estimate BSC coverage in desert regions in Xingjiang, China [28].

With the launch of the multispectral optical satellite mission Sentinel-2 (S2) within the ESA Copernicus earth observation program, global data of high spatial and temporal resolution is freely available, allowing for continuous global monitoring. The mission is designed to be comparable to other predecessor missions such as Landsat 8 [69], and was used to estimate BSC coverage and characteristics within large-scale areas in the northern and southern hemispheres [28,30,53].

While previous RS-based analyses of BSCs were predominantly performed in semi-arid and arid drylands (Table 1), investigations in temperate regions focused on small-scale proximal methods using digital camera systems for temporally observing small segments of BSCs [70,71,72].

Our study aims to transfer RS-based methods and indices to the temperate climate regions and marks the first time that BSC coverage is estimated and mapped in Central Europe. As a representative area, we chose a dry acid grassland, where the characteristics of BSCs are well known [70,71,72]. Our objective is to estimate the BSC cover using RF classification incorporating multispectral reflectance bands and indices combined with meaningful topographic information for an extended observation period of more than five years (2015–2020). Furthermore, we aim at identifying the most important parameters (i.e., predictors) in spatiotemporally estimating the BSC coverage to enable better understanding of the BSCs’ sensitivity to them. We aim at subsequentially examining the BSC-covered areas to monitoring the changes in metabolic activity arising from dry and wet cycles caused by rainfall events by employing NDVI time series and daily precipitation data. Our study marks the first time that BSC coverage is estimated and mapped in Central Europe. Simultaneously, it is the first time that an extended dataset of multi-annual time series of BSC related is investigated. Therefore, we consider our study to provide a wide range of impulses for future research concerning biocrusts in temperate climates.

2. Materials and Methods

2.1. Study Area

Our study area (51°56′18″ N, 14°21′54″ E) is an open landscape located approximately 90 km southeast of the German capital Berlin and about 20 km north of the city of Cottbus in the Federal State of Brandenburg (Figure 1). It covers an area of 530 ha, commonly known to locals as ‘Lieberoser Heide’ referring to a generally flat (altitudes between 76 m and 90 m.a.s.l.) heath landscape. The region is characterized by the transitional Atlantic to continental climate. Long-term (1991–2020) average annual precipitation at the climate station Cottbus (ID 880; 51°46′33″ N, 14°19′3″ E, 69 m.a.s.l.) is 566 mm, with a clear peak during the summer season. The average annual temperature for the same period is 10.1 °C, whereby during the winter months, temperatures reach the point of freezing while the summer is warm.

Late Pleistocene influences shaped the region as it comprises parts of the southernmost German top moraine zones and features the complete glacial series, including large sand dune areas made up of gravel and sand. Sandy brown earth and podzol fine-grained soil with mainly silt and clay are prevailing [73]. The study site is completely unsettled. It is surrounded by forest with mainly Scots pine (Pinus sylvestris) stands and Silver birch (Betula pendula) trees to a lesser extent. The site was initially cleared by a forest fire in 1942 and subsequently used as central Europe’s largest military training area from 1945 to 1994 [74]. Later, the natural vegetation and soil were persistently disturbed and exposed to wind erosion, causing higher vegetation to settle very sparsely only [15]. Wind erosion caused the development of an initial mobile sand dune on the Pleistocene sand [17], which is the morphologically dominant feature in the otherwise flat area. Ever since the troops’ withdrawal in 1994, the study site has been under strict nature conservation, protecting it from human activities and ceded to natural succession. This led to the development of dry acidic grassland together with BSCs [17].

In the study area, the dry acidic grassland is dominated by Corynephorus canescens [12] (grey hair-grass) and heath vegetation dominated by Calluna vulgaris [16], and succession and spatial patterns of different BSC communities can be found. According to the various successional stages, these are: (i) initial green-algae BSCs (e.g., Zygogonium ericetorum and Klebsormidium flaccidum), (ii) moss-dominated BSCs (e.g., Polytrichum piliferum), and (iii) soil lichen (with mosses)-dominated biocrusts with various Cladonia species (e.g., Cladonia fimbriata and Cladonia coccifera) [15,75]. Different stages (e.g., photosynthetic activity, community respiration, and nitrogen fixation cycles) were characterized in various studies, e.g., [70,76]. Since the poikilohydric organisms, mosses, and lichens largely depend on the precipitation input, the changes in their hydration during dry–wet cycles influence their ecophysiological performance and structure, subsequently influencing the spectral reflectance [71,77].

2.2. Dataset Construction and Data Processing

Our analysis was built on a dataset constructed from in situ data, gridded time series of multispectral indices and bands derived from Sentinel-2 (S2) satellite images, and gridded topographic attributes (TAs) to parameterize the Random Forest (RF) classification approach for spatiotemporally estimating and mapping the BSC coverage in the study area (Figure 2). By incorporating daily precipitation data, the analysis of the seasonal and rainfall event-based BSC activity was driven.

2.2.1. In Situ Data Collection

General reference data for training and testing the RF classification model are the most meaningful when collected in situ. However, due to its status as a natural conservation area and ammunition contamination resulting from the time when the study area was used as a military training site, access for in situ data collection was not entirely possible. Therefore, reference areas of the four dominant land cover classes (i) biological soil crust (BSCs), (ii) trees, (iii) grey hair-grass, and (iv) bare soil were sampled based on expert knowledge from previous field trips and Google Earth Pro high-resolution satellite imagery from 2018 and 2019 (Figure 3). This was considered as feasible as the land cover in the study site is assumed to not show variations during our observation period (2015–2020). Only polygons featuring one clearly distinguishable single land cover class homogenously were used as reference areas (Figure 3). Due to the 10 m spatial resolution of the S2 data, only reference polygons featuring one single land cover class within 20 m × 20 m areas and an additional 10 m buffer to areas of the other classes were sampled to avoid mixed S2 pixels. Using the R Project for Statistical Computing (version 4.1.0) [78], these polygons were converted to points (Figure 2) where each represented one S2 pixel. Thus, mixed pixels were reduced, and splitting of the reference data into more training and testing pixels via the k-folding approach (c.f., Section 2.3) was enabled.

2.2.2. Meteorological Data

We retrieved hourly recorded precipitation data from the Climate Data Center portal of the German Weather Service DWD (https://cdc.dwd.de/portal/, access on 8 February 2021). The DWD station Lieberose (ID 2997; 51°59′10″ N, 14°18′18″ E, 47 m.a.s.l.) is located closest to our study site with a distance of 4.8 km in the northwest direction. Therefore, we acquired data from this station and calculated cumulative daily precipitation (mm) to be incorporated in our further analysis.

In addition, information on the hourly wind direction 10 m above ground (in degrees) recorded at the DWD station Cottbus (ID 880) at a distance of 17.3 km from our study area was retrieved. It was used to detect the average wind direction between 1983 and 2020 (190°, corresponding to south), which was incorporated in the calculation of the TA ‘wind effect’ (c.f., Section 2.2.4).

2.2.3. Sentinel-2 Data Preprocessing and Derivation of Multispectral Indices

We used Sentinel-2 (S2) data from the European Space Agency’s Copernicus mission for our analysis. With the twin constellation of the two polar-orbiting satellites S2A (launched in June 2015) and S2B (launched in March 2017) placed in the same sun-synchronous orbit, the temporal resolution of the system is five days. Both identical satellites carry an optical payload with visible (VIS), near-infrared (NIR), and shortwave infrared sensors (SWIR) [79]. The multispectral sensors operate with 13 spectral bands with a geometric resolution ranging between 10 m (four bands), 20 m (six bands), and 60 m (three bands). The swath width of S2 is 290 km. Data were retrieved from the ESA Copernicus Open Access Hub (https://scihub.copernicus.eu, accessed on 4 March 2021) and from the CREODIAS Repository platform (https://finder.creodias.eu, accessed on 4 March 2021).

S2A and S2B jointly captured a total of 358 images covering our study area, spanning the period from 4 June 2015 (S2A) to 29 December 2020. The S2 products were retrieved from descending sensing nodes as Level-1C (L1C) top-of-atmosphere (TOA) ortho-images, distributed in 100 × 100 km² tiles in WGS84 projection [80]. Our study area is located at the center of the tile 33 UVT. The L1C products were converted to Level-2A (L2A) products featuring bottom-of-atmosphere (BOA) reflectance corrected from atmospheric effects and terrain distortions by using the Sen2Cor preprocessor, implemented in the Sen2r R package [81]. Afterwards, the L2A tiles were cropped to the extent of the study area and reprojected to WGS84 UTM zone 33N (epsg code 32633). The further preprocessing aimed at cloud removal and cloud masking to minimize noise and enable a derivation of robust multispectral indices. We applied the Function of Mask (FMask, version 4.2) algorithm [82] to the L1C data, as it was proven to provide the most accurate results for S2 data [83]. It returned a classification for each scene, where areas covered by clouds, cloud shadow, and snow are depicted. These were subsequently masked out to only proceed with genuine land surface information under clear-sky conditions.

From a total of 358 S2 acquisitions, 220 (61%) images exhibited full cloud coverage, including shadows. From the remaining 138 scenes, 49 (13%) of the images were acquired under full clear-sky and snow-free conditions. By applying a threshold of 20% valid surface information within the study area, 107 S2A/B L2A products, spanning the years from 2015 to 2020 (Figure 4), were identified as appropriate for our analysis. Eventually, the number of “analysis-ready” scenes for our analysis were as follows: 5 scenes in 2015, 10 scenes in 2016, 9 scenes in 2017, 31 scenes in 2018, 29 scenes in 2019, and 23 scenes in 2020 (Figure 4).

While preprocessing the S2 scenes, we detected an offset between the data recorded by S2A and S2B. While S2A itself is very precise, S2B data deviates to S2A as well as among itself. The deviation between the two sensors could be quantified to up to 10 m, equivalent to one pixel, mainly in the north and east directions. This, of course, introduces inaccuracies when making pixel-based comparisons. However, as the majority of the scenes are congruent and we constantly incorporated a wide range of pixels, we chose to proceed without co-registering the images. As this is only doable for satellite scenes with low cloud coverage, this would not have been applicable to all images. Of course, we have to keep this in mind for the upcoming analysis.

Only the four bands available at the highest resolution of 10 m (i.e., blue, green, red, and NIR) were used in this study, as they sufficed for generating all relevant indices. We computed the following four different multispectral indices that were reported to be meaningful for RS-based analysis of BSCs for each individual scene using Equations (1)–(4): NDVI [49], CI [26], mCI [28], and BSCI [27]:

NDVI = (R_NIR − R_red)/(R_NIR + R_red),

(1)

CI = 1 − (R_red − R_blue)/(R_red + R_blue),

(2)

mCI = 1 − (R_green − R_red)/(R_green + R_red)

(3)

BSCI = (1 − L × |R_red − R_green|)/R^mean_greenredNIR

(4)

where R_blue, R_green, R_red, and R_NIR are the surface reflectance ranging from 0 to 1 at the blue, green, red, and NIR bands, respectively, corresponding to S2 band numbers 2, 3, 4, and 8. L is an adjustment parameter to amplify the absolute difference between R_green and R_red. It was set to the value of 2, and R^mean_greenredNIR is the mean reflectance of the green, red, and NIR bands [27].

2.2.4. Topographic Attributes

Topographic attributes (TAs) are commonly applied to quantitatively describe landform characteristics and structures [84] and as proxy-drivers in environmental modeling, e.g., [85], and are readily available as computed derivates from Digital Elevation Models (DEMs) [86]. In Brandenburg, the Federal State’s Department of State Survey and Geospatial Base Information provides DEMs in a high spatial resolution of 1 m (post-processed airborne LiDAR data acquired on 19 March 2010) free of charge through its web application (https://geobroker.geobasis-bb.de, access on 13 June 2020). A mosaic of six DEM tiles (each with 1 × 1 km²) was converted to GeoTiff and cropped to a bounding box capturing the study area. From this, a total of 13 TAs (

Table 2. Topographic attributes calculated in SAGA GIS and their relevance as predictors used in the Random Forest classification model for mapping the biological soil crusts in the study area.

Attribute (unit)	Relevance
Aspect (degree)	the landscape’s spatial heterogeneity affecting the surface energy balance and through this, the soil water retention capacity and water availability [88]
Elevation (m.a.s.l.)	the physical landscape properties and their spatial patterns, and key attributes for the further derivation of terrain-shape indices
Flow accumulation (-)	the effects of the depth and velocity of flow (here by integrating the multiple flow direction based on the maximum downslope gradient; top-down processing) [89]
Insolation (W m−2)	the annual sum of direct and diffuse potential incoming solar radiation calculated according to times of sunrise and sunset [90]
Topographic openness (positive and negative; rad m−1)	the dominance (positive) or enclosure (negative) of a landscape location [91]
Plan curvature (rad m−1)	the contour line formed by intersecting a horizontal plane with the surface, thus a proxy on the convergence or divergence of water during downslope flow [88]
Profile curvature (rad m−1)	the surface in the direction of the steepest slope, thus a proxy on the acceleration/deceleration of surface water [88]
Slope angle (degree)	the landscape’s spatial heterogeneity and catchment-related hydrological processes (e.g., flow direction, water runoff velocity and accumulation) [88]
Topographic Position Index (TPI; -)	the topographic slope positions and for landform classifications [92]
Terrain Ruggedness Index (TRI; -)	the surface heterogeneity of the landscape; averaged from the absolute differences in elevation between a focal cell and its eight neighboring DEM cells [93]
Topographic Wetness Index (TWI; -)	the topography’s spatial scale effect on hydrological processes and proxy on the terrain-related soil moisture potential [94]
Wind Exposition Index (WEI1; -)	the wind shadowed pixels/areas (values < 1) and pixels/area that are exposed to wind (values > 1) [95]
Wind Effect Index (WEI2; -)	The direction in which the wind is going (windward/leeward index) [95]

) were calculated representing continuous primary local TAs (i.e., elevation, aspect, slope angle, plan, and profile curvature) and secondary compound TAs (i.e., flow accumulation, topographic position, ruggedness, wetness index, and wind exposition and effect). Their DEM-based computation was performed using SAGA GIS (version 2.3.2) [87]. To meet the same pixel size of the S2 input data, we resampled the gridded TAs to a joint geometric resolution of 10 m to be incorporated into our modeling framework supporting the classification (Figure 2).

2.3. Biocrust Mapping and Response to Rainfall

2.3.1. Estimating Biological Soil Crust Coverage

We applied an RF classification analysis [67] using the R-package randomForest version 4.6-14 [96] to estimate the spatial distribution of the BSC cover by classifying the study area into the four dominating land cover classes (c.f., Section 2.2.1). RF is an ensemble ML method for classification and regression, combining numerous tree predictors voting for the most popular class, growing by selecting random features in the training dataset [67]. It also works well using high-dimensional input variables [97].

The model requires two input parameters: reference areas and predictors (Figure 2). Our sampled reference areas (c.f., Section 2.2.1) included 1,085 S2 pixels partitioned into 322 (29.7%) for the class ‘biological soil crusts’ (BSCs), 73 (6.7%) for the class ‘trees’, 564 (52%) for the class ‘grey hair-grass’, and 126 (11.6%) for the class ‘bare soil’ while the entire study area is captured by 53,139 S2 pixels. The predictors are the parameters the classification was performed on. In our study, these were the multispectral data (i.e., the calculated indices and the individual bands in 10 m spatial resolution; c.f., Section 2.2.3) and the TAs (Table 2). As we aimed at one classification valid for the full timeframe, all multispectral parameters were integrated over the whole observation period (2015–2020) using all 107 analysis-ready acquisitions (c.f., Section 2.2.3) via the pixel-based mean. This was considered feasible, as the land cover in the study site is not subject to short-term variations. The RF algorithm resulted in one probability map for each class (c.f., Section 3.2), containing each pixel’s probability of matching the respective class. In the discrete classification map (c.f., Section 3.3.1), pixels were strictly assigned to the class with the highest probability.

Model training and testing and the accuracy assessment were performed using k-fold cross-validation [98], a commonly used method where the reference data is randomly split into k (in this study ten) disjunct subsets. The model was then run ten times, each time using nine of the subsets for training and one of the subsets for testing (Figure 2). The results of each run were assessed and averaged across all ten replicates. This approach was superior to simply splitting the data into strictly one portion for training and one for testing, as this could lead to overfitting and biased accuracy estimations.

Our RF-based approach also provided measurements for evaluating the importance of the predictors to the performance of the model: (i) the mean decrease in accuracy (MDA), measuring how much accuracy is lost in percent when the respective variable was permuted from the model and (ii) the mean decrease in Gini coefficient (MDG), measuring how much each variable contributes to the homogeneity of the nodes in the resulting RF model [99].

To obtain information on the classification results corresponding to reality, we conducted an accuracy assessment [100] by applying the trained model on the test data, i.e., the portions derived from the reference dataset by the k-folding approach. Therefore, a confusion matrix was used to compare the classified pixels with the test data for each fold and then averaged. Four measures of accuracy were calculated for each class based on this matrix [100]: (i) overall accuracy (OA), revealing the overall model performance, (ii) user’s accuracy (UA), providing information on how well the classification results match the classes in reality, and (iii) producer’s accuracy (PA), indicating the proportion of correctly classified pixels. The kappa coefficient of agreement (K^) returns the extent to which the result differs from a random classification.

2.3.2. Analyzing the Seasonal and Rainfall Event-Based Activity of Biological Soil Crusts

The areas classified as BSCs were subsequently used to retrieve information on dry and wet cycles. We used the NDVI for assessing the change in metabolic activity of BSCs, i.e., under wet conditions after precipitation events and under dry conditions, when the activities are reported to be suspended (c.f., Section 1). We compared the S2-based NDVI data from the available 107 acquisitions (c.f., Section 2.2.3) to daily precipitation data measured at the weather station Lieberose (c.f., Section 2.2.2). Hereby, we evaluated whether the postulated effect can be depicted in our study area using the characteristics of the S2 sensor in terms of spatiotemporal resolution. A one-way ANOVA was tested for explanatory rainfall response examples to investigate if the NDVI values for the class ‘biological soil crust’ (BSC) differ significantly under dry and wet conditions. Additionally, the response behavior was evaluated using change detection between NDVI grids representing dry and wet conditions, respectively, to visualize spatial patterns of change in BSC’s activity.

Moreover, we investigated seasonal trends of BSCs. To increase the temporal density of our time series and obtain meaningful information on the average seasonal trend, we merged all NDVI grids from five years into one. This was especially important for the earlier years, where only one sensor (i.e., S2A) was operating, causing a low temporal resolution (Figure 4). A polynomial trend line was calculated to illustrate the temporal course of activity through one average year. The range of NDVI values was compared for each season, tested for significant differences using one-way ANOVA. The seasons are hereby defined as winter (December–February), spring (March–May), summer (June–August), and autumn (September–November).

3. Results

3.1. Performance of the Random Forest Classification Model

The accuracy assessment returned an OA of 98.9% and a Κ^ of 0.9822 when validated against the 10-fold based test dataset. The UA is 96% for the class bare soil, 98.6% for trees, 99.1% for grey hair-grass, and 99.6% for BSCs. The values of PA are 96% for bare soil, 98.9% for grey hair-grass, 99.6% for BSCs, and 100% for trees. The higher the classes’ accuracies, the more the algorithm is capable of distinguishing them, resulting in less misclassification. The classification can therefore be categorized as highly accurate.

This also provides legitimacy concerning the predictor importance information (Figure 5). Hereby, the predictors are ranked according to their importance in estimating and differentiating the land cover classes as determined by the RF classification model. Higher values refer to a higher importance of the corresponding parameter within the classification.

In general, the S2-based multispectral predictors, i.e., the spectral bands and indices, were detected as more important than the TAs employed. Hereby, the NDVI is by far the most crucial variable (Figure 5). Concerning the spectral bands, those sensing in the visible spectra (i.e., red, green, and blue) were identified as more important than the NIR band (Figure 5). With regard to the multispectral indices, the mCI ranks second, followed by the CI and the BSCI. Among the TAs, elevation was identified as the most important parameter (Figure 5), revealing approximately the same parameter importance as the green spectral band and the NIR (Figure 5). Descending importance was revealed for the TAs representing the exposition to and the effect of wind (WEI₁ and WEI₂), the insolation, and the topographic openness in the study area. Flow accumulation, TPI, and curvature contributed the least to the classification (Figure 5).

3.2. Estimated Biological Soil Crust Coverage

As revealed in the RF-based classification, our study area (531.4 ha in total) classifies into 10.07 ha (1.9%) of bare soil, 287.3 ha (54.1%) of BSCs, 219.25 ha (41.3%) of grey hair-grass, and 14.77 ha (2.8%) of trees (Figure 6 and Figure 8d). Hence, BSCs were identified as predominating land cover within the ‘Lieberoser Heide’. Tree cover is mainly present in the northwest of the study area and south of the sand dune (Figure 6). The sand dune was classified as the most extensive bare soil feature. With a small transitional zone to the bare soil, grey hair-grass was classified, creating a buffer surrounding the dune. Grey hair-grass also dominates the central and eastern study area (Figure 6c), where biocrust-dominated patches gradually blend in (Figure 6c), creating large transitions. BSCs were classified as expansively and continuously present in the western and northern study area (Figure 6b). In the very northwest, they mix with the class ‘trees’.

In the final classification map (Figure 8d), pixels were strictly assigned towards their dominant land cover class. This implies the presence of minor patches of other land cover types, and the assigned one has the highest share and therefore dominates the spectral signature. In the case of the study area, this is especially true regarding the intersection of the classes ‘biological soil crust’ (BSC) and ‘grey hair-grass’, and the fluctuating density of grey hair-grass on the sandy bare soil areas. Therefore, information on the transitional zones of the classes can be derived from the probability maps showing the probability for each pixel to be assigned to an individual class (Figure 6a–d).

3.3. Dry and Wet Cycles of Biological Soil Crusts

3.3.1. Activity of Biological Soil Crusts in Response to Rainfall Events

In order to obtain information on the direct response of the areas classified as BSCs to rainfall, at least two cloud-free scenes close to rainfall events were required, i.e., to capture the BSC’s metabolic activity represented by the NDVI before and after rainfall. However, only 107 adequate S2 scenes over the period of 2000 days led to an average temporal resolution of around 18.7 days (Figure 4), making the time series relatively coarse, thus likely constrained to accurately map short-term changes. Due to these inconsistencies, a meaningful analysis over the whole observation period was not feasible. Therefore, explanatory periods where the temporal resolution was assumed sufficient to depict rainfall-induced changes were selected. We focused on the two most significant ones, located in October/November 2018 and March/April 2020 (Figure 7).

The first example period shows a longer dry period between DOY 277 and DOY 294, whereby the NDVI of the BSC areas exhibit median values of 0.3 to 0.35 (Figure 7a). Within the nine consecutive days, 21.3 mm of rain was recorded at the climate station in Lieberose (c.f., Section 2.2.2), after which the median NDVI increased significantly up to 0.52 (Figure 7a). During the following eleven dry days (DOY 304 to 314 with a cumulative rainfall of 0 mm), the NDVI decreased to its initial level at DOY 294 (Figure 7a). A similar behavior of NDVI (i.e., increase in NDVI after a rainfall event) can be observed between DOYs 314 and 319 (Figure 7a). Therefore, two levels of NDVI values can be observed: NDVI values around 0.55 for wet and NDVI values around 0.33 for dry conditions (

Table 3. Descriptive statistics of the Sentinel-2-based NDVI observations for the two example periods from day-of-year (DOY) 273 to 321 in October/November 2018 and from DOY 68 to 120 in March/April 2020 for each land cover class separated for dry and wet conditions. Hereby, the NDVI values for DOY 84 in 2020 were omitted, as they represent intermediate conditions.

Statistics	Total	Biological Soil Crusts		Bare Soil		Grey Hair-Grass		Trees
		Dry	Wet	Dry	Wet	Dry	Wet	Dry	Wet
Nr. of observations (i.e., pixels)	535,875	208,967	86,190	6860	3021	149,595	65,775	11,036	4431
Minimum	0.073	0.138	0.248	0.073	0.078	0.090	0.099	0.278	0.402
1st quartile	0.273	0.301	0.520	0.108	0.124	0.231	0.392	0.464	0.636
Median	0.329	0.329	0.548	0.123	0.159	0.255	0.451	0.514	0.681
3rd quartile	0.450	0.364	0.575	0.141	0.199	0.275	0.494	0.575	0.736
Maximum	0.894	0.645	0.759	0.439	0.517	0.544	0.680	0.805	0.894
Average	0.361	0.335	0.546	0.126	0.168	0.251	0.438	0.522	0.686
Standard deviation (SD)	0.121	0.047	0.045	0.026	0.054	0.037	0.076	0.077	0.070
Coefficient of variation (%)	33.2	14.1	8.3	20.8	32.2	14.6	17.3	14.8	10.2

).

The second example period (Figure 7b) confirms this observation, as the two different NDVI levels can be seen at DOY 74 and from 99 to 120, respectively. Additionally, DOY 84 shows an intermediate value of decreasing NDVI, indicating the process of gradual drying. Interestingly, the precipitation event of 3.4 mm at DOY 104 did not affect the NDVI values recorded five days later, which are at the same level as before the rainfall (Figure 7b).

Table 3 depicts the descriptive statistics of the NDVI values for dry and wet conditions over all land cover classes within the example periods (Figure 7). It can be observed that bare soil values are the least subject to change from dry to wet conditions, with the median changing by approximately 0.04 (Table 3). The median changes from dry to wet for grey hair-grass amount to approximately 0.17, and to approximately 0.18 for trees and 0.22 for BSCs, respectively (Table 3). Thus, BSCs are observed to exhibit the highest changes in median NDVI between dry and wet conditions.

The rainfall response was also analyzed on a spatial basis. Figure 8 shows the associated NDVI grids for the observations at DOY 74 and 99 in 2020 for the entire study area, including all land cover classes. Here, DOY 74 refers to wet conditions, whereas the observation at DOY 99 represents dry conditions, i.e., after 10 days without rainfall (Figure 7). Additionally, the change between the two NDVI observations is depicted indicating that areas in the western and northern (Figure 8c) areas are subject to the highest changes in NDVI.

3.3.2. Seasonal Trends in the Activity of Biological Soil Crusts

The number of observations (i.e., S2 scenes) during spring and summer is much higher than during autumn and winter (c.f., Section 2.2.3). However, when merging all available 107 NDVI scenes capturing the entire observation period into one single aggregated year, a significant trend and seasonal differences showing a distinct phenological (i.e., referring to the metabolic activity) pattern can be observed (Figure 9). The NDVI values are significantly higher during the winter season (median 0.49, SD 0.08) when compared to the spring (median 0.34, SD 0.06) and summer (median 0.3, SD 0.03) seasons (Figure 9b). In autumn, the NDVI increases again significantly (median 0.36, SD 0.08) (Figure 9b). When comparing the autumn to the spring season, no significant differences in NDVI can be observed (Figure 9).

4. Discussion

4.1. Land Cover Classification and Biological Soil Crust Coverage Estimation

The very high overall RF-based classification accuracy (98.9%) indicates that our applied classification approach is meaningful. It corresponds well to the reality, and hence is considered successful in delineating BSCs in our study area. Particularly, the class ‘trees’ exhibited very high producer’s accuracy (c.f., Section 3.1) as they represent distinct features when compared to BSCs, bare soil, and grey hair-grasses, thus are clearly distinguishable from all other. There is only a tiny portion of tree pixels wrongly assigned to the BSC class, originating from the northwestern study area. Here, they are presently isolated within and on top of BSCs. Bare soil exhibits the lowest accuracy—which is still very high—as it tends to be confused with grey hair-grass. This originates from sampling the latter primarily in areas where bare soil is the background substrate, like in the immediate surroundings of the sandy top of the dune. Grey hair-grass is never present alone, as it grows in clusters on top of mainly sand and BSCs. When it mixes with BSCs, it is not distinguishable as the small isolated tufts do not highly contribute to the spectral signature on these scales. This is also why grey hair-grass shows the most spectral similarities with the bare soil values for all of the indices. It is surprising that BSCs and grey hair-grass do not get confused more often, as they resemble themselves concerning CI and mCI indices and show large transitional zones in the area of interest. Additionally, regarding the NDVI, they reveal similar values. The BSCI, however, is able to differentiate them quite well.

When classifying the different land cover types, the multispectral predictors reveal higher importance than the TAs (Figure 5). Hereby, the NDVI is by far the most crucial predictor. This is not very surprising, as it is a known meaningful index concerning plant physiognomy and activity, and coincides with previous observations, where the index band combination NIR and red was found to be the most crucial predictor in estimating BSC cover [28]. The mCI was identified as the second most important index, as it is able to differentiate BSCs from other types of higher plants [28]. Additionally, the CI, targeting cyanobacteria [26], and the BSCI, developed especially for lichen-dominated biocrusts [27], performed well (Figure 5). Concerning the individual spectral bands, red, blue, and green show importance comparable to the indices, with only NIR being a little less important (Figure 5). This confirms that adding the raw bands and not only relying on indices is able to increase the classification accuracy.

For our study, we cannot conclude that the distribution of the different land cover classes is highly influenced by the topography, as the TAs were identified as not distinctly important in the classification. However, from all TAs, elevation reveals the highest importance (Figure 5). This originates from the dune on the one hand, as most of the bare soil is located on the ridge. On the other hand, the fact that grey hair-grass is mainly located in the lower central-eastern area at elevations from 75 to 80 m.a.s.l., and biocrust coverage is highest in the northwestern plateau, located at higher elevations from 85 to 90 m.a.s.l., certainly also influences the importance. There might also be some bias in the training data, as the distribution of the biocrust samples is concentrated on this plateau. However, this was inevitable, as this area is the only one where the BSC cover was known to be high from previous investigations. Similar to the elevation, the relatively high-ranked wind attributes WEI₁ and WEI₂ (Table 2) are influenced by these effects, as the highest wind effects are also detected for the plateau and the dune, which the bare soil class is mainly located on top of. This confirms that the wind attributes represent reality quite well, as, without the wind effect, no bare sand areas would exist at the dune.

It is arguable whether including the TAs into the classification is beneficial, as the parameter importance implies otherwise (Figure 5). Therefore, we ran the RF model a second time using multispectral parameters only, and the resulting classification did not show high differences in terms of accuracy. Thus, we conclude that including TAs of 10 m spatial resolution into the RF-based classification was not highly beneficial in our study site, as the prevalent very small-scale topographic changes are severely smoothed when resampling from 1 to 10 m and therefore cannot be recognized anymore. The probability for each class was tested for correlation to the topography using the Pearson correlation coefficient. While other classes showed high correlations, BSCs only show small correlations to a few TAs, with 0.27, 0.16, and −0.12 being the highest numbers for Slope, WEI1, and Elevation, respectively. This confirms the statement that the biocrust coverage is not highly dependent on topography on a 10 m scale.

The land cover classes of the derived classification have to be interpreted as areas where the respective class is dominant in the spectral signature, as they are not strictly separable in reality and show transitional zones. This also implies that pixels assigned to one class may feature shares of other classes influencing the signal, which has to be considered when interpreting the rainfall response of classes. The approach of differentiating the four selected land cover classes is ideal for obtaining an uncomplicated and fast delineation of the BSC-covered areas. Using the temporal average for the multispectral predictors performed well, as confirmed by the high classification accuracy. Comparing the results to reality showed a plausible coincidence to expectations and could be evaluated on-site in a short field visit in the outermost northwest study area in Sep. 2020. We would like to reiterate that access to the study area is generally prohibited due to its natural conservation status (c.f., Section 2.2.1). We could observe that even though grey hair-grass is present in areas of high BSC probability, their percentage share can be neglected with regard to its spectral signature at a spatial resolution of 10 m, as the fascicles are considered as not distinctly contributing to the spectral signature.

4.2. Dry and Wet Cycles of Biological Soil Crusts

While the precipitation data is continuously measured at the climate station, the S2 NDVI time series is limited by two factors: the revisiting time of Sentinel 2 and the quality of the available scenes. Revisiting time is five days since the launch of S2B in March 2017. Before that, images are available only every ten days, severely limiting temporal resolution for the first two operational years. The highest drawback originates from clouds and snow, as affected areas cannot be used for the analysis. With the high probability of clouds, the most considerable portion (251/358, 70%) of satellite acquisitions is not appropriate for analysis (c.f. Section 2.2.3). This is one of the main factors that distinguishes areas in temperate climates from the semi-arid and arid drylands, as these show a much lower probability for cloud cover, making optical imagery available on a much more consistent basis. With the remaining 107 NDVI observations, the theoretical temporal resolution drops to 18.7 days, making it hard to connect them to precipitation events, as BSCs react quickly to them by recommencing their metabolic activities within a few hours. Often the observations cannot be easily assigned to dry or wet periods, as precipitation and drying events of different magnitudes alternate between them. Additionally, during summer the high temperatures and solar irradiation lead to higher evaporation rates, making the BSCs dry out much faster than during the other seasons and leading to no S2 observation able to detect the activity changes. Therefore, spring and autumn are the best periods to identify the rainfall response for our study area.

4.2.1. Activity of the Biological Soil Crusts in Response to Rainfall Events

The two selected explanatory periods of the time series show that the NDVI differs significantly (p < 0.001) for biocrusts under wet and dry conditions, and the change can be affiliated to precipitation events. When wet, BSCs tend to match the spectral characteristics and NDVI of higher green vegetation, such as fully grown crop fields [45], as the median reaches values of 0.55 (Table 3). After dry periods, the median NDVI drops down to 0.33, getting closer to bare soil values, which are below 0.16. These values and the magnitude of change (0.22) match observations from areas with high biocrust coverage within drylands [39].

We can also confirm the observation of studies performed in the Negev [2,37] that activity of BSCs persists for up to five days after precipitation events when the conditions are right. For example, the image for DOY 84 in 2020 was acquired five days after a small precipitation event of only 0.4 mm. The NDVI is still significantly higher than on DOY 99 after ten dry days, but also significantly lower than on DOY 74 after multiple days of larger rainfall events. This means that even such small precipitation events can influence the activity of BSCs over a period of five days, at least in spring, when the evaporation rates are low due to moderate solar irradiation and temperatures. It is arguable that the two minimal rainfall events of 0.4 mm each might just have prolonged the photosynthetically active time from DOY 69 to 73, where the biocrusts go into a “hold position” after large precipitation events and can be recommenced by small rainfall events [2,37]. In November 2018, in phenologically similar conditions, ten dry days were enough for the biocrusts to return to the same level of activity as before the rainfall event. The two observations of wet conditions were recorded one and two days after the precipitation event in this example, which means that during this phase, the activity has not declined.

The initiation phase of biocrusts after rainfall events, where the activity increases, could not be evaluated due to the lack of NDVI scenes immediately after precipitation events. We, however, could register that one day is enough for an activity resumption (c.f., Section 3.3.1). Here, the inconsistencies of the S2 NDVI time series are once again the limiting factor, as the most interesting parts are on the days directly before and after precipitation events. However, these days are more likely to be cloud covered, reducing the chance of obtaining an appropriate observation. In this study, there are hardly any images found to be acquired on the day of a large precipitation event or the next day.

The effect of rainfall events on the NDVI is not consistently clearly visible even though an observation might be available. The observation from DOY 2020–109, for example, is located just five days after a larger rainfall event amounting to 3.4 mm. While we already showed that five days are enough to detect activity even after small rainfall events, this is not the case for this observation, as the NDVI stays at the same level. It is unclear whether here the activity has already declined again due to a change in temperature and higher evaporation, or whether the rainfall was only recorded by the station but did not occur in the study area. Rainfall can exhibit high local variability on a small spatial scale, and we would like to stress that the climate station Lieberose (c.f., Section 2.1) is located at a distance of 4.8 km northwest and at a distance of 8.8 km southeast of our study area. This is assumed to not be close enough to ensure that each event occurring within the study site was detected by the station. The best example confirming this is the 18 mm of rainfall measured on 9 September 2019 (DOY 252), which only led to a very minor change in NDVI values of 0.05 (median) detected for the class ‘biological soil crusts’ two days later.

The median NDVI values of bare soil are constantly below 0.16 (Table 3), even under wet conditions, confirming that our classification performed well in delineating this class. Trees show the highest NDVI values as expected, but also a very high value range. This is most certainly because the trees within the study area are isolated and therefore the pixels may be influenced by signals from the surrounding land cover. The values of grey hair-grass are always lower than the ones of BSCs under both dry and wet conditions, with this difference being more pronounced after rainfall. While bare soil medians stay stable under dry and wet conditions and only their value range gets a little higher, grey hair-grass and trees show the same dry/wet pattern as BSCs. This originates from the discrete classification featuring dominant land cover (c.f., Section 4.2). As BSCs can exist below both grey hair-grass and trees, they of course influence the NDVI values. The fact that the highest difference of NDVI between dry and wet conditions was detected for areas classified as BSCs confirms the validity of our classification (Table 3), as the rainfall response patterns of NDVI are known to differ for various degrees of sub-pixel biocrust cover and are more significant with higher coverage [39]. When we transfer this hypothesis to the NDVI change map (Figure 8c), the pixels with the highest magnitude of change are the ones with the highest sub-pixel biocrust coverage, which are the ones in the very northern and central-western areas. In the very northwest, where the biocrust cover is known to be very pure and high, however, the change magnitude is not that high, which is the result of the trees scattered within this area returning a mixed signal.

4.2.2. Seasonal Trends in the Activity of Biological Soil Crusts

Due to the higher cloud and snow cover probability during the autumn and winter seasons in temperate climates, data gaps (resulting from fewer satellite observations available for analysis) are more likely than during the spring and summer seasons (Figure 4). However, aggregating all multitemporal data into one single year made it possible to generate a significant average seasonal trend for BSCs (Figure 9). It reveals patterns of BSCs with significantly higher average NDVI values during the winter season when compared to the summer season. This seems unusual for any vegetation within the temperate climates, as vascular plants usually show higher activity during the warm summer months before dying off or hibernating during the cold winter. However, BSCs are considered poikilohydric, i.e., being only active when sufficient water is available [2]. As in temperate regions, incoming solar radiation, and therefore air and soil temperatures, are higher during the summer [101]; higher evaporation rates are supposed to lead to less water available for, and therefore less metabolic activity of, the BSCs (Figure 9). This proposes a higher probability to detect active BSCs during the winter; however, this cannot be proven, as no agreeing studies could be identified.

For the years 2018 and 2019, a higher amount of adequate (i.e., no to only a few clouds) RS-based observations were available due to the occurrence of extraordinary hot and dry summers, being ranked 2nd and 3rd warmest since 1881 by the DWD for all of Germany [102]. The effects in 2018 were even more pronounced in the northern and eastern parts of Germany, including our study area [103]. When analyzing the rainfall data from the Lieberose weather station (c.f., Section 2.2.2), the long-term (1991–2020) average annual precipitation was 557 mm, while in 2018 and 2019, only 398 and 467 mm of rainfall was recorded. This certainly had an influence on the seasonality of BSC activity as less rainfall and higher solar irradiation, and therefore higher evaporation, lead to less photosynthetic activity in the BSCs, subsequently lowering the NDVI values. Thus, the averaged trendline and significant differences between the summer and the winter seasons (c.f., Section 3.2) might be much more pronounced. With global warming, the temperature around the study area is expected to rise in the future, leading to more precipitation in the form of rainfall than snow and fewer days below freezing during winter [104], implying an increase in the activity of BSCs during the cold term. On the other hand, BSCs are expected to show fewer and shorter periods of metabolic activity during the summer as they become exposed to higher evaporative stress due to lower precipitation rates and higher temperatures [3].

4.3. Future Research Potential to Enhance the Remote Sensing-Based Monitoring of Biological Soil Crusts

The main limiting factor of this study lies in the availability of the S2 observations. To realize a more regular future monitoring, the NDVI time series needs to be more consistent. Using the present NDVI data in harmonizing or gap-filling algorithms utilizing seasonal trends is not an option, as the individual response to rainfall events would be lost. Therefore, expanding the current database with additional data is more reasonable. For example, simulating S2 NDVI on 10 m resolution using NDVI values derived from comparable optical sensors such as Landsat 8 [69] (30 m resolution) based on image fusion techniques [105] can be a relatively straightforward method to increase the database while maintaining the high resolution. However, cloud coverage would still be a problem as the data gaps during cloudy periods will still be present. To eliminate this limitation, SAR sensors are beneficial due to their ability to penetrate clouds. SAR backscatter information is used to estimate the near-surface soil moisture and surface roughness [106], both parameters changing for BSCs under dry and wet conditions [107], and is thus capable of detecting the response of BSC-covered areas to rainfall [57]. SAR-based information could also support the delineation of classes, and hence improve the classification of BSCs.

A consistent monitoring would allow quantifying and characterizing the rainfall response based on the frequency and magnitude of rainfall events, as well as the influence of different weather conditions (evaporation rates) on BSC activity. Moreover, more precise, and spatially explicit data, such as on rainfall, would be needed to reduce uncertainties. Additionally, meteorological data (e.g., precipitation, temperature, insolation) and surface soil moisture information recorded by an observational network within or in the immediate surrounding of the study area is expected to increase the significance of the estimation of the response of BSC’s activity to rainfall.

Another limitation lies in the characteristics of the study area. Variations in land cover occur on small scales that are simply not detectable using a 10 m pixel size. Therefore, optimizing the workflow towards the implementation of data of higher spatial resolution (i.e., <10 m) would be an interesting option. With this, the topography is assumed to become more important in estimating the coverage of BSCs, as all the small topographic features would be recognized (c.f., Section 4.1).

With possible spatial resolutions in the range of centimeters, UAV-based multispectral and thermal observations of the study would enable for distinguishing individual grey hair-grass tussocks from BSCs, or even different types of biological soil crusts, e.g., lichen- and moss-dominated [58,59]. This information could be used for generating more precise reference data, as the protection status and the ammunition contamination do not allow accessing the study area. Additionally, information on the share of biocrust cover below each satellite pixel can be generated using UAV-based information [28]. This would be useful to validate the hypothesis that higher biocrust cover leads to higher rainfall response (c.f., Section 4.3), hence allowing a more detailed analysis. Additionally, multitemporal UAV-based observations can be useful for identifying and validating the rainfall response [28]. This approach can also be used to model the photosynthetic characteristics and the ecological functions of biocrusts and their spatial variance on a larger scale, as recent studies located around our study area typically derived these measures from very small-scale (<1 m²) NDVI plots [70,72].

5. Conclusions

The ultimate aim of our study was to evaluate the potential of Sentinel-2-based time series toward biocrust research within the temperate regions in an area of limited access. Our results prove that even with scarce data availability, but rather with the combination of expert knowledge for gathering training and testing data and high-resolution optical satellite imagery, we were able to successfully detect BSCs and map their activity in response to rainfall by transferring methods and meaningful BSC indices previously developed for arid and semi-arid drylands to a temperate dry acid grassland in Central Europe. Our RF-based classification algorithm provides highly accurate results when using multispectral data, biocrust-specific indices, and topographic information. Particularly, the spectral indices (e.g., NDVI) employed were identified as crucial for estimating and mapping the BSCs. Satellite-based NDVI measurements were also revealed as supportive for quantifying the response of BSCs to rainfall events and their phenological cycle. Using Sentinel-2 acquisitions with a spatial resolution of 10 m for our small-scale study area (~530 ha) was found to be generally suitable. Our multitemporal, spatial high-resolution investigations mark the first time that rainfall response of BSCs was examined over an extended observation period of around five years. However, the high cloud cover probability leads to a considerable loss in data availability for observing the quick response of BSCs to rainfall and differences between dry and wet cycles, constraining a continuous monitoring. Here, our applied workflow provides potential to incorporate further meaningful data. Moreover, impulses for its transfer to similar regions in Germany and Central Europe for estimating and mapping BSCs coverage can be derived. This is of particular relevance for linking our spatiotemporal analysis to those focusing on characteristics and ecological functions of BSCs and their variance on a spatial scale. Moreover, our study emphasizes the potential for future studies that aim at refining the modelling and enhancing the understanding of spatiotemporal dynamic of BSCs, not only in the global drylands but also beyond. The results might be of particular relevance for stakeholders in the natural conservation domain and show the high value of unique areas such as the ‘Lieberoser Heide’ to ecological research.