Received: 9 November 2020 Revised: 16 June 2021 Accepted: 17 June 2021 DOI: 10.1002/esp.5189 RESEARCH ARTICLE UAV RGB, thermal infrared and multispectral imagery used to investigate the control of terrain on the spatial distribution of dryland biocrust Javier Blanco-Sacristán1 | Cinzia Panigada1 | Rodolfo Gentili1 | Giulia Tagliabue1 | Roberto Garzonio1 | M. Pilar Martín2 | Mo nica Ladro  n de Guevara3,4 | Roberto Colombo1 | Thomas P. F. Dowling5 | 1 Micol Rossini 1 Remote Sensing of Environmental Dynamics Lab, University of Milano-Bicocca, Milan, Italy Abstract 2 Environmental remote sensing and Biocrusts (topsoil communities formed by mosses, lichens, bacteria, fungi, algae, and spectroscopy laboratory (SpecLab), Spanish National Research Council (CSIC), Madrid, cyanobacteria) are a key biotic component of dryland ecosystems. Whilst climate Spain patterns control the distribution of biocrusts in drylands worldwide, terrain and soil 3 stoles, Spain Universidad Rey Juan Carlos, Mo attributes can influence biocrust distribution at landscape scale. Multi-source 4 Centre for Ecological Research and Forestry unmanned aerial vehicle (UAV) imagery was used to map and study biocrust ecology Applications, CREAF-CSIC-UAB, Barcelona, Spain in a typical dryland ecosystem in central Spain. Red, green and blue (RGB) imagery 5 United Nations Environment Programme was processed using structure-from-motion techniques to map terrain attributes World Conservation Monitoring Centre, Cambridge, UK related to microclimate and terrain stability. Multispectral imagery was used to pro- duce accurate maps (accuracy > 80%) of dryland ecosystem components (vegetation, Correspondence bare soil and biocrust composition). Finally, thermal infrared (TIR) and multispectral Cinzia Panigada, Remote Sensing of Environmental Dynamics Lab, University of imagery was used to calculate the apparent thermal inertia (ATI) of soil and to evalu- Milano-Bicocca, Milan, Italy. ate how ATI was related to soil moisture (r2 = 0.83). The relationship between soil Email: cinzia.panigada@unimib.it properties and UAV-derived variables was first evaluated at the field plot level. Then, Funding information the maps obtained were used to explore the relationship between biocrusts and ter- European Research Council, Grant/Award Number: 647038; Marie Skłodowska-Curie rain attributes at ecosystem level through a redundancy analysis. The most significant Actions, Grant/Award Number: 721995; MIUR variables that explain biocrust distribution are: ATI (34.4% of variance, F = 130.75; – Dipartimenti di Eccellenza 2018–2022; European Union’s Horizon 2020 p < 0.001), Elevation (25.8%, F = 97.6; p < 0.001), and potential solar incoming radia- tion (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. Twitter: Multi-source UAV imagery was used to study biocrust distribution in a typical dryland ecosystem in central Spain and how it is affected by terrain and soil properties. This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. © 2021 The Authors. Earth Surface Processes and Landforms published by John Wiley & Sons Ltd. 2466 wileyonlinelibrary.com/journal/esp Earth Surf. Process. Landforms. 2021;46:2466–2484.  ET AL. BLANCO-SACRISTAN 2467 KEYWORDS apparent inertia, thermal, biocrusts, biological soil crusts, drylands, lichen, moss, multispectral, thermal, UAV 1 | I N T RO DU CT I O N (e.g., Blanco-Sacristán et al., 2019; Panigada et al., 2019; Rodríguez- Caballero, Escribano, et al., 2017a; Rodríguez-Caballero, Paul, Drylands constitute one of the largest biomes on Earth, covering et al., 2017b; Román et al., 2019) have exploited the subtle differ- 47% of the terrestrial surface (Koutroulis, 2019). In these ences in this absorption features for biocrust identification using the environments, topography determines the redistribution of scarce continuum removal (CR) algorithm (Clark & Roush, 1984). This allows precipitation, controls water content and the availability of soil nutri- for the normalization of reflectance to a common baseline, which ents and organic matter affected by runoff (e.g., Aguiar & Sala, 1999; enables the analysis of individual absorption features. Manzoni et al., 2006; Puigdefábregas & Sánchez, 1996; UAVs have been commonly used to obtain digital surface models Puigdefábregas et al., 1999; Puigdefábregas, 2005). This modifies not (DSMs) from RGB imagery using digital photogrammetry, from which only vegetation distribution, but also the allocation of components digital terrain models (DTMs) of a very fine spatial resolution can be that appear in plant interspaces, such as biocrust communities derived. This information can be used to calculate topography-related (Rodríguez-Caballero et al., 2019). Biocrusts are a combination of top- variables (Westoby et al., 2012). Detailed DTMs allow reseachers to soil communities including mosses, lichens, liverworts, bacteria, fungi, apply hydrological models and therefore to study the impact of algae and cyanobacteria. They play an important role in nutrient changing terrain properties on a landscape’s hydrology (e.g., Lucieer cycling (Elbert et al., 2012; Weber et al., 2015), soil carbon (C) fluxes et al., 2014). Thus, using these models it is possible to study the (Castillo-Monroy et al., 2011; Tucker et al., 2019) and runoff/runon components of hydrological systems, such as surface and subsurface dynamics (Chamizo et al., 2012a; Faist et al., 2017; Rodríguez- flows, which are key to understanding nutrient and sediment trans- Caballero et al., 2015), playing a key role in dryland ecosystems’ ser- port in the landscape (Stieglitz et al., 2003). However, site-specific vice maintenance (Rodríguez-Caballero, Escribano, et al., 2017a; relationships between soil physical properties and topography Rodríguez-Caballero, Paul, et al., 2017b). The effect that these generate deviations from larger-scale climate patterns, which favours communities have on dryland surfaces depends on the dominant and the formation of a heterogeneous community assembly (Rossi developmental stage of the crust (Belnap & Lange, 2013; et al., 2014). In this context, as soil physical properties are a significant Concostrina-Zubiri et al., 2013; Faist et al., 2017; Tucker et al., 2019). control on biocrust distribution in drylands (Bowker et al., 2016; For example, biocrusts can increase water availability for plants by Rodríguez-Caballero et al., 2019), it is important to explore the link augmenting water retention in topsoil (Eldridge et al., 2020) and between DTM-derived variables and soil properties. However, not all reducing soil evaporation (e.g., Adessi et al., 2018; Chamizo physical soil properties can be evaluated through DTMs. Key variables et al., 2016). They can also modify the erosion rate (Gao et al., 2020), for biocrust development, such as soil moisture, need additional data which affects sediment accumulation, C and nutrient content in soils sources in order to be estimated. (e.g., Chamizo et al., 2017) and surface roughness (e.g., Rodríguez- The knowledge of soil moisture distribution at a high spatio- Caballero et al., 2015; Wang et al., 2017). For this reason, mapping temporal resolution plays an important role in the understanding of the spatial patterns of different biocrust-dominated surfaces and their ecological and hydrological processes at the basin scale. Soil moisture extent is important to understand their role in the ecosystem. plays a key role in surface water flow by controlling transport The heterogeneous, mixed structure of dryland ecosystems is a processes in the soil–plant–atmosphere system (Campbell & challenge for biocrust mapping using remotely-sensed imagery (Smith Norman, 1988). However, estimation of the soil moisture involves in et al., 2019). Previous studies have used unmanned aerial vehicles situ measurements by standard point-based techniques (UAVs) to investigate dryland vegetation (e.g., Cunliffe et al., 2016; (e.g., thermogravimetric method or time domain reflectometry) and Sankey et al., 2018; Milling et al., 2018), proving that this is an achiev- requires interpolation techniques to obtain spatially explicit able task. However, remote sensing of biocrusts has traditionally been information. Different methods have been applied to analyse the spa- relegated to other platforms, such as airborne (e.g., Weber tial distribution of soil moisture based on thermal inertia, a physical et al., 2008; Rodríguez-Caballero et al., 2014) and satellite sensors property of the surfaces, which determines resistance to temperature (e.g., Panigada et al., 2019). Only recently the potential of red, green change under seasonal and diurnal heating (Price, 1985; Short & and blue (RGB) imagery acquired from UAVs to identify dryland Stuart, 1982). The dynamics of soil thermal properties and water biocrusts has been explored (Havrilla et al., 2020; Jung et al., 2020; content have also been analysed using proximal sensing (Krzeminska Sevgi et al., 2020). As biocrusts are difficult to distinguish from the et al., 2012; Nearing et al., 2012; Negm et al., 2017). In arid and background due to their small size and similar optical properties to semi-arid environments, the spatial distribution of soil surface water vegetation and soil (Weber & Hill, 2016; Smith et al., 2019), UAVs content in bare soils has been evaluated from high resolution visible/ integrated with multispectral cameras could greatly improve their near-infrared (VIS/NIR) and thermal infrared (TIR) airborne and satel- identification. Multispectral sensors, in fact, usually have at least one lite data (Van Doninck et al., 2011) through the computation of the spectral band in the red region at 660–680 nm, that corresponds to apparent thermal inertia (ATI). However, the coupling between ATI the chlorophyll-a absorption feature which is present in all and soil moisture is not straightforward in heterogeneous surfaces chlorophytic biocrusts (Weber & Hill, 2016). Several authors since ATI might be directly related to soil moisture only in 2468  ET AL. BLANCO-SACRISTAN homogeneous areas with a single land cover type, with limitations 2 | M A T E R I A L S A N D M ET H O D S observed for vegetated surfaces (Maltese et al., 2013; Van Doninck et al., 2011). For these reasons, the synergistic use of TIR and multi- 2.1 | Study area spectral cameras on board UAVs allows the derivation of high spatial resolution maps of ATI, with the possibility to mask the vegetation This study was conducted in a typical dryland ecosystem at the component, helping to decipher the problem of mixed surfaces. Aranjuez Experimental Station, located in central Spain (40 020 N, The objective of this study was to evaluate the effect of terrain  n de Guevara 3 320 W; 590 m above sea level, Figure 1) (Ladro attributes in the distribution of dryland ecosystem components’ using et al., 2018). The climate in this area is semi-arid Mediterranean, with multisensor UAV-based imagery by( i) mapping dryland constituents a mean annual temperature of 15 C and mean annual rainfall of (i.e., vegetation, biocrusts, and bare soil) through multispectral imag- 349 mm. The soil is classified as Gypsiric Leptosol (IUSS Working ery; (ii) estimating terrain attributes from structure-from-motion (SfM) Group WRB, 2006). The vascular vegetation cover in this area is techniques applied to RGB imagery and evaluating their relation with < 40% and is dominated by Macrochloa tenacissima tussocks and, in a soil properties; (iii) estimating soil moisture using maps of ATI derived lesser amount, small shrubs such as Helianthemum squamatum and from TIR imagery; (iv) evaluating through a multivariate statistical Gypsophila struthium. On the surface area that is not covered by approach how terrain attributes and soil moisture affect biocrust vascular vegetation, a rich biocrust community dominated by lichens distribution in the study area. develops, including Diploschistes diacapsis, Squamarina lentigera and F I G U R E 1 (a) Location of the study area in Aranjuez, central Spain; (b) study Area A and Zones 1 and 2; (c) study Area B and Zones 3 and 4. Black squares represent the ground control points (GCPs) used for the generation of the digital terrain model (DTM) and the orthomosaics and black circles are the ground validation points (GVPs) used to carry out an independent validation of the generated models [Color figure can be viewed at wileyonlinelibrary.com]  ET AL. BLANCO-SACRISTAN 2469 Psora decipiens among others. A moss-dominated crust also develops, a camera mounted on a tripod at 1 m height. Three representative soil with species such as Pleurochaete squarrosa and Tortula revolvens. sub-samples were taken in each plot from 0 to 10 cm depth (after Cyanobacteria genera, such as Microcoleus, Tolypothrix and Nostoc can removing the litter layer) and pooled to one sample, which represen- also be found in the area (Cano-Díaz et al., 2018). Refer to Maestre ted the average soil characteristics of the plot. The samples were et al. (2013) for a complete list of species of the visible biocrusts in chemically and physically analysed in the laboratory applying standard the study area. We worked in two different areas within the study site techniques (Colombo & Miano, 2015). The samples were dried and (Figure 1b): Area A, which covers an extension of 5.4 ha and is further sieved (2 mm mesh) and analysed to determine the percentage of sub-divided in two zones (Zones 1 and 2). Zone 1 with a moderate soluble salts (SolSal), total phosphorus (Ptot), percentage of organic slope and a dominant west exposure and Zone 2 with a smaller slope carbon (OrgCarb) and total nitrogen (Ntot). Soil moisture was and a dominant east exposure. And Area B, which covers 1.2 ha and is measured in the field at 0–10 cm depth at three representative also sub-divided in two zones (Zones 3 and 4). Zone 3 presents an locations within each plot using a FieldScout TDR-100 soil moisture abrupt hill, mainly exposed to the south and Zone 4 a minor slope meter (Spectrum Technologies, Fort Worth, TX, USA). exposed to the north. 2.2.2 | Unmanned aerial vehicle (UAV) flights 2.2 | Field campaign The UAV flights (Table 1) were conducted using an eBee SQ In March 2018 a field campaign took place to characterize the biocrust, (senseFly, Cheseaux-sur-Lausanne, Switzerland) fixed-wing UAV oper- vascular vegetation cover and soils. We collected soil moisture samples ated by the company 3D Scanner (Zaragoza, Spain), on 22 March and characterized the biocrust and vascular vegetation of the study 2018 under clear-sky conditions. These flights were conducted at the area as described in Section 2.2.1. RGB, multispectral and TIR UAV- end of the rainy season in the study area because biocrusts are more based imagery was concurrently acquired to map the study area, derive colourful when are well hydrated by the rain. When colourful they are terrain attributes and evaluate soil moisture. A general workflow with easier to map and differentiate using spectral information. In the main actions of this study is represented in Figure 2. Supporting Information Figure S1, the daily precipitation as well as the minimum and maximum temperatures registered in the study area during 2018 are reported. Three different spectral cameras, i.e., S.O. 2.2.1 | Biocrust characterization and soil sampling D.A (senseFly), Parrot Sequoia (Parrot SA, Paris, France) and ther- moMAP (senseFly) were flown to acquire RGB, multispectral and TIR On 21 March 2018, we established 23 2  2 m2 observational field images, respectively. Two sets of TIR imagery were acquired: the first plots within the two areas. Within each plot the relative abundance of before the sunrise, and the second at noon. During this second flight each type of surface component (biocrust, bare soil and vegetation) multispectral imagery was simultaneously acquired. A third flight, was visually assessed and the species composition of biocrust and conducted immediately after the previous one, was performed to vascular vegetation was characterized. Pictures were also taken with acquire RGB images on the site. F I G U R E 2 Workflow with the main processes carried out in this study. UAV, unmanned aerial vehicle; TIR, thermal infrared; RGB, red, green and blue; SVM, support vector machine; ATI, apparent thermal inertia; DTM, digital terrain model; RDA, redundancy analysis 2470  ET AL. BLANCO-SACRISTAN T A B L E 1 Technical information of the unmanned aerial vehicle (UAV) flights using the thermoMAP, Parrot Sequoia and S.O.D.A. sensors over the Areas A and B Area Start time Approximate Number of Average altitude above Average ground Flight Sensor flown (UTC) duration (min) images collected ground level (m) sampling distance (cm) 1 thermoMAP A 06:13 13 2397 69 13.09 2 thermoMAP B 06:28 5 660 84 15.00 3 Sequoia A 11:45 19 297 82 7.74 4 Sequoia B 12:06 9 115 90 8.48 5 thermoMAP A 11:52 13 2211 69 13.09 6 thermoMAP B 12:06 5 597 84 15.00 7 S.O.D.A. A 13:17 19 361 66 1.55 8 S.O.D.A. B 13:06 9 135 76 1.79 All flights followed a single grid flight geometry at constant to ten spectra were taken on each target and were averaged to obtain height, with flight lines parallel to the slope and with the cameras in the final spectrum for each target. These ground spectral measure- nadir-view position. Flights were done at average flight speeds ments were acquired around solar noon, contemporaneously to the between 10 and 12 m/s. Lateral and longitudinal overlaps of at least multispectral images collected by the Parrot Sequoia camera (UAV 75% and 70%, respectively, were achieved between the RGB images, flight 2). Their exact location was registered using a differential global with an average distance of 22 m between the flight lines and 20 m navigation satellite system (dGNSS) including one Topcon HiPer Pro between images. Lateral and longitudinal overlaps of at least 85% and antennas and one Topcon GPR-1 (Topcon, Tokyo, Japan). The spectra 70%, respectively, were obtained between the multispectral images, of all ground measurements were resampled to the Sequoia spectral with a distance of 12 m between the flight lines and 24 m between resolution using the response function of the sensor (Fawcett et al., images. Lateral and longitudinal overlaps of at least 80% and 95%, 2020). These spectra were used to apply an empirical line correction respectively, were achieved between the TIR images, with a distance (Smith & Milton, 1999) to the multispectral images to guarantee the of 16 m between the flight lines and 3 m between images. The multi- radiometric quality in the final reflectance orthomosaics. spectral Parrot Sequoia camera is a compact bundle of four cameras Simultaneous to the drone flights, the coordinates of 10 reference with complementary metal-oxide-semiconductor sensors, each one targets were taken in each area using the dGNSS. The targets were equipped with an individual narrow-band filter (Assmann et al., 2018). deployed in clear areas as black and white grids printed on thick paper This results in a multi-band image composed of four bands: Green sheets (297 mm  420 mm). To optimize the geometric accuracy of (550  40 nm), Red (660  40 nm), Red Edge (735  10 nm) and NIR the UAV-derived products in both planimetry (XY direction) and (790  10 nm). The Red Edge band is located in the spectral region altimetry (Z direction), the targets were evenly distributed across the characterized by the rapid change in vegetation reflectance between whole study area of each zone, since the target distribution might red absorption due to pigment content and NIR reflectance due to the affect later model generation accuracy (Gindraux et al., 2017; ® leaf structure. A white reference Micasense panel was measured Martínez-Carricondo et al., 2018) (Figure 1). In each area, six of the with the Sequoia camera immediately before and directly after the targets were used as ground control points (GCPs) for the generation drone survey to calibrate the images. of the orthophoto and DTM, while the remaining four were used as ground validation points (GVPs). Giving that the size of Area A and Area B were around 5.4 ha and 1.2 ha, respectively, the density of 2.2.3 | Ground data targets in the study site (111 and 600 GCP/km2 for Areas A and B, respectively) was enough to guarantee high accuracy in the derived Reflectance factors were collected in the field to calibrate and validate digital models (Gindraux et al., 2017; Martínez-Carricondo the spectral data acquired from the UAV sensors using an ASD et al., 2018). We used a rapid-static measurement technique and the FieldSpec® FR3 field spectroradiometer (Analytical Spectral Devices, coordinates were corrected during post-processing (Hofmann- Longmont, CO, USA) which features a spectral range from 350 nm to Wellenhof et al., 2007). One of the two GNSS devices was used as 2500 nm, an optical resolution of 1 nm and a nominal field of view of master and the other as rover. The master station was placed in a  25 . Two types of reference targets were measured during the flight clear flat area for a cumulative time of 5 h and was used as a static campaign: white and black 2  2 m2 targets made of polyvinyl chloride reference point. The rover was used to measure the centre point of (PVC)-coated canvas material (Odyssey trademark material; each GCP and GVP, with the rover keept over the point for at least Kayospruce Ltd, Fareham, UK) assumed to be lambertian surfaces, and 2 min, taking a measurement every second with at least six satellites selected homogeneous plots from the main biocrust communities in in the field of view of the receiver, according to the protocol proposed the study site. These measurements covered homogeneous surfaces by Hegarty and Chatre (2008). The master and rover data were post- of at least 30 cm in diameter. Before each target measurement, the processed using Topcon Tools software (Topcon Positioning Systems, downwelling solar radiance was measured on a calibrated 99% reflec- Inc.). The dGNSS system used during the surveys has a horizontal ® tive Spectralon panel (Labsphere Inc., Sutton, NH, USA) to optimize precision of 3 mm + 0.5 ppm and a vertical precision of 5 mm the instrument parameters and calculate the reflectance factors. Eight + 0.5 ppm according to the Topcon HiPer Pro specification data.  ET AL. BLANCO-SACRISTAN 2471 2.3 | Unmanned aerial vehicle (UAV) data each cell of the DTM to the North) was calculated as the cosine of processing each cell of the aspect model. Since this variable ranges from 1 to 1, +1 was added to each pixel value to avoid problems in later 2.3.1 | Red, green and blue (RGB) images: calculations. Orthophoto, digital terrain model (DTM) and terrain The contributing catchment area (CCA) of each pixel of the DTMs attributes was modelled using the TauDEM toolbox (Tarboton, 1997; Tesfa et al., 2011) in ArcMap v.10.5. For this, the pits were first removed RGB images acquired over Areas A and B of the study site by the from the DTM and the D-Infinity Contributing Area function was used S.O.D.A. camera were processed into orthophotos and grid-based to calculate it. The D-Infinity algorithm (Tarboton, 1997) calculates DSMs following a SfM workflow (Lucieer et al., 2014; Westoby the flow direction using the steepest downwards slope on the eight et al., 2012) in Agisoft PhotoScan v.1.4 (Agisoft, St. Petersburg, triangular facets formed in a 3  3 pixel window centred on the grid Russia). A detailed description of the SfM algorithms used in Agisoft cell of interest. The topographic wetness index (TWI; Beven & PhotoScan can be found in Verhoeven (2011). The images captured Kirkby, 1979), which quantifies the topographic control on hydrologi- during the flights were aligned using an image feature recognition cal processes, was calculated for each cell according to Equation 1: algorithm similar to Lowe’s (2004) scale-invariant feature transform   method, which automatically detects and matches unique image A TWI ¼ ln ð1Þ feature characteristics that are stable under variations in view per- tanβ spective and illumination across input photographs. This alignment results in a sparse three-dimensional (3D) point cloud, which is used where A is the upstream area for each pixel and β is the slope in to create a dense point cloud by an iterative bundle adjustment algo- degrees. rithm. In turn this recreates its 3D geometry and camera positions The length slope factor (LSF; Renard et al., 1997), an indicator of from a sequence of two-dimensional images acquired from multiple the potential sediment transport or erosion risk under specific slope viewpoints (Ullman, 1983). The onboard navigation sensors allowed conditions, was calculated according to Equation 2: the camera positions and 3D point cloud to be automatically  n   georeferenced within the precision of the instrument. The GCPs were A sinβ m LSF ¼ ðn þ 1Þ   ð2Þ then manually identified on the images, and their dGNSS coordinates 22:13 0:0896 were imported to optimize the spatial accuracy and improve the geometry of the 3D point cloud. A multi-view stereo image matching where LSF is the length slope factor, A the CCA, β the slope gradient algorithm was applied to the point clouds to increase the density of and n and m are constant parameters set to 0.4 and 1.3, respectively. the point cloud and to convert them into the orthomosaics of each The potential solar incoming radiation (PSIR) was calculated using area. The GVPs were also identified in the images and they were only the vegetation filtered DTMs of each area and the geometric solar used for the model accuracy assessment in post-processing and not radiation model implemented in ArcGIS Solar Analyst tool by assum- for their creation. The geometric accuracy of the models was assessed ing a uniform clear sky condition with a constant transmissivity of 0.5 by calculating the root mean square error (RMSE), the mean error and (Fu & Rich, 1999) and a diffuse radiation proportion of 0.3. The map standard deviation of error of the GCPs and the GVPs. of the whole year PSIR (Wh/m2/yr) was used in this study as an indi- To generate free-of-vegetation DTMs, the points belonging to cator of the potential evapotranspiration (Monteith & Szeicz, 1962). vegetation were removed from the point cloud (Figure 3) using the CANUPO plug-in (Brodu & Lague, 2012) in CloudCompare. For each area, two-point clouds representing each class (i.e., vegetation and 2.3.2 | Multispectral Parrot Sequoia images bare soil) were manually segmented from the point cloud. These sub- clouds were used for training the classifier, computing a dimensional- The multispectral images were processed in Pix4DMapper (Pix4D ity descriptor on the original point cloud with a regular ramp of scale S.A., Lausanne, Switzerland) to obtain the final reflectance for each values (i.e., sampled from 0.1 to 1 with steps of 0.1). The classifier was spectral band of the image. First, the calibration coefficient K for con- applied to the point cloud (using selected core points) classifying and verting measured pseudo-radiance R to reflectance was derived for removing vegetation points. Finally, isolated points were removed by each band as in Equation 3: using the noise filter tool (i.e., specifying a radius of 0.3 m and a rela- ρref tive error) and the DTMs were exported from CloudCompare with a K¼ ð3Þ Rref resolution of 10 cm, interpolating empty cells with average values. Free-of-vegetation DTMs were smoothed in ArcMap by applying an average filter with the Focal statistics tool by considering a circular where ρref is the known reflectance of the panel for the Sequoia band window of 1 m diameter. This filter was applied to minimize potential and Rref represents the measured pseudo-radiance averaged over the errors that can occur with the removing vegetation procedure. pixels of the reference panel. Values for ρref were provided by Besides elevation, the slope of each cell was also calculated from the manufacturer. The digital number (P) of each pixel were converted the filtered DTM, applying the nine parameter second-order to pseudo-radiance in arbitrary units R according to Equation 4: polynomial method (Zevenbergen & Thorne, 1987) within the Slope function in QGIS v.2.18.12. The aspect was derived using the Aspect P-B R ¼ f2 ð4Þ function in QGIS v.2.18.12. The Northernness (i.e. the exposure of Aγε þ C 2472  ET AL. BLANCO-SACRISTAN F I G U R E 3 Example of vegetation classification derived from the point cloud used to derive the vegetation free digital terrain model (DTM) in Area A of the study site. Left: Red, green and blue (RGB) true colour composite with vegetation and bare soil areas. Right: Classification of vegetation (red) and bare soil (blue) [Color figure can be viewed at wileyonlinelibrary.com] where γ is the ISO number, ε is the exposure time in seconds, f is the drylands (Panigada et al., 2019) and vegetation parameters (Chauhan f-number and A, B and C are camera-specific calibration coefficients et al., 2020). The absorption feature related to chlorophyll-a was which model the non-linear behaviour of the sensor (Fawcett computed at 660 nm by extracting the value of the CR in the red band et al., 2020). (CRred). The CRred was used to improve later classification of Finally, the reflectance per pixel (ρ) of the scene was derived as in the images. Equation 5: Seven different surface covers were identified in the study site during the field surveys (Figure 4): bare soil (Soil), bright lichens ρ ¼ KR ð5Þ (BL; patches of Squamarina spp., Diploschistes spp., Buellia spp.), bright lichens with moss (BLM; Squamarina spp., Diploschistes spp., Buellia To correct the reflectance measured by the multispectral sensor, spp. and Tortula revolvens [Zone 1] or Syntrichia ruralis [Zone 2]), moss the empirical line method (Smith & Milton, 1999) was applied. This (Moss; Tortula revolvens [Zone 1] and Syntrichia ruralis [Zone 2]), method was based on an empirical relationship between the drone Fulgensia spp. and moss (Fulg; Fulgensia spp. and Tortula revolvens reflectance and ground reflectance of the reference targets (i.e., black [Zone 3] or Pleurochaete squarrosa [Zone 4]), green vegetation and white reference panels) collected with the ASD Field (GreenVeg; green tussocks of Macrochloa tenacissima) and dry vegeta- Spectroradiometer in the field. tion (DryVeg; dry individuals of Helianthemum squamatum, Gypsophila The white panel saturated on the Sequoia images, therefore struthium and Macrochloa tenacissima). These surface covers were homogeneous biocrust targets were used as bright targets instead. used as reference for the digital classification of the multispectral The linear equation obtained for each multispectral band was applied orthomosaic. to obtain the final reflectance of the multispectral images. At least 15 polygons (around 100 pixels) per class were selected An orthomosaic was built from the multispectral images covering on each orthomosaic to train the support vector machine (SVM) the study area with an original spatial resolution of 7.74 cm and algorithm used in the classification. These training areas were 8.48 cm in Areas A and B, respectively, further resampled to selected using ancillary information such as field notes, global 10 cm/pixel. The geometric accuracy of the orthomosaic was assessed positioning system (GPS) points and high resolution RGB by calculating the RMSE, the mean error and standard deviation error orthomosaics. The four bands of the multispectral images were used of the GCPs and the GVPs. as input in the SVM. The CR algorithm was computed on the geometrically and SVM is a statistical learning based supervised classification radiometrically corrected multispectral orthomosaic. This method method (Cortes & Vapnik, 1995; Vapnik, 2005). It uses training normalizes the reflectance spectrum to a common baseline at samples of the target categories and separates them using a decision specified wavelengths and allows the comparison of absorption surface called the ‘hyperplane’, which maximizes the margin between features in the spectra that are produced by target components, for each category. The support vectors are the closest training samples to example pigments for vegetation covers. This is achieved by approxi- this hyperplane, and are the ones used by the algorithm. In this study, mating the continuum between local spectral maxima through a radial basis kernel function was used. The SVM was used to classify straight-line segments; a value of 1 is assigned to the local maxima, the images and derive thematic maps with the spatial distribution of and a value between 0 and 1 is obtained in correspondence of the the main land covers in the two areas. CRred was used to improve the absorption features. Interpolating the reflectance between 550 nm discrimination of some categories with similar spectral behaviour such (Green band) and 735 nm (Red Edge band) as the continuum baseline. as DryVeg and Moss, but different pigment absorption. Pixels classi- Although this algorithm is more typically applied to hyperspectral fied as DryVeg with CRred values lower than 0.75 were assigned to data, it has already been used in multispectral imagery to map Moss and pixels of BL with CRred values of 1 were assigned to Soil.  ET AL. BLANCO-SACRISTAN 2473 F I G U R E 4 Orthomosaic of Area B and details of the classes used to classify the multispectral images. BL, bright lichens; BLM, bright lichens and moss; Fulg, Fulgensia spp. and moss; GreenVeg, green vegetation; DryVeg, dry vegetation [Color figure can be viewed at wileyonlinelibrary.com] These thresholds were defined from the CRred mean and standard T ¼ 0:01 R  100 ð6Þ deviation values calculated in correspondence with the training polygons of soil, DryVeg and Moss classes. where T is the absolute temperature in Celsius degrees and R is the The classification accuracy was assessed through confusion matri- radiometric value of thermoMap thermal images. The emissivity of ces (Story & Congalton, 1986). In the confusion matrix the ground field sampled biocrusts and soils was measured in the laboratory and truth data (i.e., columns) are compared to the classified data found to be close to 1 for both targets (see Supporting Information, (i.e., rows). The major diagonal represents the agreement between Table S1) in the range of the thermoMAP (8.5–13.5 μm). Since the these two data sets, and the overall accuracy (OA) of the classification aim of using these images in this study was to measure properties is calculated by dividing the sum of the entries of this diagonal by the related to biocrusts, emissivity was therefore neglected and thus not total number of samples taken. In addition, the Kappa coefficient included in this calculation. (κ, Cohen, 1960) was calculated to asses the classification accuracy. Two orthomosaics (one captured before sunrise and one at noon) The accuracy of the classification over Area A and Area B was of each target area with a spatial resolution of 15 cm were generated assessed with two validation datasets, each consisting in 10 polygons using Pix4Dmapper (Pix4D). The geometric accuracy of these TIR of nine pixels for each class (90 pixels per class). Polygons were manu- orthomosaics was assessed by calculating the RMSE, the mean error ally defined using ancillary information such as field pictures, field and standard deviation of error of the GCPs and the GVPs. notes, GPS points and high-resolution RGB orthomosaics. 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 2.3.3 | Thermal imagery and apparent thermal et al., 2006). We applied the apparent thermal inertia (ATI; in K1), inertia (ATI) calculation which is an approximation of thermal inertia and is derived directly from multi-spectral remote sensing imagery (Mitra & Majumdar, 2004; TIR images acquired by the thermoMAP camera were calibrated Van Doninck et al., 2011) according to Equation 7: during the flight using an integrated shutter for in-flight radiometric calibration. This shutter automatically closes after every picture is 1α ATI ¼ C ð7Þ ΔT captured and self-calibrates by comparing the grey level of each pho- tograph with the temperature measured by the built-in temperature sensor of the camera. The raw images were converted using Pix4D where ΔT is the amplitude of the diurnal temperature range calculated software (Pix4D) to temperature in Celsius degrees according to as the difference between the maximum and the minimum daily sur- Equation 6: face temperature, α is the surface spectral albedo and C is the solar 2474  ET AL. BLANCO-SACRISTAN correction factor that changes over space and time to normalize for of observations in the x terrain category, yx is the mean of the cate- solar flux variations with latitude and solar declination changes gory x and y is the mean of the whole population. between seasons. The ΔT was calculated as the difference between The relationship between the different biocrust covers and the the surface temperature of the orthomosaic captured at noon and the terrain variables in the whole study area was evaluated through a one captured before the sunrise. The value of C was calculated and redundancy analysis (RDA) performed using the vegan R package had a value of 1.19. We approximated the albedo as in Equation 8: (Oksanen et al., 2017) v.2.4. To test the significance of the selected variables in the RDA, permutation tests (N = 999) were performed using the marginal effect of the terms in the anova.cca ρGREEN þ ρRED þ ρREG þ ρNIR α¼ ð8Þ function. 4 The CCA variable was removed from the analysis since it was used to calculate TWI and LSF. Based on this, a set of seven terrain where ρGREEN, ρRED, ρREG and ρNIR are the reflectance values in attributes were considered in the analysis: (i) PSIR as an indicator of the Green, Red, Red Edge and NIR bands, respectively. Multispectral the potential evapotranspiration (Monteith & Szeicz, 1962) and as a bands were calibrated so that the brightness map can reasonably be proxy of thermal microclimate (Durham et al., 2018; Suggitt considered to be a proxy of the broadband visible albedo. The et al., 2018); (ii) LSF as an index of surface stability (Renard et al., correlation between thermal inertia and soil moisture, observed in the 1997); (iii) TWI as an index of the topographic control on hydrological literature (Minacapilli et al., 2009), was here tested through a linear processes (Sörensen et al., 2006); (iv) elevation; (v) Slope gradient, as regression model between the mean ATI value extracted from it has implications for dryland components distribution not orthomosaics and the mean soil moisture measured on each plot captured by other variables (Rodríguez-Caballero et al., 2019); during the field campaign by traditional TDR methods. The average (vi) Northernness, as indicator of facing of the slope and (vii) ATI, as a value of ATI on each 2  2 m2 plot was obtained in QGIS v.2.18.12 proxy of soil moisture. The vegetation (combination of DryVeg and by averaging the ATI values of the corresponding polygons. GreenVeg) and soil fractional covers were also included as variables in Vegetation pixels were excluded from our analysis. The high- this analysis. The dataset used in this analysis was obtained through a spatial resolution of TIR and multispectral images allowed us to use random sampling of 1147 points in the study area, with a minimum the classification of the multispectral images to create a distance of 3 m between them to avoid spatial autocorrelation. This vegetation-free mask that was applied to the ATI orthomosaics. distance was calculated through a semi-variogram analysis (Curran, 1988). A circle of 3 m was drawn around each point to extract the terrain attributes and the fractional cover of the different components 2.4 | Statistical analysis (i.e., soil, vegetation and different biocrust compositions). Several models with all combinations of variables were computed to identify To evaluate if there were statistically significant differences between the best-fitting RDA model based on maximum adjusted r2, but the four study zones, a post hoc Dunn’s test followed by a Kruskal– avoiding multicollinearity. This was estimated by the variance inflation Wallis test was performed, using the dunnTest function from the FSA factor (VIF), which was calculated through the vif.cca function of package (Ogle et al., 2021) in R. To analyse the effect of the terrain vegan R package. Variables with a VIF higher than 10 were excluded attributes in relation to the spatial distribution of soil properties mea- from the analysis. sured in the sampling plots, the correlation ratio (Pearson, 1926) was calculated between the terrain attributes obtained from the RGB images and the ATI and the soil characteristics retrieved from ground 3 | RE SU LT S measurements. The average value of the terrain attributes on each plot was obtained in QGIS v.2.18.12 by averaging the terrain attribute 3.1 | Unmanned aerial vehicle (UAV) imagery values of the corresponding polygons, excluding vegetation pixels. geometric accuracy The correlation ratio was calculated as in Equation 9. It can assume values in the interval (0,1), and gives an indication of how much the The geometric accuracy of the products generated from the UAV data variance (in our case the soil property) is explained by the terrain imagery was assessed based on the RMSE, mean and standard factors; a high correlation ratio (close to 1) means that the factor deviation of the GCPs used for model generation and the GVPs explains most of the data variance. The weighted variance of the used for the assessment of model accuracy. A summary of the category means divided by the variance of all the samples. In our case statistics is reported in Table 2. The total RMSE computed on the the values of each terrain attribute (category) were divided in three GVPs was lower than 3.2 cm for both RGB models, about 22 cm classes standardized using their maximum and minimum values for the multispectral images and lower than 25 cm for the TIR (i.e. low, medium, high): images. The mean error over the whole study area was lower than 1 cm in the RGB models. For the multispectral and the TIR images, P  2 P  2 the mean error was always lower than 2 cm, in many cases lower σ 2y x nx y x  y nx yxi  y η2 ¼ , where σ 2y ¼ P and σ 2y ¼ x,i ð9Þ than 1 cm. Focusing on the RGB survey used to derive the DSMs, σ 2y x nx n DTMs and terrain attributes, the mean vertical error and standard deviation were very low. This suggests that the where yxi is the single soil property observation, x indicates the terrain generated models and terrain attribute maps were not affected by attribute category, and i indicates an observation. If nx is the number systematic error.  ET AL. BLANCO-SACRISTAN 2475 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) GCP GCP GCP GVP GVP GVP Area Imagery RMSE XY (cm) RMSE Z (cm) Total RMSE (cm) RMSE XY (cm) RMSE Z (cm) Total RMSE (cm) A RGB 2.17 4.2 4.7 2.2 2.3 3.2 B 1.58 1.3 2.1 1.9 1.3 2.4 A Multispectral 18.8 21.8 22.4 19.4 13.7 21.8 B 11.1 17.1 18.2 13.7 21.4 22.3 A TIR 21.2 23.4 24.7 23.1 22.4 24.5 B 16.6 17.3 19.8 16.8 15.9 19.9 MEjSTDEV X (cm) MEjSTDEV Y (cm) MEjSTDEV Z (cm) MEjSTDEV X (cm) MEjSTDEV Y (cm) MEjSTDEV Z (cm) A RGB 0.09j 0.49 0.02j 0.91 0.25j 1.40 0.09j 0.51 0.01j 0.49 0.17j 0.95 B 0.08j 0.16 0.01j 0.13 0.03j 0.41  0.1j 0.41 0.01j 0.13 0.04j 0.47 A Multispectral 0.83j 4.28 0.15j 4.72 1.2j 6.7 1.27j 8.45 0.25j 2.89 1.8j 2.27 B 0.6j 1.1 0.2j 1.7 0.3j 3.5 0.9j 2.6 0.3j 2.5 0.4j 6.1 A TIR 0.93j 4.82 0.16j 5.32 1.32j 5.62 1.01j 5.25 0.15j 5.09 1.31j 5.57 B 0.89j 1.7 0.26j 1.77 0.3j 2.03 0.91j 1.72 0.24j 1.63 0.33j 2.04 F I G U R E 5 Mean reflectance and standard deviation spectra (a) and mean continuum removed (CR) and standard deviation spectra (b) of the vegetation, biocrust and soil classes used as training for the support vector machine (SVM) classification. DryVeg, dry vegetation; GreenVeg, green vegetation; BLM, bright lichens with moss; Fulg, Fulgensia spp. with moss [Color figure can be viewed at wileyonlinelibrary.com] 3.2 | Classification of the multispectral similar in the spectral shape, DryVeg and Moss presented orthomosaics and evaluation of the classification differences in their CRred and reflectance values. The red absorption accuracy feature was present in all classes but Soil, and allowed to better differentiate it from biocrust classes. These subtle differences Mean reflectances and standard deviations of the training classes allowed an improvement in the classification carried out by the used for the classification procedure are shown in Figure 5(a). SVM, which slightly confused these two classes, by using thresholds Classes dominated by lichens had higher reflectance, with BL the in the CRred. highest reflectance class among them. Classes dominated by mosses The SVM classification generated land cover maps of the two had the lowest reflectance among biocrusts, only DryVeg target areas in the study site (Figure 6). The SVM classifications reflectance was lower. The spectra CR shown in Figure 5(b) high- were improved using the CRred thresholds on the Soil and DryVeg lights the absorption due to chlorophyll at 660 nm. Moss was the classes. Final accuracies are reported in Table S2, OAs higher than biocrust class with the highest absorption at 660 nm due to 80% in most cases were achieved in both areas, with κ equal to 0.8 chlorophyll (i.e., lowest CRred value), while GreenVeg was the class and 0.93 in areas A and B, respectively. Only BLM in Area A had an showing the highest absorption of all the evaluated classes. While OA below 80%. As shown in the confusion matrices (Table S2), BL 2476  ET AL. BLANCO-SACRISTAN and GreenVeg were the most accurately classified covers in this et al., 2018), confirming that the moisture status we captured with the area, while there was 10% of Moss misclassified as BLM. There was drone flight was representative of the winter rainy season. ATI values some confusion between Moss and DryVeg, with 10% of Moss correlated well with the soil moisture (r2 = 0.83; Figure 7), indicating misclassified as DryVeg in Area A, and between BLM and Soil, with that ATI is a good estimator of soil moisture, even when soils are cov- 14% of BLM classified as Soil. In both areas, the most dominant ered by biocrusts. ATI was then used as a proxy of soil moisture for class was DryVeg, while the least ones were the classes with lichens vegetation-free areas in the subsequent analysis. and mosses mixed (i.e., BLM and Fulg in Area A and Area B, respectively; Table 3). 3.4 | Analysis of soil properties 3.3 | Thermal data and soil moisture Soils in the study area are Gypsiric Leptosols and are characterized by the very high content of SalSol, reaching values higher than 90% in The soil moisture values measured in the 23 field plots were similar to some of the plots, and the low presence of OrgCarb (Table 4). Soil that of the winter long-term soil moisture of this area (Lafuente samples showed similar low nutrient contents in all plots, reaching F I G U R E 6 Support vector machine classifications (top) of the multispectral images (bottom), represented as false-colour composite (near infrared [NIR], green and red bands), in the six principal surface components. Bare soil, patches of bare soil; Bright lichens, patches of Squamarina spp., Diploschistes spp., Buellia spp.; Bright lichens with moss, patches of Squamarina spp., Diploschistes spp., Buellia spp. and moss; Moss, Tortula revolvens, and/or Syntrichia ruralis and/or Pleurochaete squarrosa; Fulg, Fulgensia spp. and moss; Dry vegetation, dry individuals of Helianthemum squamatum, Gypsophila struthium and Macrochloa tenacissima; Green vegetation, green tussocks of Macrochloa tenacissima. Left: study Area A; right: study Area B [Color figure can be viewed at wileyonlinelibrary.com] TABLE 3 Fractional cover of each class in the study areas as resulted from the classification of the multispectral images Surface cover (%) Soil DryVeg GreenVeg BL BLM Fulg Moss Area A 4.01 59.88 10.23 1.32 10.21 — 14.32 Area B 8.94 43.11 8.03 0.86 — 14.84 24.19 Note: — indicates class not present. Soil, bare soil; DryVeg, dry vegetation; GreenVeg, green vegetation; BL, bright lichens; BLM, bright lichens and moss; Fulg, Fulgensia spp. and moss; Moss, moss.  ET AL. BLANCO-SACRISTAN 2477 values of Ptot of 0% and values of Ntot close to 0% in some of them. Among the terrain attributes investigated, Northernness, ATI and Moisture presented high variability among the plots, with mean values PSIR are the ones that explained most of the variance of the soil prop- around 14%. erties (Table 5). Soils more exposed to north and characterized by higher ATI show a lower content of SolSal and higher content of OrgCarb, Ptot and Ntot. Soils characterized by higher values of PSIR 3.5 | Statistical analysis show an opposite trend, with higher content of SolSal and lower content of OrgCarb, Ptot and Ntot. The variance of soil properties is 3.5.1 | Relationship between terrain attributes and similarly explained by Northernness, ATI and PSIR with values of η2 soils at plot level close to 0.4 when related to SalSol, OrgCarb, Ptot and Ntot, with slightly higher values of η2 driven by PSIR. 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 3.5.2 | Relationships between topography, higher in Zone 3, reaching values of 35% of slope gradient. As biocrusts and vegetation distribution at landscape level expected, LSF presented similar patterns as slope. The TWI presented values around 3.5 with values not statistically different in the four The best fitted RDA model excluded Northernness and Slope variables, zones. PSIR showed an opposite pattern compared to Northernness, which showed multicollinearity, with VIF values around 20. The RDA with higher values in Zones 1 and 3, which are mainly exposed to results indicated that terrain variables significantly explained variation in southwest and south, respectively, indicating higher evapotranspira- the biocrust distribution in our study area (39.8% of the total inertia; tion. PSIR presented lower values in Zone 2 mainly exposed to the F = 92.5; p < 0.001). The first two RDA-axes accounted for 34.0% of east and minimum values in Zone 4, mainly exposed to the north. the variation for biocrust data (Figure 9). The most significant variables Similar to Northernness, ATI, that is, soil moisture, presented the were: ATI (F = 115.2; p < 0.001), Elevation (F = 100.8; p < 0.001), highest values in Zone 4 and lower values in Zone 1, evidencing that Northernness (F = 32.0; p > 0.001) and Veg (F = 29.7; p < 0.001) the soil moisture is related to the aspect. (Table S3). For those variables, differences were found between areas dominated by lichens and mosses. Lichen-dominated biocrusts appeared on the left side of the first axis, which was characterized by low values of vegetation presence and ATI (i.e., low moisture) and high values of LSF, which indicates high potential soil erosion, from a topo- graphic point of view and high values of PSIR and Elevation. As evidenced in Table 5, low ATI and high PSIR values are related to high SolSal content and low OrgCarb, Ptot and Ntot. Among lichen-dominated biocrusts, Fulg class development was mostly driven by the gradient of elevation, BL was controlled by high LSF and BLM was driven by low ATI and low elevation. Conversely, the distribution of mosses typical of shaded areas (i.e., Pleurochaete squarrosa and Syntrichia ruralis) appeared on the left side of the first axis, where vegetation presence and soil moisture (high ATI values) were higher and soil was richer in OrgCarb, Ptot and Ntot. Nevertheless, Tortula revolvens develops in soils with simi- lar characteristics to those where lichens dominate, characterized by low elevation and low ATI and vegetation presence. 4 | DI SCU SSION F I G U R E 7 Linear regression between apparent thermal inertia (ATI) and soil moisture in the 23 field plots used in this study. Shaded areas represent 95% symmetrical confidence interval [Color figure Previous studies have evaluated the effect of the terrain on biocrust can be viewed at wileyonlinelibrary.com] spatial distribution (Durham et al., 2018; Rodríguez-Caballero TABLE 4 Soil attributes retrieved from the 23 2  2 m2 plots sampled in the field Soil attribute Mean Maximum Minimum Standard deviation Moisture (%) 14.3 22.6 4.44 5.3 SolSal (%) 81 96 19 15.02 OrgCarb (kg/m2) 0.92 3.96 0.14 0.74 Ptot (%) 0.6 7.26 0 1.44 Ntot (%) 0.08 0.37 0.01 0.07 Note: SolSal, soluble salts; OrgCarb, organic carbon; Ptot, total phosphorus; Ntot, total nitrogen. 2478  ET AL. BLANCO-SACRISTAN F I G U R E 8 Box plots of the terrain attributes for the different aspects of the plots deployed in the field. The different lowercase letters indicate statistically significant differences between aspects (post hoc Dunn’s test following a Kruskal–Wallis test). 1, 2, 3 and 4 indicate the zones of the study area, which are mainly exposed to west, east, south and north, respectively. LSF, length slope factor; TWI, topographic wetness index; PSIR, potential solar incoming radiation; ATI, apparent thermal inertia [Color figure can be viewed at wileyonlinelibrary.com] TABLE 5 Descriptive statistics (i.e., mean and standard deviation) of soil properties grouped in terciles of terrain attributes L M H η2 North SolSal 86.50  4.74 68.60  4.94 62.50  29.21 0.43 OrgCarb 1.24  0.57 1.39  0.62 3.27  2.49 0.35 Ptot 112.72  30.51 167.01  65.05 219.50  151.14 0.31 Ntot 0.06  0.03 0.07  0.05 0.17  0.13 0.35 ATI SolSal 86.09  3.80 83.30  12.02 56.00  32.04 0.42 OrgCarb 1.18  0.49 1.19  0.53 3.30  3.04 0.32 Ptot 112.81  28.91 119.00  48.69 251.66  167.50 0.37 Ntot 0.06  0.02 0.06  0.03 0.18  0.15 0.35 PSIR SolSal 56.00  30.04 78.25  12.31 86.31  4.67 0.45 OrgCarb 3.30  3.04 1.28  0.38 1.25  0.61 0.35 Ptot 251.66  167.50 125.00  68.01 116.43  27.38 0.37 Ntot 0.18  0.15 0.06  0.03 0.06  0.03 0.36 2 Note: North, northernness; ATI, apparent thermal inertia, PSIR, potential solar incoming radiation. L, low; M, medium; H, high. Correlation ratio (η ) among soil variables and topography variables is also reported. Only variables with η2 higher than 0.3 are represented. SolSal, soluble salts (%); OrgCarb, organic carbon (km/m2); Ptot, total phosphorus (%); Ntot, total nitrogen (%) et al., 2019). We have extended this to provided a reproducible frame- between typical semi-arid environmental components (i.e., vascular work based only on drone remotely-sensed information; integrating vegetation, biocrust with different composition and soils) reinforces RGB, multispectral and TIR imagery to evaluate how terrain attributes the ability of the SVM algorithm to distinguish between similar drive biocrust distribution. In this study, the fine spatial resolution of spectral classes (Plaza et al., 2009) and of the CR algorithm’s ability to UAV-based multispectral imagery allowed us to produce very accurate improve identification of biocrusts (e.g., Blanco-Sacristán et al., 2019; maps of dryland biocrusts, with a higher degree of detail compared to Panigada et al., 2019; Román et al., 2019). Although differences in the the ones produced in previous studies using airborne and satellite absorption feature caused by chlorophyll at 660 nm between bio- imagery (e.g., Panigada et al., 2019; Rodríguez-Caballero et al., 2014; crusts and bare soils are subtle (Weber & Hill, 2016), we showed that Weber et al., 2008). The high accuracy reached in discriminating they can be detected with multispectral data. CRred allowed us to  ET AL. BLANCO-SACRISTAN 2479 F I G U R E 9 Redundancy analysis (RDA) of the biocrust and vegetation cover with respect to terrain attributes, soil and vegetation cover in the four zones of the study area. The red arrows indicate the increase along terrain attribute gradients, while the black arrows indicate the response variables (biocrust classes) and their directions represent the explanatory relationships with the axis. LSF, length slope factor; ATI, apparent thermal inertia; TWI, topographic wetness index; Veg, presence of vascular vegetation; Soil, presence of bare soil; PSIR, potential solar incoming radiation; BL, bright lichens; BLM, bright lichens with moss; Fulg, Fulgensia spp. with moss; TorRev, Tortula revolvens; Moss, Pleurochaete squarrosa and Syntrichia ruralis [Color figure can be viewed at wileyonlinelibrary.com] improve the classification results by exploiting the differences impact on soil properties. The mean error of these DTMs was lower between the biocrust classes. than 1 cm, which is comparable to previous studies applying SfM It should be noted that UAV data collection took place during the techniques for topography reconstruction (Bakker & Lane, 2017; wet season, when the water content in biocrusts is high and pigments Smith & Vericat, 2015). These results, together with the rough and are metabolically active (Weber & Hill, 2016). This enhances the dif- steep topography of the study areas, suggest that these models were ferences in the spectral properties of biocrusts and improves their not affected by systematic error. We found that PSIR, ATI and identification (Blanco-Sacristán et al., 2019). While the classifications Northernness were the terrain attributes mostly related to soil proper- from the multispectral images allowed mapping of the different envi- ties. These attributes explained in similar magnitudes the distribution ronmental components of the study area, the values of ATI derived of SolSal, OrgCarb, Ptot and Ntot. Increased shadows and higher soil from TIR imagery correlated very highly with the soil moisture mea- moisture, indicated by low PSIR and high ATI, are related to a higher sured in the field. We demonstrated that the correlation between ATI content of soil nutrients, as shown in plots located in Zones 2 and and soil moisture is also maintained when biocrust cover is present. 4, predominantly east- and north-facing plots, respectively. Increased To our knowledge, these are the first maps of ATI derived from UAV shadows and soil moisture control the survival and activity of microor- in drylands with presence of biocrust cover. Previous studies have ganisms (Borken & Matzner, 2009; Drenovsky et al., 2004) and lead to used hyperspectral thermal sensors to discriminate different types of better nutrient cycling and higher activity of microbial communities biocrusts from bare sand based on emissivity values and to character- (e.g., Xue et al., 2018). This increased moisture favours not only the ize biocrusts’ maturity (Rozenstein & Karnieli, 2015). In this study the preservation of OrgCarb and its association with other mineral com- ATI maps derived from TIR imagery were used to estimate how ponents (Plaza et al., 2013), but also the presence of vegetation and the soil moisture affects biocrust type and distribution, but ATI maps moss-dominated biocrusts on the most humid slopes. Conversely, could also be used to better understand water distribution in drylands, plots from Zones 1 and 3 presented the highest PSIR, which deter- especially when there is controversy about the final role of biocrust mines the lower content of nutrients (lower microbial activity-lower covers on soil water balance due to their influence on several hydro- nutrients) and increased content of SolSal, which is observed in other logical processes and soil properties with antagonist effects (Chamizo areas with high evapotranspiration rates (Rodríguez-Caballero et al., 2016). We found that topographic characteristics can partially et al., 2019). explain the soil moisture gradient estimated by ATI (ATI and PSIR are The variance in soil properties not explained by the terrain inversely correlated with r = 0.46). This moisture gradient promotes 2 attributes is probably explained by the presence of vegetation and a diversification of ecological niches that also explains the distribution biocrust cover. Vegetation patches enhance the accumulation of pattern of biocrust communities. However, long-term studies are still water and nutrient capture, which increase biological activity under necessary to better understand the dynamics of soil moisture and the and close to their cover (e.g., Okin et al., 2015). In addition, vegetation feedback effects on biocrust cover, and to try to decouple the increases shadows which promotes microbial activity (e.g., Huang biocrust effect on water from the water effect on biocrusts, which are et al., 2015; Xue et al., 2018) and thus soil nutrients content. strongly interrelated. ATI could be a useful variable to expand our Macrochloa tenacissima in particular, the dominant species in the study knowledge in this direction. area, modifies stocks of soil carbon and organic matter (e.g. Gauquelin The UAV RGB images were used to produce fine DTMs from et al., 1996; Kaouthar & Chaieb, 2009; Maestre et al., 2001). While its which terrain attributes were derived and utilized to evaluate their effect on nirogen (N) content is less clear, some studies have found 2480  ET AL. BLANCO-SACRISTAN negative effects on this parameter (Armas & Pugnaire, 2011) while dominant species of BL, was also found in areas with high LSF. This others have described positive effects (Castillo-Monroy et al., 2010) species is very versatile and can adapt its physiology depending on when compared to non-vegetated surfaces. where it grows (Pintado et al., 2005). We observed that Diploschistes Together with the effect of vegetation on soil properties, diacapsis can develop directly on outcrop rocks with very few biocrusts have been found to take up significant amounts of atmo- millimetres of soil beneath them in areas where gypsum substrates spheric C and N by photosynthesis and N fixation (Elbert et al., 2012). with fine soil texture favour stability, but where water availability is They are therefore an important pool and source of organic inputs scarse. In these areas with outcrop rock, high concentrations of gyp- into the soil in drylands (Castillo-Monroy et al., 2011; Chamizo sum, soluble salts and higher solar radiation (i.e., high values of PSIR), et al., 2012b; Concostrina-Zubiri et al., 2013). Furthermore, biocrusts Diploschistes diacapsis is more adapted to establish than other lichens, strengthen soil structure by interacting with mineral particles and representing here the early stage of biocrust development. forming aggregates (Belnap & Lange, 2013; Eldridge & Leys, 2003), We found instead mixed patches of lichens and moss (BLM) thus protecting soil from C loss. Biocrusts promote microbial (thought to be of the last stages of development) in the left axis of the community growth where they appear and the increased shadows RDA plot, where water availability is still scarces but the terrain created by vegetation increase this growth indirectly (Huang presentes lower elevation and LSF and less concentrations of gypsum. et al., 2015). However, the role of biocrusts on soil nutrients content Moss-dominated biocrusts become dominant in areas with higher soil depends on biocrust composition and patch-size distribution (e.g., moisture and north-facing areas, where soil and micro-environmental Bowker et al., 2013; Delgado-Baquerizo et al., 2015; Sedia & conditions are less selective. In the Tabernas desert (southern Spain, Ehrenfeld, 2006). For example, mosses have greater photosynthetic Almería; Rodríguez-Caballero et al., 2019), a nearby site with 220 mm capacity compared to other types of biocrust (Weber & Hill, 2016), of mean annual rainfall, similar lichen-dominated biocrusts appear in and thus can incorporate higher levels of C. In addition, positive rela- north-faced slopes (low evapotranspiration, low PSIR), while in our tionships between N and the presence of mosses have already been study area they appear in south-facing slopes (high evapotranspira- observed (Delgado-Baquerizo et al., 2016). This might partially explain tion, high PSIR). Biocrust distribution in drylands might be more con- the increased soil nutrients content in the north- and east-facing plots strained by the effect of terrain attributes on local hydric availability, in the study area, where mosses appear to dominate plant rather than rainfall water inputs. For this reason, it is difficult to make interspaces. generalizations in the development of biocrusts’ communities and Variations in biocrust distribution and composition in dryland developing monitoring methodologies that allow up-scaling local biomes are primarily controlled by climatic differences (Bowker et al., relationships of biocrusts with the terrain is key for dryland ecology. 2016). However, vegetation and terrain properties generate variations from these climatic controls and modify the local distribution of nutri- ents in dryland landscapes (e.g., Manzoni et al., 2006; 5 | CONC LU SIONS Puigdefábregas, 2005; Puigdefábregas & Sánchez, 1996). Thus, the different sensitivity of biocrusts to the distribution of resources condi- The synergistic information obtained from RGB, multispectral and tions the composition of these ecosystems (Bowker et al., 2016; Dur- thermal UAV imagery has allowed us to better understand the ham et al., 2018; Rodríguez-Caballero et al., 2019). The RDA complex factors influencing the distribution of dryland ecosystem conducted using data from 1147 points showed that Fulgensia spp. components. We produced accurate, high-resolution maps that can presented a strong positive relationship with unstable zones (high help to up-scale the local effect of biocrust components on ecosystem values of elevation and LSF). Fulgensia spp. has already been observed functioning by improving the accuracy of erosion and infiltration in unstable terrains several times, highlighting the pioneering behav- modelling. These aspects of dryland ecology are key to understanding  n et al., 2020; Miralles et al., 2020; iour of this genus (Canto the distinctive role that biocrust-dominated surfaces have on water Rodríguez-Caballero et al., 2013). Mosses, typically found in shaded runoff and erosion depending on the predominant type of biocrust areas (i.e. Pleurochaete squarrosa and Syntrichia ruralis), dominated (e.g., Rodríguez-Caballero et al., 2015; Wang et al., 2017). Producing areas with low PSIR, high ATI, and a high presence of vascular vegeta- high-resolution maps will help to monitor these communities in space tion, where soils were richer in organic carbon and nutrients. A posi- and time, a key task to understand the compositional changes they tive plant-biocrust relationship is common for bryophytes (Zhou are already experiencing due to climate change (Escolar et al., 2012; et al., 2020) and close correspondence with Macrochloa tenacissima  n de Guevara et al., 2018). Ladro presence has already been observed (Martínez-Sánchez et al., 1994). Since field-based data collection has many drawbacks, mainly The RDA showed a lower explanatory power of the terrain related to time and costs (Palmer et al., 2002), exploiting UAV-based attributes regarding Tortula revolves, a moss that develops in arid methodologies will help to build standardized procedures in dryland environments, and BLM class dominated by the association of monitoring programmes, while providing very detailed information of Diploschistes diacapsis and Squamarina lentigera (BLM). They appeared these environments. These are needed because drylands may be in more stable zones (i.e., low values of elevation and LSF), as already characterized by different structural and functional organization and n observed in geographic areas with similar characteristics (e.g., Ladro present a wide spectrum of compositions, depending on local climate de Guevara et al., 2018) that have high solar radiation and soils rich in and terrain properties. Our results have shown some discrepancies soluble salts. In particular, Squamarina lentigera can be physiologically compared to the traditional biocrust development models, in which adapted to light-exposed environments, and requires high more advanced stages (i.e., mixed patches of lichens and moss) are temperatures for optimal photosynthesis while being well hydrated assumed to appear only in stable soils. Thus, we suggest integrating (Lange et al., 1997). Nevertheless, Diploschistes diacapsis, the more traditional approaches with new methodologies such as the one  ET AL. BLANCO-SACRISTAN 2481 developed in this study, which can provide the very fine spatial resolu- Bowker, M.A., Maestre, F.T. & Mau, R.L. (2013) Diversity and patch-size tion maps of dryland composition and terrain attributes needed for a distributions of biological soil crusts regulate dryland ecosystem mul- tifunctionality. Ecosystems, 16(6), 923–933. Available from: https:// detailed description and monitoring of these complex and threatened doi.org/10.1007/s10021-013-9644-5 ecosystems. Brodu, N. & Lague, D. (2012) 3D terrestrial lidar data classification of com- plex natural scenes using a multi-scale dimensionality criterion: ACKNOWLEDGEMEN TS Applications in geomorphology. ISPRS Journal of Photogrammetry and Remote Sensing, 68(1), 121–134. Available from: https://doi.org/10. This research has received funding from the European Union’s 1016/j.isprsjprs.2012.01.006 Horizon 2020 research and innovation programme under the Marie Campbell, G.S. & Norman, J.M. (1988) Introduction to Environmental Bio- Sklodowska-Curie grant agreement No 721995. It also received physics, second edition. New York: Springer.  & Maestre, F.T. (2018) funding by MIUR – Dipartimenti di Eccellenza 2018–2022. The work Cano-Díaz, C., Mateo, P., Muñoz-Martín, M.A. Diversity of biocrust-forming cyanobacteria in a semiarid gypsiferous of M.L.d.G was supported by the European Research Council grant site from central Spain. Journal of Arid Environments, 151, 83–89. agreement no. 647038 (BIODESERT). The authors acknoledge Fabio Available from: https://doi.org/10.1016/j.jaridenv.2017.11.008 Moia and Davide Abu El Khair (University of Milano-Bicocca) for soil Canto n, Y., Chamizo, S., Rodriguez-Caballero, E., Lázaro, R., Roncero- laboratory measurements. Ramos, B., Román, J.R. & Solé-Benet, A. (2020) Water regulation in cyanobacterial biocrusts from drylands: negative impacts of anthro- pogenic disturbance. Water, 12(3), 720. Available from: https://doi. CONF LICT OF IN TE RE ST org/10.3390/w12030720 The authors declare there is no conflict of interest. Castillo-Monroy, A.P., Maestre, F.T., Delgado-Baquerizo, M. & Gallardo, A. (2010) Biological soil crusts modulate nitrogen availability in semi- DATA AVAILABILITY STATEMENT arid ecosystems: Insights from a Mediterranean grassland. Plant and Soil, 333(1–2), 21–34. Available from: https://doi.org/10.1007/ The data that support the findings of this study are available from the s11104-009-0276-7 corresponding author upon reasonable request. Castillo-Monroy, A.P., Maestre, F.T., Rey, A., Soliveres, S. & García- Palacios, P. (2011) Biological soil crust microsites are the main ORCID contributor to soil respiration in a semiarid ecosystem. Ecosystems, Cinzia Panigada https://orcid.org/0000-0002-3745-8037 14(5), 835–847. Chamizo, S., Canto  n, Y., Lázaro, R., Solé-Benet, A. & Domingo, F. (2012a) Crust composition and disturbance drive infiltration through biologi- RE FE R ENC E S cal soil crusts in semiarid ecosystems. Ecosystems, 15(1), 148–161. Adessi, A., de Carvalho, R.C., De Philippis, R., Branquinho, C. & da Silva, J. Chamizo, S., Canto n, Y., Miralles, I. & Domingo, F. (2012b) Biological soil M. (2018) Microbial extracellular polymeric substances improve crust development affects physicochemical characteristics of soil water retention in dryland biological soil crusts. Soil Biology and surface in semiarid ecosystems. Soil Biology and Biochemistry, 49, Biochemistry, 116, 67–69. Available from: https://doi.org/10.1016/j. 96–105. Available from: https://doi.org/10.1016/j.soilbio.2012. soilbio.2017.10.002 02.017 Aguiar, M.R. & Sala, O.E. (1999) Patch structure, dynamics and implications Chamizo, S., Canto  n, Y., Rodríguez-Caballero, E. & Domingo, F. (2016) for the functioning of arid ecosystems. Trends in Ecology & Evolution, Biocrusts positively affect the soil water balance in semiarid ecosys- 14(7), 273–277. Available from: https://doi.org/10.1016/S0169- tems. Ecohydrology, 9(7), 1208–1221. Available from: https://doi. 5347(99)01612-2 org/10.1002/eco.1719 Armas, C. & Pugnaire, F.I. (2011) Belowground zone of influence in a tus- Chamizo, S., Rodríguez-Caballero, E., Román, J.R. & Canto n, Y. (2017) sock grass species. Acta Oecologica, 37(3), 284–289. Available from: Effects of biocrust on soil erosion and organic carbon losses under https://doi.org/10.1016/j.actao.2011.02.013 natural rainfall. Catena, 148, 117–125. Available from: https://doi. Assmann, J.J., Kerby, J.T., Cunliffe, A.M. & Myers-Smith, I.H. (2018) org/10.1016/j.catena.2016.06.017 Vegetation monitoring using multispectral sensors—best practices Chauhan, S., Darvishzadeh, R., Lu, Y., Boschetti, M. & Nelson, A. (2020) and lessons learned from high latitudes. Journal of Unmanned Vehicle Understanding wheat lodging using multi-temporal Sentinel-1 and Systems, 7(1), 54–75. Sentinel-2 data. Remote Sensing of Environment, 243(March), 1–14. Bakker, M. & Lane, S.N. (2017) Archival photogrammetric analysis of river- Available from: https://doi.org/10.1016/j.rse.2020.111804 floodplain systems using Structure from Motion (SfM) methods. Clark, R.N. & Roush, T.L. (1984) Reflectance spectroscopy. Quantitative Earth Surface Processes and Landforms, 42(8), 1274–1286. Available analysis techniques for remote sensing applications. Journal of from: https://doi.org/10.1002/esp.4085 Geophysical Research, 89(B7), 6329–6340. Available from: https:// Belnap, J. & Lange, O.L. (2013) Biological Soil Crusts: Structure, Function, doi.org/10.1029/JB089iB07p06329 and Management. Berlin: Springer. Cohen, J. (1960) A coefficient of agreement for nominal scales. Educational Beven, K.J. & Kirkby, M.J. (1979) A physically based, variable and Psychological Measurement, 20(1), 37–46. Available from: contributing area model of basin hydrology. Hydrological Sciences https://doi.org/10.1177/001316446002000104 Bulletin, 24(1), 43–69. Available from: https://doi.org/10.1080/ Colombo, C. & Miano, T. (2015) Metodi di analisi chimica del suolo. Società 02626667909491834 italiana della scienza del suolo; Associazione Italiana dei Laboratori Blanco-Sacristán, J., Panigada, C., Tagliabue, G., Gentili, R., Colombo, R., Pubblici di Agrochimica. Modugno (Bari): Pubblicità & Stampa, Ladro n de Guevara, M. et al. (2019) Spectral diversity successfully pp. 1–470. estimates the α-diversity of biocrust-forming lichens. Remote Sensing, Concostrina-Zubiri, L., Huber-Sannwald, E., Martínez, I., Flores Flores, J. 11(24), 2942. Available from: https://doi.org/10.3390/rs11242942 L. & Escudero, A. (2013) Biological soil crusts greatly contribute to Borken, W. & Matzner, E. (2009) Reappraisal of drying and wetting effects small-scale soil heterogeneity along a grazing gradient. Soil Biology on C and N mineralization and fluxes in soils. Global Change Biology, and Biochemistry, 64, 28–36. Available from: https://doi.org/10. 15(4), 808–824. Available from: https://doi.org/10.1111/j.1365- 1016/j.soilbio.2013.03.029 2486.2008.01681.x Cortes, C. & Vapnik, V. (1995) Support-vector networks. Machine Learning, Bowker, M.A., Belnap, J., Büdel, B., Sannier, C., Pietrasiak, N., Eldridge, D. 20(3), 273–297. Available from: https://doi.org/10.1007/ J. & Rivera-Aguilar, V. (2016) Controls on distribution patterns of bf00994018 biological soil crusts at micro to global scales. In: Biological Soil Crusts: Cunliffe, A.M., Brazier, R.E. & Anderson, K. (2016) Ultra-fine grain An Organizing Principle in Drylands. Cham: Springer, pp. 173–197. landscape-scale quantification of dryland vegetation structure with 2482  ET AL. BLANCO-SACRISTAN drone-acquired structure-from-motion photogrammetry. Remote Hegarty, C.J. & Chatre, E. (2008) Evolution of the global navigation satel- Sensing of Environment, 183, 129–143. Available from: https://doi. lite system (GNSS). Proceedings of the IEEE, 96(12), 1902–1917. org/10.1016/j.rse.2016.05.019 Available from: https://doi.org/10.1109/JPROC.2008.2006090 Curran, P.J. (1988) The semivariogram in remote sensing: An introduction. Hofmann-Wellenhof, B., Lichtenegger, H. & Wasle, E. (2007) GNSS–global Remote Sensing of Environment, 24(3), 493–507. Available from: navigation satellite systems: GPS, GLONASS, Galileo, and more. Berlin: https://doi.org/10.1016/0034-4257(88)90021-1 Springer Science & Business Media. Delgado-Baquerizo, M., Gallardo, A., Covelo, F., Prado-Comesa.a, A., Huang, Y.M., Liu, D. & An, S.S. (2015) Effects of slope aspect on soil nitro- Ochoa, V. & Maestre, F.T. (2015) Differences in thallus chemistry are gen and microbial properties in the Chinese Loess region. Catena, related to species-specific effects of biocrust-forming lichens on soil 125, 135–145. Available from: https://doi.org/10.1016/j.catena. nutrients and microbial communities. Functional Ecology, 29(8), 2014.09.010 1087–1098. Available from: https://doi.org/10.1111/1365-2435.12403 IUSS Working Group WRB. (2006) World Reference Base for Soil Resources Delgado-Baquerizo, M., Maestre, F.T., Eldridge, D.J., Bowker, M.A., 2006. World Soil Resources Reports No. 103. Rome: FAO. Ochoa, V., Gozalo, B. et al. (2016) Biocrust-forming mosses mitigate Jung, P., Baumann, K., Lehnert, L.W., Samolov, E., Achilles, S., the negative impacts of increasing aridity on ecosystem Schermer, M. et al. (2020) Desert breath—how fog promotes a novel multifunctionality in drylands. New Phytologist, 209(4), 1540–1552. type of soil biocenosis, forming the coastal Atacama Desert’s living Available from: https://doi.org/10.1111/nph.13688 skin. Geobiology, 18(1), 113–124. Available from: https://doi.org/10. Drenovsky, R., Vo, D., Graham, K. & Scow, K. (2004) Soil water content 1111/gbi.12368 and organic carbon availability are major determinants of soil micro- Kaouthar, J. & Chaieb, M. (2009) The effect of Stipa tenacissima tussocks bial community composition. Microbial Ecology, 48(3), 424–430. on some soil surface properties under arid bioclimate in the southern Available from: https://doi.org/10.1007/s00248-003-1063-2 Tunisia. Acta Botanica Gallica, 156(2), 173–181. Available from: Durham, R.A., Doherty, K.D., Antoninka, A.J., Ramsey, P.W. & Bowker, M. https://doi.org/10.1080/12538078.2009.10516149 A. (2018) Insolation and disturbance history drive biocrust biodiver- Koutroulis, A.G. (2019) Dryland changes under different levels of global sity in western Montana rangelands. Plant and Soil, 430(1), 151–169. warming. Science of the Total Environment, 655, 482–511. Available Available from: https://doi.org/10.1007/s11104-018-3725-3 from: https://doi.org/10.1016/j.scitotenv.2018.11.215 Elbert, W., Weber, B., Burrows, S., Steinkamp, J., Büdel, B., Andreae, M. Krzeminska, D.M., Steele-Dunne, S.C., Bogaard, T.A., Rutten, M.M., O. & Pöschl, U. (2012) Contribution of cryptogamic covers to the Sailhac, P. & Geraud, Y. (2012) High-resolution temperature global cycles of carbon and nitrogen. Nature Geoscience, 5(7), observations to monitor soil thermal properties as a proxy for soil 459–462. Available from: https://doi.org/10.1038/ngeo1486 moisture condition in clay-shale landslide. Hydrological Processes, Eldridge, D.J. & Leys, J.F. (2003) Exploring some relationships between 26(14), 2143–2156. Available from: https://doi.org/10.1002/ biological soil crusts, soil aggregation and wind erosion. Journal of hyp.7980 Arid Environments, 53(4), 457–466. Available from: https://doi.org/  n de Guevara, M., Gozalo, B., Raggio, J., Lafuente, A., Prieto, M. & Ladro 10.1006/jare.2002.1068 Maestre, F.T. (2018) Warming reduces the cover, richness and even- Eldridge, D.J., Reed, S., Travers, S.K., Bowker, M.A., Maestre, F.T., Ding, J. ness of lichen-dominated biocrusts but promotes moss growth: et al. (2020) The pervasive and multifaceted influence of biocrusts on Insights from an 8 yr experiment. The New Phytologist, 220, water in the world’s drylands. Global Change Biology, 26(10), 811–823. 6003–6014. Available from: https://doi.org/10.1111/gcb.15232 Lafuente, A., Berdugo, M., Ladron de Guevara, M., Gozalo, B. & Escolar, C., Martínez, I., Bowker, M.A. & Maestre, F.T. (2012) Warming Maestre, F.T. (2018) Simulated climate change affects how biocrusts reduces the growth and diversity of biological soil crusts in a semi- modulate water gains and desiccation dynamics after rainfall events. arid environment: implications for ecosystem structure and Ecohydrology, 11(6), e1935. Available from: https://doi.org/10.1002/ functioning. Philosophical Transactions of the Royal Society of London. eco.1935 Series B: Biological Sciences, 367, 3087–3099. Lange, O.L., Belnap, J., Reichenberger, H. & Meyer, A. (1997) Photosynthe- Faist, A.M., Herrick, J.E., Belnap, J., Van Zee, J.W. & Barger, N.N. (2017) sis of green algal soil crust lichens from arid lands in southern Utah, Biological soil crust and disturbance controls on surface hydrology in USA: Role of water content on light and temperature responses of a semi-arid ecosystem. Ecosphere, 8(3), e01691. Available from: CO2 exchange. Flora, 192(1), 1–15. Available from: https://doi.org/ https://doi.org/10.1002/ecs2.1691 10.1016/S0367-2530(17)30749-1 Fawcett, D., Panigada, C., Tagliabue, G., Boschetti, M., Celesti, M., Lowe, D.G. (2004) Distinctive image features from scale-invariant Evdokimov, A. et al. (2020) Multi-scale evaluation of drone-based keypoints. International Journal of Computer Vision, 60(2), 91–110. multispectral surface reflectance and vegetation indices in opera- Available from: https://doi.org/10.1023/B:VISI.0000029664. tional conditions. Remote Sensing, 12(3), 514. 99615.94 Fu, P. & Rich, P.M. (1999) Design and implementation of the Solar Analyst: Lucieer, A., Turner, D., King, D.H. & Robinson, S.A. (2014) Using an an ArcView extension for modeling solar radiation at landscape unmanned aerial vehicle (UAV) to capture micro-topography of Ant- scales. In: Proceedings of the nineteenth annual ESRI user conference. arctic moss beds. International Journal of Applied Earth Observation Vol. 1. USA: San Diego. and Geoinformation, 27, 53–62. Available from: https://doi.org/10. Gao, L., Bowker, M.A., Sun, H. & Zhao, Y. (2020) Linkages between bio- 1016/j.jag.2013.05.011 crust development and water erosion and implications for erosion Maestre, F.T., Bautista, S., Cortina, J. & Bellot, J. (2001) Potential for using model implementation. Geoderma, 357, 113973. Available from: facilitation by grasses to establish shrubs on a semiarid degraded https://doi.org/10.1016/j.geoderma.2019.113973 steppe. Ecological Applications, 11(6), 1641–1655. Available from: Gauquelin, T., Jalut, G., Iglesias, M., Valle, F., Fromard, F. & Dedoubat, J.J. https://doi.org/10.1890/1051-0761(2001)011[1641:PFUFBG]2.0. (1996) Phytomass and carbon storage in the Stipa tenacissima steppes CO;2 of the Baza basin, Andalusia, Spain. Journal of Arid Environments, 34(3), Maestre, F.T., Escolar, C., Ladro n de Guevara, M., Quero, J.L., Lázaro, R., 277–286. Available from: https://doi.org/10.1006/jare.1996.0109 Delgado-Baquerizo, M. et al. (2013) Changes in biocrust cover drive Gindraux, S., Boesch, R. & Farinotti, D. (2017) Accuracy assessment of carbon cycle responses to climate change in drylands. Global Change digital surface models from unmanned aerial vehicles’ imagery on Biology, 19(12), 3835–3847. Available from: https://doi.org/10. glaciers. Remote Sensing, 9(2), 186. Available from: https://doi.org/ 1111/gcb.12306 10.3390/rs9020186 Maltese, A., Capodici, F., Ciraolo, G. & La Loggia, G. (2013) Mapping soil Havrilla, C.A., Villarreal, M.L., DiBiase, J.L., Duniway, M.C. & Barger, N.N. water content under sparse vegetation and changeable sky condi- (2020) Ultra-high-resolution mapping of biocrusts with unmanned tions: Comparison of two thermal inertia approaches. Journal of aerial systems. Remote Sensing in Ecology and Conservation, 6(4), Applied Remote Sensing, 7(1), 073548. Available from: https://doi. 441–456. Available from: https://doi.org/10.1002/rse2.180 org/10.1117/1.JRS.7.073548  ET AL. BLANCO-SACRISTAN 2483 Manzoni, S., Porportato, A. & D’Odorico, P. (2006) Modelling of carbon Plaza, A., Benediktsson, J.A., Boardman, J.W., Brazile, J., Bruzzone, L., and nitrogen cycling in arid and semiarid ecosystems. In: Dryland Camps-Valls, G. et al. (2009) Recent advances in techniques for Ecohydrology. Dordrecht: Springer, pp. 183–199. hyperspectral image processing. Remote Sensing of Environment, 113, Martínez-Carricondo, P., Agüera-Vega, F., Carvajal-Ramírez, F., Mesas- S110–S122. Available from: https://doi.org/10.1016/j.rse.2007. Carrascosa, F.J., García-Ferrer, A. & Pérez-Porras, F.J. (2018) Assess- 07.028 ment of UAV-photogrammetric mapping accuracy based on variation Plaza, C., Courtier-Murias, D., Fernández, J.M., Polo, A. & Simpson, A.J. of ground control points. International Journal of Applied Earth Obser- (2013) Physical, chemical, and biochemical mechanisms of soil vation and Geoinformation, 72, 1–10. Available from: https://doi.org/ organic matter stabilization under conservation tillage systems: A 10.1016/j.jag.2018.05.015 central role for microbes and microbial by-products in C sequestra- Martínez-Sánchez, J.J., Casares-Porcel, M., Guerra, J., Gutiérrez- tion. Soil Biology and Biochemistry, 57, 124–134. Available from: Carretero, L., Ros, R.M., Hernández-Bastida, J. & Cano, M.J. (1994) A https://doi.org/10.1016/j.soilbio.2012.07.026 special habitat for bryophytes and lichens in the arid zones of Spain. Price, J.C. (1985) On the analysis of thermal infrared imagery: The limited Lindbergia, 19(2/3), 116–121. utility of apparent thermal inertia. Remote Sensing of Environment, Milling, C.R., Rachlow, J.L., Olsoy, P.J., Chappell, M.A., Johnson, T.R., 18(1), 59–73. Available from: https://doi.org/10.1016/0034-4257 Forbey, J.S. et al. (2018) Habitat structure modifies microclimate: An (85)90038-0 approach for mapping fine-scale thermal refuge. Methods in Ecology Puigdefábregas, J. (2005) The role of vegetation patterns in structuring and Evolution, 9(6), 1648–1657. Available from: https://doi.org/10. runoff and sediment fluxes in drylands. Earth Surface Processes and 1111/2041-210X.13008 Landforms, 30(2), 133–147. Available from: https://doi.org/10.1002/ Minacapilli, M., Iovino, M. & Blanda, F. (2009) High resolution remote esti- esp.1181 mation of soil surface water content by a thermal inertia approach. Puigdefábregas, J. & Sánchez, G. (1996) Geomorphological implications of Journal of Hydrology, 379(3–4), 229–238. Available from: https://doi. vegetation patchiness on semi-arid slopes. In: Anderson, M.G. & org/10.1016/j.jhydrol.2009.09.055 Brooks, S.M. (Eds.) Advances on Hillslope Processes. Chichester: John Miralles, I., Lázaro, R., Sánchez-Maraño  n, M., Soriano, M. & Ortega, R. Wiley & Sons, pp. 1027–1060. (2020) Biocrust cover and successional stages influence soil bacterial Puigdefábregas, J., Sole, A., Gutierrez, L., Del Barrio, G. & Boer, M. (1999) composition and diversity in semiarid ecosystems. Science of the Scales and processes of water and sediment redistribution in dry- Total Environment, 709, 134654. Available from: https://doi.org/10. lands: Results from the Rambla Honda field site in southeast Spain. 1016/j.scitotenv.2019.134654 Earth-Science Reviews, 48(1–2), 39–70. Available from: https://doi. Mitra, D.S. & Majumdar, T.J. (2004) Thermal inertia mapping over the org/10.1016/S0012-8252(99)00046-X Brahmaputra basin, India using NOAA-AVHRR data and its possible Renard, K.G., Foster, G.R., Weesies, G.A., McCool, D.K. & Yoder, D.C. geological applications. International Journal of Remote Sensing, (1997) Predicting Soil Erosion by Water: A guide to conservation plan- 25(16), 3245–3260. ning with the revised universal soil loss equation (RUSLE). US Depart- Monteith, J.L. & Szeicz, G. (1962) Radiative temperature in the heat bal- ment of Agriculture, Agriculture Handbook No. 703. Washington, DC: ance of natural surfaces. Quarterly Journal of the Royal Meteorological United States Government Printing. Society, 88(378), 496–507. Available from: https://doi.org/10.1002/ Rodríguez-Caballero, E., Aguilar, M.A.,  Castilla, Y.C., Chamizo, S. & qj.49708837811 Aguilar, F.J. (2015) Swelling of biocrusts upon wetting induces Nearing, G.S., Moran, M.S., Scott, R.L. & Ponce-Campos, G. (2012) Cou- changes in surface micro-topography. Soil Biology and Biochemistry, pling diffusion and maximum entropy models to estimate thermal 82, 107–111. Available from: https://doi.org/10.1016/j.soilbio.2014. inertia. Remote Sensing of Environment, 119, 222–231. Available 12.010 from: https://doi.org/10.1016/j.rse.2011.12.012 Rodríguez-Caballero, E., Canto  n, Y., Chamizo, S., Lázaro, R. & Escudero, A. Negm, A., Capodici, F., Ciraolo, G., Maltese, A., Provenzano, G. & Rallo, G. (2013) Soil loss and runoff in semiarid ecosystems: A complex (2017) Assessing the performance of thermal inertia and Hydrus interaction between biological soil crusts, micro-topography, and models to estimate surface soil water content. Applied Sciences, hydrological drivers. Ecosystems, 16(4), 529–546. Available from: 7(10), 975. Available from: https://doi.org/10.3390/app7100975 https://doi.org/10.1007/s10021-012-9626-z Ogle, D.H., Wheeler, P. & Dinno, A. (2021) FSA: Fisheries Stock Analysis. Rodríguez-Caballero, E., Escribano, P. & Canto  n, Y. (2014) Advanced image R package version 0.8.32, https://github.com/droglenc/FSA processing methods as a tool to map and quantify different types of Okin, G.S., Heras, M.M.L., Saco, P.M., Throop, H.L., Vivoni, E.R., biological soil crust. ISPRS Journal of Photogrammetry and Remote Parsons, A.J. et al. (2015) Connectivity in dryland landscapes: Shifting Sensing, 90, 59–67. Available from: https://doi.org/10.1016/j. concepts of spatial interactions. Frontiers in Ecology and the Environ- isprsjprs.2014.02.002 ment, 13(1), 20–27. Available from: https://doi.org/10.1890/140163 Rodríguez-Caballero, E., Escribano, P., Olehowski, C., Chamizo, S., Hill, J., Oksanen, J., Blanchet, F.G., Kindt, R., Legendre, P., O’Hara, R.B., Canto  n, Y. & Weber, B. (2017a) Transferability of multi-and hyper- Simpson, G.L. et al. (2017) Vegan: Community Ecology Package. R spectral optical biocrust indices. ISPRS Journal of Photogrammetry and Package Version 2.4–5. Vienna: R Foundation for Statistical Comput- Remote Sensing, 126, 94–107. Available from: https://doi.org/10. ing Available online: https://CRAN.R-project.org/package=vegan 1016/j.isprsjprs.2017.02.007 (accessed. Rodríguez-Caballero, E., Paul, M., Tamm, A., Caesar, J., Büdel, B., Palmer, M.W., Earls, P., Hoagland, B.W., White, P.S. & Wohlgemuth, T. Escribano, P. et al. (2017b) Biomass assessment of microbial surface (2002) Quantitative tools for perfecting species lists. Environmetrics, communities by means of hyperspectral remote sensing data. Science 13(2), 121–137. Available from: https://doi.org/10.1002/env.516 of the Total Environment, 586, 1287–1297. Available from: https:// Panigada, C., Tagliabue, G., Zaady, E., Rozenstein, O., Garzonio, R., Di doi.org/10.1016/j.scitotenv.2017.02.141 Mauro, B. & Rossini, M. (2019) A new approach for biocrust and Rodríguez-Caballero, E., Román, J.R., Chamizo, S., Roncero Ramos, B. & vegetation monitoring in drylands using multi-temporal Sentinel-2 Canto  n, Y. (2019) Biocrust landscape-scale spatial distribution is images. Progress in Physical Geography: Earth and Environment, strongly controlled by terrain attributes: Topographic thresholds for 43(4), 496–520. Available from: https://doi.org/10.1177/ colonization in a semiarid badland system. Earth Surface Processes 0309133319841903 and Landforms, 44(14), 2771–2779. Available from: https://doi.org/ Pearson, E.S. (1926) Review of statistical methods for research workers 10.1002/esp.4706 (RA Fisher). Science Progress, 20, 733–734. Román, J.R., Rodríguez-Caballero, E., Rodríguez-Lozano, B., Roncero- Pintado, A., Sancho, L.G., Green, T.G.A., Blanquer, J.M. & Lázaro, R. (2005)  Ramos, B., Chamizo, S., Aguila-Carricondo, P. & Canto n, Y. (2019) Functional ecology of the biological soil crust in semiarid SE Spain: Spectral response analysis: An indirect and non-destructive sun and shade populations of Diploschistes diacapsis (Ach.) Lumbsch. methodology for the chlorophyll quantification of biocrusts. Remote The Lichenologist, 37(5), 425–432. Available from: https://doi.org/10. Sensing, 11(11), 1350. Available from: https://doi.org/10.3390/ 1017/s0024282905015021 rs11111350 2484  ET AL. BLANCO-SACRISTAN Rossi, G., Ferrarini, A., Dowgiallo, G., Carton, A., Gentili, R. & Tomaselli, M. Ullman, S. (1983) The interpretation of structure from motion. Proceedings (2014) Detecting complex relations among vegetation, soil and geo- of the Royal Society of London B, 203, 405–426. morphology. An in-depth method applied to a case study in the Van Doninck, J., Peters, J., De Baets, B., De Clercq, E.M., Ducheyne, E. & Apennines (Italy). Ecological Complexity, 17, 87–98. Available from: Verhoest, N.E.C. (2011) The potential of multitemporal Aqua and Terra https://doi.org/10.1016/j.ecocom.2013.11.002 MODIS apparent thermal inertia as a soil moisture indicator. Interna- Rozenstein, O. & Karnieli, A. (2015) Identification and characterization of tional Journal of Applied Earth Observation and Geoinformation, 13(6), biological soil crusts in a sand dune desert environment across Israel- 934–941. Available from: https://doi.org/10.1016/j.jag.2011.07.003 Egypt border using LWIR emittance spectroscopy. Journal of Arid Vapnik, B.V. (2005) Universal learning technology: Support vector Environments, 112(PA), 75–86. Available from: https://doi.org/10. machines. Journal of Advances in Technology, 2, 137–144. 1016/j.jaridenv.2014.01.017 Verhoeven, G. (2011) Taking computer vision aloft–archaeological three- Sankey, T.T., McVay, J., Swetnam, T.L., McClaran, M.P., Heilman, P. & dimensional reconstructions from aerial photographs with Nichols, M. (2018) UAV hyperspectral and lidar data and their fusion photoscan. Archaeological Prospection, 18(1), 67–73. Available from: for arid and semi-arid land vegetation monitoring. Remote Sensing in https://doi.org/10.1002/arp.399 Ecology and Conservation, 4(1), 20–33. Available from: https://doi. Verstraeten, W.W., Veroustraete, F., van der Sande, C.J., Grootaers, I. & org/10.1002/rse2.44 Feyen, J. (2006) Soil moisture retrieval using thermal inertia, deter- Sedia, E.G. & Ehrenfeld, J.G. (2006) Differential effects of lichens and mined with visible and thermal spaceborne data, validated for mosses on soil enzyme activity and litter decomposition. Biology European forests. Remote Sensing of Environment, 101(3), 299–314. and Fertility of Soils, 43(2), 177–189. Available from: https://doi.org/ Available from: https://doi.org/10.1016/j.rse.2005.12.016 10.1007/s00374-006-0077-6 Wang, L., Zhang, G., Zhu, L. & Wang, H. (2017) Biocrust wetting induced Sevgi, E., Yilmaz, O.Y. & Sevgi, O. (2020) Small-scale spatial patterns of change in soil surface roughness as influenced by biocrust type, cov- two terricolous lichens in a conifer plantation. Applied Ecology and erage and wetting patterns. Geoderma, 306, 1–9. Available from: Environmental Research, 18(2), 3557–3576. Available from: https:// https://doi.org/10.1016/j.geoderma.2017.06.032 doi.org/10.15666/aeer/1802_35573576 Weber, B. & Hill, J. (2016) Remote Sensing of Biological Soil Crusts at Dif- Short, N.M. & Stuart, L.M. (1982) The heat capacity mapping mission ferent Scales. In: Weber, B., Büdel, B. & Belnap, J. (Eds.) Biological Soil (HCMM) anthology, Vol. 465. Washington, DC: Scientific and Techni- Crust: An Organizing Principle in Drylands, Ecological Studies, 1st edi- cal Information Branch, National Aeronautics and Space tion, Vol. 226. Berlin: Springer, pp. 215–234. Administration. Weber, B., Olehowski, C., Knerr, T., Hill, J., Deutschewitz, K., Wessels, D.C. Smith, G.M. & Milton, E.J. (1999) The use of the empirical line method to J. et al. (2008) A new approach for mapping of biological soil crusts calibrate remotely sensed data to reflectance. International Journal of in semidesert areas with hyperspectral imagery. Remote Sensing of Remote Sensing, 20(13), 2653–2662. Available from: https://doi.org/ Environment, 112(5), 2187–2201. Available from: https://doi.org/10. 10.1080/014311699211994 1016/j.rse.2007.09.014 Smith, M.W. & Vericat, D. (2015) From experimental plots to experimental Weber, B., Wu, D., Tamm, A., Ruckteschler, N., Rodriguez-Caballero, E., landscapes: topography, erosion and deposition in sub-humid bad- Steinkamp, J. & Cheng, Y. (2015) Biological soil crusts accelerate the lands from Structure-from-Motion photogrammetry. Earth Surface nitrogen cycle through large NO and HONO emissions in drylands. Processes and Landforms, 40(12), 1656–1671. Available from: Proceedings of the National Academy of Sciences, 112(50), 15384– https://doi.org/10.1002/esp.3747 15389. Available from: https://doi.org/10.1073/pnas.1515818112 Smith, W.K., Dannenberg, M.P., Yan, D., Herrmann, S., Barnes, M.L., Westoby, M.J., Brasington, J., Glasser, N.F., Hambrey, M.J. & Reynolds, J.M. Barron-Gafford, G.A. et al. (2019) Remote sensing of dryland ecosys- (2012) Structure-from-Motion’photogrammetry: A low-cost, effective tem structure and function: Progress, challenges, and opportunities. tool for geoscience applications. Geomorphology, 179, 300–314. Remote Sensing of Environment, 233, 111401. Available from: https:// Available from: https://doi.org/10.1016/j.geomorph.2012.08.021 doi.org/10.1016/j.rse.2019.111401 Xue, R., Yang, Q., Miao, F., Wang, X. & Shen, Y. (2018) Slope aspect influ- Sörensen, R., Zinko, U., Seibert, J., Sørensen, R., Zinko, U. & Seibert, J. ences plant biomass, soil properties and microbial composition in (2006) On the calculation of the topographic wetness index: evalua- alpine meadow on the Qinghai-Tibetan plateau. Journal of Soil Science tion of different methods based on field observations. Hydrology and and Plant Nutrition, 18(1), 1–12. Earth System Sciences, 10(1), 101–112. Available from: https://doi. Zevenbergen, L.W. & Thorne, C.R. (1987) Quantitative analysis of land org/10.5194/hess-10-101-2006 surface topography. Earth Surface Processes and Landforms, 12(1), Stieglitz, M., Shaman, J., McNamara, J., Engel, V., Shanley, J. & Kling, G.W. 47–56. Available from: https://doi.org/10.1002/esp.3290120107 (2003) An approach to understanding hydrologic connectivity on the Zhou, X., Ke, T., Li, S., Deng, S., An, X., Ma, X. & Chen, L. (2020) Induced hillslope and the implications for nutrient transport. Global Biogeo- biological soil crusts and soil properties varied between slope aspect, chemical Cycles, 17(4), 1–15. slope gradient and plant canopy in the Hobq desert of China. Catena, Story, M. & Congalton, R.G. (1986) Accuracy assessment: a user’s perspec- 190(March), 104559. Available from: https://doi.org/10.1016/j. tive. Photogrammetric Engineering and Remote Sensing, 52(3), catena.2020.104559 397–399. Suggitt, A.J., Wilson, R.J., Isaac, N.J.B., Beale, C.M., Auffret, A.G., August, T. SUPPORTING INFORMATION et al. (2018) Extinction risk from climate change is reduced by micro- Additional supporting information may be found online in the climatic buffering. Nature Climate Change, 8(8), 713–717. Available from: https://doi.org/10.1038/s41558-018-0231-9 Supporting Information section at the end of this article. Tesfa, T.K., Tarboton, D.G., Watson, D.W., Schreuders, K.A., Baker, M.E. & Wallace, R.M. (2011) Extraction of hydrological proximity measures from DEMs using parallel processing. Environmental Modelling & Soft- How to cite this article: Blanco-Sacristán, J., Panigada, C., ware, 26(12), 1696–1709. Available from: https://doi.org/10.1016/j. Gentili, R., Tagliabue, G., Garzonio, R., Martín, M.P. et al. envsoft.2011.07.018 (2021) UAV RGB, thermal infrared and multispectral imagery Tarboton, D.G. (1997) A new method for the determination of flow direc- tions and upslope areas in grid digital elevation models. Water used to investigate the control of terrain on the spatial Resources Research, 33(2), 309–319. Available from: https://doi.org/ distribution of dryland biocrust. Earth Surface Processes and 10.1029/96wr03137 Landforms, 46(12), 2466–2484. Available from: https://doi. Tucker, C.L., Ferrenberg, S. & Reed, S.C. (2019) Climatic sensitivity of dry- org/10.1002/esp.5189 land soil CO2 fluxes differs dramatically with biological soil crust suc- cessional state. Ecosystems, 22(1), 15–32. Available from: https://doi. org/10.1007/s10021-018-0250-4