UAV RGB, thermal infrared and multispectral imagery used to investigate the control of terrain on the spatial distribution of dryland biocrust

Abstract Biocrusts (topsoil communities formed by mosses, lichens, bacteria, fungi, algae, and cyanobacteria) are a key biotic component of dryland ecosystems. Whilst climate patterns control the distribution of biocrusts in drylands worldwide, terrain and soil attributes can influence biocrust distribution at landscape scale. Multi‐source unmanned aerial vehicle (UAV) imagery was used to map and study biocrust ecology in a typical dryland ecosystem in central Spain. Red, green and blue (RGB) imagery was processed using structure‐from‐motion techniques to map terrain attributes related to microclimate and terrain stability. Multispectral imagery was used to produce accurate maps (accuracy > 80%) of dryland ecosystem components (vegetation, bare soil and biocrust composition). Finally, thermal infrared (TIR) and multispectral imagery was used to calculate the apparent thermal inertia (ATI) of soil and to evaluate how ATI was related to soil moisture (r 2 = 0.83). The relationship between soil properties and UAV‐derived variables was first evaluated at the field plot level. Then, the maps obtained were used to explore the relationship between biocrusts and terrain attributes at ecosystem level through a redundancy analysis. The most significant variables that explain biocrust distribution are: ATI (34.4% of variance, F = 130.75; p < 0.001), Elevation (25.8%, F = 97.6; p < 0.001), and potential solar incoming radiation (PSIR) (52.9%, F = 200.1; p < 0.001). Differences were found between areas dominated by lichens and mosses. Lichen‐dominated biocrusts were associated with areas with high slopes and low values of ATI, with soil characterized by a higher amount of soluble salts, and lower amount of organic carbon, total phosphorus (Ptot) and total nitrogen (Ntot). Biocrust‐forming mosses dominated lower and moister areas, characterized by gentler slopes and higher values of ATI with soils with higher contents of organic carbon, Ptot and Ntot. This study shows the potential to use UAVs to improve our understanding of drylands and to evaluate the control that the terrain has on biocrust distribution.

For example, biocrusts can increase water availability for plants by augmenting water retention in topsoil (Eldridge et al., 2020) and reducing soil evaporation (e.g., Adessi et al., 2018;Chamizo et al., 2016). They can also modify the erosion rate (Gao et al., 2020), which affects sediment accumulation, C and nutrient content in soils (e.g., Chamizo et al., 2017) and surface roughness (e.g., Rodríguez-Caballero et al., 2015;Wang et al., 2017). For this reason, mapping the spatial patterns of different biocrust-dominated surfaces and their extent is important to understand their role in the ecosystem.
The heterogeneous, mixed structure of dryland ecosystems is a challenge for biocrust mapping using remotely-sensed imagery (Smith et al., 2019). Previous studies have used unmanned aerial vehicles (UAVs) to investigate dryland vegetation (e.g., Cunliffe et al., 2016;Sankey et al., 2018;Milling et al., 2018), proving that this is an achievable task. However, remote sensing of biocrusts has traditionally been relegated to other platforms, such as airborne (e.g., Weber et al., 2008;Rodríguez-Caballero et al., 2014) and satellite sensors (e.g., Panigada et al., 2019). Only recently the potential of red, green and blue (RGB) imagery acquired from UAVs to identify dryland biocrusts has been explored (Havrilla et al., 2020;Jung et al., 2020;Sevgi et al., 2020). As biocrusts are difficult to distinguish from the background due to their small size and similar optical properties to vegetation and soil (Weber & Hill, 2016;Smith et al., 2019), UAVs integrated with multispectral cameras could greatly improve their identification. Multispectral sensors, in fact, usually have at least one spectral band in the red region at $660-680 nm, that corresponds to the chlorophyll-a absorption feature which is present in all chlorophytic biocrusts (Weber & Hill, 2016). Several authors (e.g., Blanco-Sacristán et al., 2019;Panigada et al., 2019;Rodríguez-Caballero, Escribano, et al., 2017a;Rodríguez-Caballero, Paul, et al., 2017b;Román et al., 2019) have exploited the subtle differences in this absorption features for biocrust identification using the continuum removal (CR) algorithm (Clark & Roush, 1984). This allows for the normalization of reflectance to a common baseline, which enables the analysis of individual absorption features.
UAVs have been commonly used to obtain digital surface models (DSMs) from RGB imagery using digital photogrammetry, from which digital terrain models (DTMs) of a very fine spatial resolution can be derived. This information can be used to calculate topography-related variables (Westoby et al., 2012). Detailed DTMs allow reseachers to apply hydrological models and therefore to study the impact of changing terrain properties on a landscape's hydrology (e.g., Lucieer et al., 2014). Thus, using these models it is possible to study the components of hydrological systems, such as surface and subsurface flows, which are key to understanding nutrient and sediment transport in the landscape (Stieglitz et al., 2003). However, site-specific relationships between soil physical properties and topography generate deviations from larger-scale climate patterns, which favours the formation of a heterogeneous community assembly (Rossi et al., 2014). In this context, as soil physical properties are a significant control on biocrust distribution in drylands Rodríguez-Caballero et al., 2019), it is important to explore the link between DTM-derived variables and soil properties. However, not all physical soil properties can be evaluated through DTMs. Key variables for biocrust development, such as soil moisture, need additional data sources in order to be estimated.
The knowledge of soil moisture distribution at a high spatiotemporal resolution plays an important role in the understanding of ecological and hydrological processes at the basin scale. Soil moisture plays a key role in surface water flow by controlling transport processes in the soil-plant-atmosphere system (Campbell & Norman, 1988 (Price, 1985;Short & Stuart, 1982). The dynamics of soil thermal properties and water content have also been analysed using proximal sensing (Krzeminska et al., 2012;Nearing et al., 2012;Negm et al., 2017). In arid and semi-arid environments, the spatial distribution of soil surface water content in bare soils has been evaluated from high resolution visible/ near-infrared (VIS/NIR) and thermal infrared (TIR) airborne and satellite data (Van Doninck et al., 2011) through the computation of the apparent thermal inertia (ATI). However, the coupling between ATI and soil moisture is not straightforward in heterogeneous surfaces since ATI might be directly related to soil moisture only in homogeneous areas with a single land cover type, with limitations observed for vegetated surfaces (Maltese et al., 2013;Van Doninck et al., 2011). For these reasons, the synergistic use of TIR and multispectral cameras on board UAVs allows the derivation of high spatial resolution maps of ATI, with the possibility to mask the vegetation component, helping to decipher the problem of mixed surfaces. Psora decipiens among others. A moss-dominated crust also develops, with species such as Pleurochaete squarrosa and Tortula revolvens.
Cyanobacteria genera, such as Microcoleus, Tolypothrix and Nostoc can also be found in the area (Cano-Díaz et al., 2018). Refer to Maestre et al. (2013) for a complete list of species of the visible biocrusts in the study area. We worked in two different areas within the study site ( Figure 1b): Area A, which covers an extension of 5.4 ha and is further sub-divided in two zones (Zones 1 and 2). Zone 1 with a moderate slope and a dominant west exposure and Zone 2 with a smaller slope and a dominant east exposure. And Area B, which covers 1.2 ha and is also sub-divided in two zones (Zones 3 and 4). Zone 3 presents an abrupt hill, mainly exposed to the south and Zone 4 a minor slope exposed to the north. was visually assessed and the species composition of biocrust and vascular vegetation was characterized. Pictures were also taken with a camera mounted on a tripod at 1 m height. Three representative soil sub-samples were taken in each plot from 0 to 10 cm depth (after removing the litter layer) and pooled to one sample, which represented the average soil characteristics of the plot. The samples were chemically and physically analysed in the laboratory applying standard techniques (Colombo & Miano, 2015). The samples were dried and sieved (2 mm mesh) and analysed to determine the percentage of soluble salts (SolSal), total phosphorus (P tot ), percentage of organic carbon (OrgCarb) and total nitrogen (N tot ). Soil moisture was measured in the field at 0-10 cm depth at three representative locations within each plot using a FieldScout TDR-100 soil moisture meter (Spectrum Technologies, Fort Worth, TX, USA).

| Unmanned aerial vehicle (UAV) flights
The UAV flights (Table 1)   to ten spectra were taken on each target and were averaged to obtain the final spectrum for each target. These ground spectral measurements were acquired around solar noon, contemporaneously to the multispectral images collected by the Parrot Sequoia camera (UAV flight 2). Their exact location was registered using a differential global navigation satellite system (dGNSS) including one Topcon HiPer Pro antennas and one Topcon GPR-1 (Topcon, Tokyo, Japan). The spectra of all ground measurements were resampled to the Sequoia spectral resolution using the response function of the sensor (Fawcett et al., 2020). These spectra were used to apply an empirical line correction (Smith & Milton, 1999) to the multispectral images to guarantee the radiometric quality in the final reflectance orthomosaics.
Simultaneous to the drone flights, the coordinates of 10 reference targets were taken in each area using the dGNSS. The targets were   (Ullman, 1983). Besides elevation, the slope of each cell was also calculated from the filtered DTM, applying the nine parameter second-order polynomial method (Zevenbergen & Thorne, 1987)  The contributing catchment area (CCA) of each pixel of the DTMs was modelled using the TauDEM toolbox (Tarboton, 1997;Tesfa et al., 2011) in ArcMap v.10.5. For this, the pits were first removed from the DTM and the D-Infinity Contributing Area function was used to calculate it. The D-Infinity algorithm (Tarboton, 1997) calculates the flow direction using the steepest downwards slope on the eight triangular facets formed in a 3 Â 3 pixel window centred on the grid cell of interest. The topographic wetness index (TWI; Beven & Kirkby, 1979), which quantifies the topographic control on hydrological processes, was calculated for each cell according to Equation 1: where A is the upstream area for each pixel and β is the slope in degrees.
The length slope factor (LSF; Renard et al., 1997), an indicator of the potential sediment transport or erosion risk under specific slope conditions, was calculated according to Equation 2: where LSF is the length slope factor, A the CCA, β the slope gradient and n and m are constant parameters set to 0.4 and 1.3, respectively.
The potential solar incoming radiation (PSIR) was calculated using the vegetation filtered DTMs of each area and the geometric solar radiation model implemented in ArcGIS Solar Analyst tool by assuming a uniform clear sky condition with a constant transmissivity of 0.5 (Fu & Rich, 1999) and a diffuse radiation proportion of 0.3. The map of the whole year PSIR (Wh/m 2 /yr) was used in this study as an indicator of the potential evapotranspiration (Monteith & Szeicz, 1962).

| Multispectral Parrot Sequoia images
The multispectral images were processed in Pix4DMapper (Pix4D S.A., Lausanne, Switzerland) to obtain the final reflectance for each spectral band of the image. First, the calibration coefficient K for converting measured pseudo-radiance R to reflectance was derived for each band as in Equation 3: : where γ is the ISO number, ε is the exposure time in seconds, f is the f-number and A, B and C are camera-specific calibration coefficients which model the non-linear behaviour of the sensor (Fawcett et al., 2020).
Finally, the reflectance per pixel (ρ) of the scene was derived as in Equation 5: To correct the reflectance measured by the multispectral sensor, the empirical line method (Smith & Milton, 1999) was applied. This method was based on an empirical relationship between the drone reflectance and ground reflectance of the reference targets (i.e., black and white reference panels) collected with the ASD Field Spectroradiometer in the field.
The white panel saturated on the Sequoia images, therefore homogeneous biocrust targets were used as bright targets instead.
The linear equation obtained for each multispectral band was applied to obtain the final reflectance of the multispectral images.
An orthomosaic was built from the multispectral images covering the study area with an original spatial resolution of 7.74 cm and 8.48 cm in Areas A and B, respectively, further resampled to 10 cm/pixel. The geometric accuracy of the orthomosaic was assessed by calculating the RMSE, the mean error and standard deviation error of the GCPs and the GVPs.
The CR algorithm was computed on the geometrically and radiometrically corrected multispectral orthomosaic. This method normalizes the reflectance spectrum to a common baseline at specified wavelengths and allows the comparison of absorption features in the spectra that are produced by target components, for example pigments for vegetation covers. This is achieved by approximating the continuum between local spectral maxima through straight-line segments; a value of 1 is assigned to the local maxima, and a value between 0 and 1 is obtained in correspondence of the absorption features. Interpolating the reflectance between 550 nm (Green band) and 735 nm (Red Edge band) as the continuum baseline.
Although this algorithm is more typically applied to hyperspectral data, it has already been used in multispectral imagery to map drylands  and vegetation parameters (Chauhan et al., 2020). The absorption feature related to chlorophyll-a was SVM is a statistical learning based supervised classification method (Cortes & Vapnik, 1995;Vapnik, 2005  The classification accuracy was assessed through confusion matrices (Story & Congalton, 1986). In the confusion matrix the ground truth data (i.e., columns) are compared to the classified data (i.e., rows). The major diagonal represents the agreement between these two data sets, and the overall accuracy (OA) of the classification is calculated by dividing the sum of the entries of this diagonal by the total number of samples taken. In addition, the Kappa coefficient (κ, Cohen, 1960) was calculated to asses the classification accuracy.
The accuracy of the classification over Area A and Area B was assessed with two validation datasets, each consisting in 10 polygons of nine pixels for each class (90 pixels per class). Polygons were manually defined using ancillary information such as field pictures, field notes, GPS points and high-resolution RGB orthomosaics.

| Thermal imagery and apparent thermal inertia (ATI) calculation
TIR images acquired by the thermoMAP camera were calibrated during the flight using an integrated shutter for in-flight radiometric calibration. This shutter automatically closes after every picture is captured and self-calibrates by comparing the grey level of each photograph with the temperature measured by the built-in temperature sensor of the camera. The raw images were converted using Pix4D software (Pix4D) to temperature in Celsius degrees according to Equation 6: where T is the absolute temperature in Celsius degrees and R is the radiometric value of thermoMap thermal images. The emissivity of field sampled biocrusts and soils was measured in the laboratory and found to be close to 1 for both targets (see Supporting Information, Table S1) in the range of the thermoMAP (8.5-13.5 μm). Since the aim of using these images in this study was to measure properties related to biocrusts, emissivity was therefore neglected and thus not included in this calculation.
Two orthomosaics (one captured before sunrise and one at noon) of each target area with a spatial resolution of 15 cm were generated using Pix4Dmapper (Pix4D). The geometric accuracy of these TIR orthomosaics was assessed by calculating the RMSE, the mean error and standard deviation of error of the GCPs and the GVPs.
Land surface temperature differences and albedo were used in the definition of the concept of apparent thermal inertia to assess the space variability of soil water content (Negm et al., 2017;Verstraeten et al., 2006). We applied the apparent thermal inertia (ATI; in K À1 ), which is an approximation of thermal inertia and is derived directly from multi-spectral remote sensing imagery (Mitra & Majumdar, 2004;Van Doninck et al., 2011) according to Equation 7: where ΔT is the amplitude of the diurnal temperature range calculated as the difference between the maximum and the minimum daily surface temperature, α is the surface spectral albedo and C is the solar correction factor that changes over space and time to normalize for solar flux variations with latitude and solar declination changes between seasons. The ΔT was calculated as the difference between the surface temperature of the orthomosaic captured at noon and the one captured before the sunrise. The value of C was calculated and had a value of 1.19. We approximated the albedo as in Equation 8: Vegetation pixels were excluded from our analysis. The highspatial resolution of TIR and multispectral images allowed us to use the classification of the multispectral images to create a vegetation-free mask that was applied to the ATI orthomosaics.

| Statistical analysis
To evaluate if there were statistically significant differences between the four study zones, a post hoc Dunn's test followed by a Kruskal-Wallis test was performed, using the dunnTest function from the FSA package (Ogle et al., 2021) in R. To analyse the effect of the terrain attributes in relation to the spatial distribution of soil properties measured in the sampling plots, the correlation ratio (Pearson, 1926) where y xi is the single soil property observation, x indicates the terrain attribute category, and i indicates an observation. If n x is the number of observations in the x terrain category, y x is the mean of the category x and y is the mean of the whole population.
The relationship between the different biocrust covers and the terrain variables in the whole study area was evaluated through a redundancy analysis (RDA) performed using the vegan R package (Oksanen et al., 2017)  GreenVeg) and soil fractional covers were also included as variables in this analysis. The dataset used in this analysis was obtained through a random sampling of 1147 points in the study area, with a minimum distance of 3 m between them to avoid spatial autocorrelation. This distance was calculated through a semi-variogram analysis (Curran, 1988

| Classification of the multispectral orthomosaics and evaluation of the classification accuracy
Mean reflectances and standard deviations of the training classes used for the classification procedure are shown in Figure 5 feature was present in all classes but Soil, and allowed to better differentiate it from biocrust classes. These subtle differences allowed an improvement in the classification carried out by the SVM, which slightly confused these two classes, by using thresholds in the CR red .
The SVM classification generated land cover maps of the two target areas in the study site ( Figure 6). The SVM classifications were improved using the CR red thresholds on the Soil and DryVeg classes. Final accuracies are reported in Table S2, OAs higher than 80% in most cases were achieved in both areas, with κ equal to 0.8 and 0.93 in areas A and B, respectively. Only BLM in Area A had an OA below 80%. As shown in the confusion matrices (Table S2), BL T A B L E 2 Root mean square error (RMSE), mean error (ME) and standard deviation of error (STDEV) for red, green and blue (RGB), multispectral and thermal infrared (TIR) imagery and both ground control point (GCP) and ground validation point (GVP)

| Thermal data and soil moisture
The soil moisture values measured in the 23 field plots were similar to that of the winter long-term soil moisture of this area , confirming that the moisture status we captured with the drone flight was representative of the winter rainy season. ATI values correlated well with the soil moisture (r 2 = 0.83; Figure 7), indicating that ATI is a good estimator of soil moisture, even when soils are covered by biocrusts. ATI was then used as a proxy of soil moisture for vegetation-free areas in the subsequent analysis.

| Analysis of soil properties
Soils in the study area are Gypsiric Leptosols and are characterized by the very high content of SalSol, reaching values higher than 90% in some of the plots, and the low presence of OrgCarb (Table 4) values of P tot of 0% and values of N tot close to 0% in some of them.
Moisture presented high variability among the plots, with mean values around 14%.

| Relationship between terrain attributes and soils at plot level
Terrain attributes in the field plots were found to vary across the four zones sampled (Figure 8). The elevation presented similar values in Zones 1 and 2 (Area A) and in Zones 3 and 4 (Area B). The slope was higher in Zone 3, reaching values of 35% of slope gradient. As expected, LSF presented similar patterns as slope. The TWI presented values around 3.5 with values not statistically different in the four zones. PSIR showed an opposite pattern compared to Northernness, with higher values in Zones 1 and 3, which are mainly exposed to southwest and south, respectively, indicating higher evapotranspiration. PSIR presented lower values in Zone 2 mainly exposed to the east and minimum values in Zone 4, mainly exposed to the north.
Similar to Northernness, ATI, that is, soil moisture, presented the highest values in Zone 4 and lower values in Zone 1, evidencing that the soil moisture is related to the aspect.
Among the terrain attributes investigated, Northernness, ATI and PSIR are the ones that explained most of the variance of the soil properties (Table 5). Soils more exposed to north and characterized by higher ATI show a lower content of SolSal and higher content of OrgCarb, P tot and N tot. Soils characterized by higher values of PSIR show an opposite trend, with higher content of SolSal and lower content of OrgCarb, P tot and N tot . The variance of soil properties is similarly explained by Northernness, ATI and PSIR with values of η 2 close to 0.4 when related to SalSol, OrgCarb, P tot and N tot , with slightly higher values of η 2 driven by PSIR.
It should be noted that UAV data collection took place during the wet season, when the water content in biocrusts is high and pigments are metabolically active (Weber & Hill, 2016). This enhances the differences in the spectral properties of biocrusts and improves their identification (Blanco-Sacristán et al., 2019). While the classifications from the multispectral images allowed mapping of the different environmental components of the study area, the values of ATI derived from TIR imagery correlated very highly with the soil moisture measured in the field. We demonstrated that the correlation between ATI and soil moisture is also maintained when biocrust cover is present.
To our knowledge, these are the first maps of ATI derived from UAV in drylands with presence of biocrust cover. Previous studies have used hyperspectral thermal sensors to discriminate different types of biocrusts from bare sand based on emissivity values and to characterize biocrusts' maturity (Rozenstein & Karnieli, 2015). In this study the ATI maps derived from TIR imagery were used to estimate how the soil moisture affects biocrust type and distribution, but ATI maps could also be used to better understand water distribution in drylands, especially when there is controversy about the final role of biocrust covers on soil water balance due to their influence on several hydrological processes and soil properties with antagonist effects (Chamizo et al., 2016). We found that topographic characteristics can partially explain the soil moisture gradient estimated by ATI (ATI and PSIR are inversely correlated with r 2 = 0.46). This moisture gradient promotes a diversification of ecological niches that also explains the distribution pattern of biocrust communities. However, long-term studies are still necessary to better understand the dynamics of soil moisture and the feedback effects on biocrust cover, and to try to decouple the biocrust effect on water from the water effect on biocrusts, which are strongly interrelated. ATI could be a useful variable to expand our knowledge in this direction.
The UAV RGB images were used to produce fine DTMs from which terrain attributes were derived and utilized to evaluate their impact on soil properties. The mean error of these DTMs was lower than 1 cm, which is comparable to previous studies applying SfM techniques for topography reconstruction (Bakker & Lane, 2017;Smith & Vericat, 2015). These results, together with the rough and steep topography of the study areas, suggest that these models were not affected by systematic error. We found that PSIR, ATI and Northernness were the terrain attributes mostly related to soil properties. These attributes explained in similar magnitudes the distribution of SolSal, OrgCarb, P tot and N tot . Increased shadows and higher soil moisture, indicated by low PSIR and high ATI, are related to a higher content of soil nutrients, as shown in plots located in Zones 2 and 4, predominantly east-and north-facing plots, respectively. Increased shadows and soil moisture control the survival and activity of microorganisms (Borken & Matzner, 2009;Drenovsky et al., 2004) and lead to better nutrient cycling and higher activity of microbial communities (e.g., Xue et al., 2018). This increased moisture favours not only the preservation of OrgCarb and its association with other mineral components (Plaza et al., 2013), but also the presence of vegetation and Together with the effect of vegetation on soil properties, biocrusts have been found to take up significant amounts of atmospheric C and N by photosynthesis and N fixation (Elbert et al., 2012).
They are therefore an important pool and source of organic inputs into the soil in drylands (Castillo-Monroy et al., 2011;Chamizo et al., 2012b;Concostrina-Zubiri et al., 2013). Furthermore, biocrusts strengthen soil structure by interacting with mineral particles and forming aggregates (Belnap & Lange, 2013;Eldridge & Leys, 2003), thus protecting soil from C loss. Biocrusts promote microbial community growth where they appear and the increased shadows created by vegetation increase this growth indirectly (Huang et al., 2015). However, the role of biocrusts on soil nutrients content depends on biocrust composition and patch-size distribution (e.g., Bowker et al., 2013;Delgado-Baquerizo et al., 2015;Sedia & Ehrenfeld, 2006 Variations in biocrust distribution and composition in dryland biomes are primarily controlled by climatic differences . However, vegetation and terrain properties generate variations from these climatic controls and modify the local distribution of nutrients in dryland landscapes (e.g., Manzoni et al., 2006;Puigdefábregas, 2005;Puigdefábregas & Sánchez, 1996). Thus, the different sensitivity of biocrusts to the distribution of resources conditions the composition of these ecosystems Durham et al., 2018;Rodríguez-Caballero et al., 2019). The RDA conducted using data from 1147 points showed that Fulgensia spp.
presented a strong positive relationship with unstable zones (high values of elevation and LSF). Fulgensia spp. has already been observed in unstable terrains several times, highlighting the pioneering behaviour of this genus (Cant on et al., 2020;Miralles et al., 2020;Rodríguez-Caballero et al., 2013). Mosses, typically found in shaded areas (i.e. Pleurochaete squarrosa and Syntrichia ruralis), dominated areas with low PSIR, high ATI, and a high presence of vascular vegetation, where soils were richer in organic carbon and nutrients. A positive plant-biocrust relationship is common for bryophytes (Zhou et al., 2020) and close correspondence with Macrochloa tenacissima presence has already been observed (Martínez-Sánchez et al., 1994).
The RDA showed a lower explanatory power of the terrain attributes regarding Tortula revolves, a moss that develops in arid environments, and BLM class dominated by the association of Diploschistes diacapsis and Squamarina lentigera (BLM). They appeared in more stable zones (i.e., low values of elevation and LSF), as already observed in geographic areas with similar characteristics (e.g., Ladr on de Guevara et al., 2018) that have high solar radiation and soils rich in soluble salts. In particular, Squamarina lentigera can be physiologically adapted to light-exposed environments, and requires high temperatures for optimal photosynthesis while being well hydrated (Lange et al., 1997). Nevertheless, Diploschistes diacapsis, the dominant species of BL, was also found in areas with high LSF. This species is very versatile and can adapt its physiology depending on where it grows (Pintado et al., 2005). We observed that Diploschistes diacapsis can develop directly on outcrop rocks with very few millimetres of soil beneath them in areas where gypsum substrates with fine soil texture favour stability, but where water availability is scarse. In these areas with outcrop rock, high concentrations of gypsum, soluble salts and higher solar radiation (i.e., high values of PSIR), Diploschistes diacapsis is more adapted to establish than other lichens, representing here the early stage of biocrust development.
We found instead mixed patches of lichens and moss (BLM) (thought to be of the last stages of development) in the left axis of the RDA plot, where water availability is still scarces but the terrain presentes lower elevation and LSF and less concentrations of gypsum. north-faced slopes (low evapotranspiration, low PSIR), while in our study area they appear in south-facing slopes (high evapotranspiration, high PSIR). Biocrust distribution in drylands might be more constrained by the effect of terrain attributes on local hydric availability, rather than rainfall water inputs. For this reason, it is difficult to make generalizations in the development of biocrusts' communities and developing monitoring methodologies that allow up-scaling local relationships of biocrusts with the terrain is key for dryland ecology.

| CONCLUSIONS
The synergistic information obtained from RGB, multispectral and thermal UAV imagery has allowed us to better understand the complex factors influencing the distribution of dryland ecosystem components. We produced accurate, high-resolution maps that can help to up-scale the local effect of biocrust components on ecosystem functioning by improving the accuracy of erosion and infiltration modelling. These aspects of dryland ecology are key to understanding the distinctive role that biocrust-dominated surfaces have on water runoff and erosion depending on the predominant type of biocrust (e.g., Rodríguez-Caballero et al., 2015;Wang et al., 2017). Producing high-resolution maps will help to monitor these communities in space and time, a key task to understand the compositional changes they are already experiencing due to climate change (Escolar et al., 2012; Ladr on de Guevara et al., 2018).
Since field-based data collection has many drawbacks, mainly related to time and costs (Palmer et al., 2002), exploiting UAV-based methodologies will help to build standardized procedures in dryland monitoring programmes, while providing very detailed information of these environments. These are needed because drylands may be characterized by different structural and functional organization and present a wide spectrum of compositions, depending on local climate and terrain properties. Our results have shown some discrepancies compared to the traditional biocrust development models, in which more advanced stages (i.e., mixed patches of lichens and moss) are assumed to appear only in stable soils. Thus, we suggest integrating more traditional approaches with new methodologies such as the one developed in this study, which can provide the very fine spatial resolution maps of dryland composition and terrain attributes needed for a detailed description and monitoring of these complex and threatened ecosystems.