Spatiotemporal Analysis of Vegetation Cover Change in a Large Ephemeral River: Multi-Sensor Fusion of Unmanned Aerial Vehicle (UAV) and Landsat Imagery

: Ephemeral rivers in arid regions act as linear oases, where corridors of vegetation supported by accessible groundwater and intermittent surface ﬂows provide biological refugia in water-limited landscapes. The ecological and hydrological dynamics of these systems are poorly understood compared to perennial systems and subject to wide variation over space and time. This study used imagery obtained from an unmanned aerial vehicle (UAV) to enhance satellite data, which were then used to quantify change in woody vegetation cover along the ephemeral Kuiseb River in the Namib Desert over a 35-year period. Ultra-high resolution UAV imagery collected in 2016 was used to derive a model of fractional vegetation cover from ﬁve spectral vegetation indices, calculated from a contemporaneous Landsat 8 Operational Land Imager (OLI) image. The Normalized Difference Vegetation Index (NDVI) provided the linear best-ﬁt relationship for calculating fractional cover; the model derived from the two 2016 datasets was subsequently applied to 24 intercalibrated Landsat images to calculate fractional vegetation cover for the Kuiseb extending back to 1984. Overall vegetation cover increased by 33% between 1984 and 2019, with the most highly vegetated reach of the river exhibiting the greatest positive change. This reach corresponds with the terminal alluvial zone, where most ﬂood deposition occurs. The spatial and temporal trends discovered highlight the need for long-term monitoring of ephemeral ecosystems and demonstrate the efﬁcacy of a multi-sensor approach to time series analysis using a UAV platform.

Several recent studies have used UAVs to estimate FVC, either on their own or in conjunction with coarser-resolution satellite images, in environments ranging from the northern Fennoscandian Arctic [63] to the Qinghai-Tibetan Plateau [64,65] to Australian rangelands [66]. Approaches have included random forest regression and object-based image analysis [66], spectral unmixing [66,67], and upscaling of binary vegetated-unvegetated classifications from UAV imagery to satellite imagery using vegetation indices [63], an approach similar to the one used in this study. Several of these studies show that FVC estimates are more accurate when upscaled to 30-m resolution or coarser, with lower FVC accuracy at finer spatial resolution [63,64,66].
While most prior studies have used UAVs to estimate FVC in grasslands, shrublands, or agricultural environments (e.g., Reference [63,64,66,67]), this study uses UAV imagery to quantify woody vegetation cover in a desert riparian forest. Such research contributes to the growing understanding of the extent of dryland woody cover, which recent studies have shown to be much higher than previously estimated by using very high-resolution remotely sensed imagery [68,69].
The goal of this study was to quantify change in vegetation cover along a large ephemeral river in the Namib Desert using high-resolution UAV imagery to derive fractional vegetation cover from a time series of Landsat images. Such enhancement of satellite datasets with UAV-derived products, combined with empirical intercalibration of historical and current satellite imagery, provides a novel approach to temporal analysis using multi-resolution remote sensing techniques, adding to the emerging field of UAV-based remote sensing for environmental research, monitoring, and conservation.

Study Area
One of the Namib's twelve major ephemeral rivers, the 560-km-long Kuiseb, drains a predominantly arid watershed covering 14,700 km 2 as it flows from its headwaters in the Khomas Hochland westward to its outlet to the South Atlantic Ocean near Walvis Bay. The river has several tributaries as it cuts through rock sediments to form the Kuiseb Canyon, before becoming a sandy channel approximately 150 km from the coast and eventually widening to a delta at Rooibank, approximately 31 km from the coast ( Figure 1; Reference [3]). Rainfall in the lower Kuiseb basin, ranging from 15 to 50 mm per annum, is too little to contribute to river flow; floods are largely caused by precipitation upstream in the Khomas Hochland (approx. 350 mm per annum) and occur sporadically between December and April [3,5]. In its middle reaches, the vegetation and intermittent flood pulses of the river form a boundary between two distinct and unique ecosystems of the Central Namib Desert, the gravel plains to the north and the Namib Sand Sea to the south. The riparian corridor is home to most of the vegetation in this desert region, and the riparian forest is primarily comprised of five woody species: Acacia (Vachellia) erioloba, Faidherbia albida, Euclea pseudebenus, Tamarix usneoides, and Salvadora persica. Both A. erioloba and F. albida have been proposed as conservation priorities in Southern Africa [19,70] and produce pods that serve as the primary source of fodder for the livestock of the indigenous Topnaar ( =Aonin) people who live along the banks of the river [16]. Detailed information on the past floristic composition of the Kuiseb's riparian forest are provided by Seely et al. [5] and Theron et al. [51], while the lateral and longitudinal patterns in present tree species distribution within the riparian corridor are analyzed by Morgan et al. [34].
Our study area ( Figure 1) covered a 112-km section of the lower Kuiseb, extending from the lower end of the Kuiseb Canyon (approximately 168 km from the coast) to just beyond the Swartbank settlement (approximately 56 km from the coast). UAV data were collected at twelve 500-m-wide sites spaced every five kilometers within the middle reaches  Our study area (Fig. 1)  Although all of the images were downloaded as atmospherically corrected products, direct 138 comparison between images from different dates and sensors required additional intercalibration to 139 account for variations in incident and reflected energy, atmospheric effects, sensor response, and band 140 designations, which vary slightly between the Landsat sensors [72]. We performed relative radiometric 141 calibration (or radiometric normalization; [73,74]) to standardize reflectance values according to the 142 procedure proposed by Furby and Campbell [72] involving: (1) selection of a reference image; (2)

Intercalibration of Landsat Data
The Landsat archive provides global imagery from 1982 to the present at 30 m spatial resolution and is freely available [71], making it well-suited for quantifying vegetation change in regions with historically sparse in-situ monitoring programs, such as the Namib Desert. To extend the record as far back as possible, we used imagery obtained with the Landsat 5 Thematic Mapper (TM), Landsat 7 Enhanced Thematic Mapper Plus (ETM+; Scan Line Corrector-on), and Landsat 8 Operational Land Imager (OLI) sensors. Atmospherically corrected Landsat images (surface reflectance products) obtained between late May and early September that were unobstructed by clouds were downloaded from the USGS Earth Explorer website. In total, we analyzed 24 images collected between 1984 and 2019 ( Table 1). The study area comprised a subset of a single Landsat scene (WRS Path 179, Row 76).
Although all of the images were downloaded as atmospherically corrected products, direct comparison between images from different dates and sensors required additional intercalibration to account for variations in incident and reflected energy, atmospheric effects, sensor response, and band designations, which vary slightly between the Landsat sensors [72]. We performed relative radiometric calibration (or radiometric normalization; Reference [73,74]) to standardize reflectance values according to the procedure proposed by Furby and Campbell [72], involving: (1) selection of a reference image; (2) identification of at least 20-25 invariant targets spanning a range of brightness values; (3) calculation of calibration coefficients for each individual band; (4) evaluation of calibration curves and refinement of targets as necessary; and (5) application of the calculated coefficients for each band to each image. The 2016 Landsat 8 OLI image, contemporaneous with the UAV data, was used as the reference image, and reflectance values for 217 invariant pixels from 25 locations distributed across the study area were extracted and used to generate unique regression equations for the blue, green, red, and NIR bands of each image. The spectral shapes of the calibration coefficients (slope and intercept) and coefficient of determination (r 2 ) values were examined to compare the success of calibration across the four bands and detect differences between the sensors.
The resulting calibration equations were applied to each image by band to generate reflectance values standardized to the 2016 image.

Aerial Imagery Collection and Processing
The primary function of the UAV imagery in this study was to quantify fractional vegetation cover using natural color (RGB) and near-infrared (NIR) imagery collected between 15 July and 30 August 2016. The imagery was collected at twelve 500-m-long sites along the river, each spanning the width of the riparian zone. We used a Canon Powershot D30 digital camera with an internal GPS mounted in a nadir position on a 3D-Robotics IRIS+ quadcopter UAV to collect natural color imagery and a modified Canon Powershot A4000IS digital camera to collect NIR imagery. The modified Canon A4000IS used to collect the infrared imagery had its internal IR-blocking filter removed and a red filter placed over the lens, allowing red and NIR light to the sensor. To achieve stereoscopic images with at least 65 percent endlap and sidelap [75] and maximize both spatial resolution and total coverage area, the UAV was flown at an altitude of 70 m and a speed of 5 m s −1 . Both cameras were equipped with a Canon Hack Development Kit (CHDK) intervalometer script programmed to capture images every four (NIR) or five (RGB) seconds. The flight paths were planned in advance using ArduPilot Mission Planner software; the UAV was operated manually during takeoff and landing and autonomously during the collection phase. Based on the area of each site, between one and four missions were required to cover the entire riparian zone at each site. On average, each flight lasted approximately ten minutes and covered an area of approximately 5.1 hectares.
We processed the JPEG images from each mission in Agisoft PhotoScan Professional 1.3.0 using "high" quality settings (see Supplementary Materials for details). The processing steps included sparse point cloud reconstruction, dense reconstruction, production of a digital elevation model (DEM), and generation of an orthoimage mosaic. At each site, we identified eight to twelve ground control points (GCPs) on the RGB imagery (directly georeferenced using the camera's internal GPS) for each mission and used these GCPs to align and georeference the NIR orthoimages. All orthoimages were subsequently re-georeferenced to a 2010 WorldView-2 satellite image to correct for systematic errors and slight discrepancies between images from successive missions. A total of 32 RGB orthoimages with a ground sampling distance of 1.8 cm per pixel and 30 NIR orthoimages with a ground sampling distance of 1.5 cm per pixel were generated at 12 sites within the study area.

Calculation of Fractional Vegetation Cover from UAV Imagery
The UAV orthoimages were used to produce fractional cover raster datasets based on classification of vegetated and non-vegetated pixels. Initial attempts to separate the two classes based on spectral indices including NDVI and GRVI were unsuccessful. We developed an alternative false-red vegetation index (FRVI; Equation (1)) to create vegetation masks from each image based on two bands from the NIR camera.
Most consumer-grade digital cameras use silicon-diode charge coupled detectors (CCDs) that are sensitive to light wavelengths from about 350 to 1100 nm [76]. A Bayer pattern filter and an internal IR-blocking filter are used to obtain blue, green, and red bands for a normal image. The blue and red pixels of the sensor under the Bayer array are the most sensitive to NIR light; by removing the internal IR-blocking filter and adding a red-pass filter to the lens of the modified camera, the sensor's blue channel was able to collect only NIR reflectance, while the red channel recorded light over a range of red and NIR wavelengths, as the red-pass filter blocks green and blue light and transmits red and NIR light. With calibration and extensive post-processing, the NIR contributions on the red channel can be separated, allowing solely red-band reflectance to be resolved [76]. In this case, because RGB imagery was collected by a different sensor on separate flights, we were unable to isolate the red band of the modified camera's imagery. Incorporation of the broad red-IR band into our modified normalized difference index, however, proved more successful in isolating vegetation than the traditional indices. FRVI employs the broad red-IR band and the NIR band of this modified camera: Like other spectral indices, FRVI values range from −1.0 to +1.0, with high values for green vegetation and low values for sand and bedrock.
We calculated FRVI for all 30 NIR orthoimages. To reduce the memory usage and processing times, the pixel depth of the resulting FRVI images was rescaled from 32-bit (floating point) to 8-bit (integer). The rescaled FRVI was used to generate a mask of vegetation cover, with pixels classified as vegetation when they met both of the following two criteria: (a) having FRVI values greater than 136 (0.00667 non-rescaled) and (b) having NIR values above an adaptive threshold determined by the orthoimage's mean NIR value. Pixels with values equal to or below either threshold were classified as non-vegetation.

Calibration of the FVC Model for the Landsat Time Series
The high-resolution (cm-scale) masks of vegetation cover generated from the UAV orthoimages were aggregated to 30 × 30 m, producing a grid of fractional cover aligned with Landsat pixels (Figure 2). These FVC estimates were then used to model FVC from the Landsat imagery based on a spectral index. We evaluated five spectral indices for use as the predictor variable in this FVC model.  Perhaps the most widely used spectral index is the Normalized Difference Vegetation Index

215
(NDVI), defined as follows [36]: where ρ NIR and ρ RED are the reflectance of the NIR and red bands, respectively. NDVI effectively Perhaps the most widely used spectral index is the Normalized Difference Vegetation Index (NDVI), defined as follows [36]: where ρ NIR and ρ RED are the reflectance of the NIR and red bands, respectively. NDVI effectively measures greenness; it depends on the near-infrared (NIR) and red spectral bands, based on the principle that photosynthetic vegetation has high reflectance in the NIR portion of the electromagnetic spectrum and absorbs red light. NDVI ranges from -1.0 to +1.0, with high values resulting from high NIR and low red reflectance of healthy green vegetation, and negative values produced by surfaces with low NIR and high red reflectance, such as water or snow. Among NDVI's shortcomings, however, are its sensitivity to atmospheric effects, soil brightness, and lighting geometry [77]. A number of indices have been proposed to account Remote Sens. 2021, 13, 51 8 of 24 for these sensitivities including the Modified Soil-Adjusted Vegetation Index (MSAVI), given as: developed by Qi et al. [78] as an improved version of the SAVI, which was designed to mitigate soil brightness effects, specifically in arid and semi-arid environments; the Enhanced Vegetation Index (EVI), proposed by Huete et al. [77] as where G = 2.5, C 1 = 6, C 2 = 7.5, and L = 1, which incorporates the blue band to minimize scattering due to atmospheric effects and soil brightness; the Visible Atmospherically Resistant Index (VARI) [79], given as which exhibits a strong correlation with fractional vegetation cover and incorporates only the three visible bands; and the Green-Red Vegetation Index (GRVI; also known as the Visible-Green Index or VIG), given as which also displays a linear relationship with green vegetation cover [80]. To generate FVC models over the entire Landsat time series, we used the UAV-derived FVC data to perform pairwise linear regressions against each of these five spectral indices calculated from the contemporaneous 2016 Landsat 8 OLI image (Figure 2). Initial preprocessing of the satellite imagery involved cropping the reference image to the extent of each of the 30 NIR orthoimages and manual geo-rectification of these subsets to optimize the alignment between the Landsat data and the high-resolution UAV imagery.
FVC values were calculated from the UAV orthoimages for each corresponding 30 × 30 m Landsat pixel and regressed on spectral index values extracted from the 2016 Landsat image. Based on the expected linear relationships between the vegetation indices and FVC [38,80], linear regression analysis was performed on the entire dataset (1724 cells) to generate a model for calculation of FVC from each spectral index. Residual plots were examined to assess the fit of each model, and the optimal model was selected based on the index yielding the strongest correlation with FVC. The accuracy of the Landsat FVC model was assessed by comparing the high-resolution vegetation cover observations from each UAV mission to modeled FVC from the 2016 Landsat image for the corresponding area.
We calculated the selected spectral index for the entire study area for all 24 standardized Landsat images, including the 2016 reference image and applied the FVC model to these datasets, yielding fractional cover images of the entire study area dating back to 1984. All pixels with negative fractional cover values were set to 0, and those with values greater than 1.0 to 1.0 to constrain modeled FVC values to a physically realistic range. The riparian zone, clearly distinguishable from the surrounding landscape in this desert environment, was manually delineated on the natural-color 2016 Landsat image. Using this hand-digitized polygon, we then extracted FVC values from the Landsat datasets for each of the 215 contiguous 500-m segments.

Analysis of Vegetation Cover Change
We averaged the FVC values for all segments by year to provide an estimate of overall riparian vegetation cover for each year studied. We conducted a regression analysis to determine how the mean total cover of the Kuiseb River varied from 1984 to 2019. A more Remote Sens. 2021, 13, 51 9 of 24 detailed spatiotemporal analysis was also undertaken using the data from individual segments. Bivariate regressions of arcsine-transformed FVC [81,82] over time were performed for all 215 riparian segments. The slopes of the regression equations were plotted against distance upstream using a LOWESS regression curve to examine the relationship between FVC change and longitudinal position. Figures 3 and 4 show the calibration curves for the blue and NIR bands, respectively. All bands fit the calibration curves well, with the NIR bands showing the strongest correlation to the 2016 reference data and the blue bands exhibiting the most deviation.

Intercalibration of Landsat Data
The spectral shapes of the calibration coefficients and r 2 values are shown in Figure 5. There were no overall systematic differences detected between the sensors, though Landsat 8 OLI images typically had slopes closer to 1.0 and intercepts closer to 0 than Landsat 7 ETM+ images. Landsat 8 OLI images also had the highest r 2 values, as would be expected given that they were collected from the same sensor as the 2016 reference image. Regressions of 2016 reference values on Landsat 5 TM reflectance values yielded a wide range of slope and intercept values. The 1998 image had the lowest r 2 values for all bands, which we suspect was due to an error in its georeferencing.

Modeling of Fractional Vegetation Cover
The FVC regressions for the four spectral vegetation indices with the weakest relationships to FVC are shown in Figure 6. VARI and GRVI showed the poorest relationships to FVC with the lowest r 2 values and highest root-mean-square error (RMSE), and FVC appears to become saturated at low index values.
NDVI proved to be the best indicator of FVC (r 2 = 0.83, F 1,1723 = 8140.38, p < 0.001), with NDVI predicting cover according to the equation: The residuals appear to be normally distributed (Figure 7; RMSE = 0.12). As shown in Figure 8, the modeled 2016 Landsat FVC closely matches the observed UAV FVC (RMSE = 0.037, n = 30). Though this does not represent an independent estimate of the model's accuracy, it does illustrate the performance of the model in reconstructing the aggregate area of vegetation cover at the scale of analysis. Note that the 0.037 RMSE at this scale is substantially lower than the 0.12 RMSE at the scale of individual Landsat grid cells (900 m 2 ; Figure 7). Figure 9 shows a comparison of UAV-derived FVC and Landsat-derived FVC for 2016 at one site. As shown in this example, there was good agreement between the modeled and observed FVC across the full range of cover values.       This equation was applied to NDVI images calculated from the calibrated Landsat bands to model FVC over the entire study period. Over all years, 9.3% of Landsat pixels had NDVI values that translated to FVC values outside the 0 to 1 range, with considerably more low NDVI values (7.2%) than high (2.1%).

Change in Fractional Vegetation Cover
FVC within the entire 107-km length of the riparian corridor under study increased by 33% from 0.2298 (σ = 0.1557, n = 215) in 1984 to 0.3056 (σ = 0.1915, n = 215) in 2019. This trend was significant across the entire 35-year period examined (p < 0.001, n = 24; Figure 10) with moderate interannual variation (r 2 = 0.5195). Figure 11a shows FVC of the riparian zone in 500-m segments for each year. These data indicate that the portion of the river located roughly between 95 and 145 km from the coast was consistently the most highly vegetated. This reach also appears to be the most dynamic, undergoing the greatest change in FVC during the 35-year period studied. This observation is consistent with the quantitative analysis of FVC change by 500-m segment. The greatest rates of change were seen between 95 and 145 km upstream, with FVC increasing over time in these areas (Figure 11b). On the other hand, the stretch of the river between 75 and 95 km upstream from the coast was significantly less vegetated, and FVC values appeared to have remained fairly constant. In some areas, particularly upstream, the vegetation cover of the river decreased, as shown in Figure 11.

Discussion
Moderate-resolution optical satellite imagery from systems such as Landsat is freely available worldwide and provides a historical archive extending back to the 1980s, but its 30-meter spatial resolution is too coarse to directly resolve small shrub and tree crowns typical of riparian corridors. The centimeter-scale resolution of imagery from UAVs can provide precise estimates of fractional vegetation cover, but the extent of such imagery is constrained by limitations on UAV range and flying time, and historical data are generally not available at this scale. In this study, broader spatial and temporal coverage was obtained by modeling FVC in a Landsat 8 OLI image using spatially detailed vegetation cover data acquired from high-resolution UAV imagery, radiometrically normalizing archival Landsat imagery to match the base image, and extrapolating the FVC model over space (upstream and downstream from the UAV sites) and time (to prior and subsequent years). The result is a spatiotemporal dataset of sub-pixel FVC estimates for the continuous study area length of the river over a 35-year period. This approach takes advantage of the complementary strengths of both the satellite and UAV data sources, yielding a spatially and temporally extensive dataset with limited in-situ UAV imagery.
Though the Landsat imagery used were atmospherically corrected surface reflectance products, biases in each band's spectral reflectance estimates for individual images remained due to limitations of the atmospheric correction algorithms and the atmospheric optical properties data on which they are based. Figures 3-5 show that these effects can be modeled and the images radiometrically normalized using invariant targets, a process that is optimal in a landscape with flat, bright, spectrally stable targets, like the Namib. As sensors and correction algorithms continue to improve and more spatially explicit data on the optical properties of the atmosphere become available, the biases in standard surface reflectance products should diminish, making this normalization process less critical, except in the case of comparisons with lower-quality historical imagery.
There are a number of limitations to the UAV imagery obtained and analyzed in this study. The collection of images in the nadir position along parallel flight lines does not yield very accurate SfM reconstructions, which leads to geometric errors in the resulting orthoimages [55,83]. In order to produce orthoimages and DEMs with highly accurate threedimensional coordinates, images should be collected at off-nadir positions with overlapping flight lines for repeated coverage at different angles. Collection of ground control points in the field and/or successive flights flown at different altitudes and perpendicular flight lines between missions would also improve the distortion issues. Despite their shortcomings, however, the flight plans and UAV-based collection methods employed herein did not impede the intended analysis.
The introduction of the false-red vegetation index, necessitated by the lack of success using standard indices to separate vegetated and non-vegetated pixels of the UAV imagery, raises an important point about the weaknesses of very high-resolution aerial imagery. Spectral indices such as NDVI were developed for satellite and aerial images with coarser spatial resolution [36]. Therefore, NDVI may be less useful for higher spatial resolution datasets, such as those obtained via a UAV platform. In this case, NDVI was particularly problematic in the classification of F. albida crowns, which generally appear gray or brown on the natural color UAV imagery and lack the brighter green color of other species, such as A. erioloba and S. persica. As a result, NDVI values for pixels of F. albida crowns were often low, sometimes even lower than sand. A related consequence of the lack of green leaves on F. albida trees is that the canopy is discontinuous at this scale; the sparse, divergent branches make the sand below the canopy visible from the sky. With a ground sampling distance of less than 2 cm per pixel, a single crown contains a range of NDVI values that make classification of the entire crown difficult, even with object-based approaches. This issue was exacerbated by the fact that the NIR and red bands were obtained from two different sensors with different resolutions; even downsampling of the higher resolution NIR image from 1.5 cm to the 1.8 cm per pixel ground sampling distance of the RGB imagery yielded imperfect alignment of pixel edges. FRVI outperformed the other vegetation indices in this study for use with the highresolution UAV imagery, resulting in more accurate separation of vegetated and nonvegetated pixels. This index requires additional exploration to explain its success, however, and is specific to modified color-infrared digital cameras. The broad red-IR band used in FRVI covers a wider range of wavelengths than the two distinct red and NIR bands on Landsat and other multispectral sensors; it includes the red-edge and adjacent wavelength ranges that fall between the narrower red and NIR bands. FRVI may be a useful replacement for NDVI for images of other areas obtained from a similarly modified camera, circumventing the need to calibrate the red and/or green bands of the resulting image, but it would have to be examined in other locations and with different sensors. The superior performance of FRVI over NDVI in the UAV images may be largely attributable to the fact that both component bands were obtained from the same sensor, eliminating the need for realignment or resampling and minimizing the influence of co-registration errors. The success of FRVI in this study, however, indicates that the additional bands of such imagery may actually provide useful information, despite the fact that they are not collecting light over traditionally defined ranges of wavelength.
Our approach differs from that used by Yan et al. [67], who used color mixture analysis (CMA; similar to spectral mixture analysis) to estimate vegetation cover from the Hue-Saturation-Value (HSV) data derived from the raw RGB spectral bands of UAV-obtained imagery. This methodology produced highly accurate estimates of FVC (mean absolute error < 0.01) for two agricultural sites in Hebei, China [67]. This approach may be worth exploring in natural environments as an alternative to the FRVI method presented here, when traditional multispectral data are not available, though the aforementioned complexity of naturally occurring woody vegetation, such as F. albida crowns, and the bright, variable-color desert soils may render CMA-based methods less effective in this region than in an intensively managed agricultural setting. In dryland environments, with extensive portions of the landscape covered by bare soil or non-photosynthetic or sparsely green vegetation, FVC can be difficult to estimate, with best results obtained when up-scaling ultra-high resolution UAV data to coarser resolutions [66], similar to the 30-m scale used in this study.
While it proved ineffective for vegetation classification in the ultra-high resolution UAV imagery, NDVI prevailed as the optimal index for modeling FVC from the satellite imagery, highlighting the scale-dependence of spectral vegetation indices-the most favorable index at one scale is not necessarily the most favorable index at another. Perhaps more surprisingly, NDVI outperformed the four other indices, many of which were designed to compensate for NDVI's shortcomings, particularly in arid environments. A possible explanation for this could be the unique nature of the study area, with virtually no green cover except in the riparian forest. Furthermore, canopy background in the Kuiseb is dominated by bright riparian sands and, on the fringes of the riparian zone, underlying gravel plains and the red dune sand. Thus, EVI, which is highly sensitive to albedo, and visible indices, like VARI and GRVI, may lack the ability to effectively separate vegetation from sand using band differences. Furthermore, FVC became saturated at very low values of both VARI and GRVI, suggesting that a logarithmic model may be more suitable. Nonetheless, NDVI was strongly linearly correlated with FVC and therefore provided an effective mechanism for calculation of fractional cover for the Landsat time series.
As a point of comparison, Riihimäki et al. [63] likewise assessed multiple vegetation indices for modeling FVC from satellite imagery in an Arctic tundra ecosystem, using binary classification of high-resolution UAV data into vegetated and non-vegetated classes followed by aggregation to the satellite-pixel scale. While they did not compare VARI, GRVI, or MSAVI, they similarly found NDVI to perform better than EVI, though alternatives using a shortwave infrared (SWIR) or red-edge band were superior overall, suggesting potential consideration for future application in the Namib region. Chen et al. [64] compared NDVI, EVI, and MSAVI for use in modeling FVC in an alpine grassland site using UAV data and 30-m-resolution imagery from the HJ-1A satellite, though they found no significant differences between the vegetation indices.
Our analysis revealed a 33% increase in fractional vegetation cover in the middle reaches of the Kuiseb since 1984. In principle, this increase could be due to either a change in land cover, such as the expansion of the riparian forest into previously unvegetated areas, an intensification of vegetation cover within existing land cover classes, or a combination of the two. Riparian forest cover in this region consists of a heterogeneous mix of trees, shrubs, and graminoids intermingled with extensive areas of sand. There are no existing land-cover change datasets for the Kuiseb region at scales allowing for resolution of the riparian zone; the best available is the Moderate Resolution Imaging Spectroradiometer (MODIS) Land Cover Type (MCD12Q1; Reference [84]) data product, which provides global land cover data at fairly coarse 500-m resolution yearly, beginning in 2000. More importantly, however, land cover change processes that are common in other parts of the world (e.g., urbanization, agricultural development, deforestation, etc.) are minimal in the middle reaches of the Kuiseb under study. Furthermore, in this region, the Kuiseb riparian zone is constrained by large sand dunes on the left bank and rocky outcrops, particularly on the right bank, which are geographically consistent at decadal time scales. Thus, the increase in FVC can be attributed to an increase in riparian forest cover.
Our findings in the Kuiseb are consistent with the trend of increasing vegetation greenness across the globe over the last four decades. Many recent studies have demonstrated a "greening" trend in vegetation cover both globally [85][86][87] and regionally, from the Arctic [88] to the tropics [85]. While obvious increases in vegetation greenness are associated with anthropogenic land use and land cover change (LULCC), such as expansion of irrigated agriculture and afforestation, greening is also occurring in regions with less direct human influence [85,88]. Potential causes of such increases involve nitrogen deposition [87], CO2 fertilization, and climate change [85]. In some regions, short-term "browning" episodes have intermittently reversed the longer-term trend [86]. In Southern Africa, precipitation and water availability may be a controlling factor in regional greening and browning trends [85,89], which in turn have critical implications for the water cycle through changes in evapotranspiration [28,90].
In the Kuiseb, the observed increase in FVC is likely due to changes in cover of the three dominant woody species-A. erioloba, F. albida, and T. usneoides-as cover of other woody and herbaceous species is low [34]. Furthermore, while the neighboring Swakop River has seen significant dieback of F. albida attributable to increasing aridity due to water extraction and expansion of the invasive Prosopis glandulosa [18], few P. glandulosa individuals have been observed in the Kuiseb, and the current abundance of both F. albida and A. erioloba remain comparable to those in documented in the 1970s and 1980s [5,34,51]. The smaller T. usneoides was substantially more abundant in 2016 than four decades prior [5,34], though its overall cover is low compared to the larger A. erioloba and F. albida, so an increase in abundance of this species is unlikely to explain much of the increase in overall vegetation cover.
Given the hydrological controls on riparian vegetation in ephemeral systems, a change in flow and precipitation patterns would be the most obvious explanation for the increase in vegetation cover. All of the Kuiseb species rely on consistent surface flow and associated aquifer recharge for survival, as groundwater is their dominant water source [91]. While regular flooding is critical for the sustaining the riparian forest, strong floods can increase mortality via uprooting or waterlogging vegetation [92], and decreases in woody biomass have been observed following large flood events [93]. The increase in vegetation cover in the Kuiseb is consistent with a recent analysis of repeat photography across the Namib Desert that showed an increase in vegetation cover over the 20th century, potentially associated with an increase in precipitation in the latter half of the century throughout the region [94]. However, regional climate and hydrological records do not show a consistent increase in precipitation or flooding over our study period (1984-2019) that would explain the observed increase in vegetation cover. Data recorded at Gobabeb over the last three decades indicate that flood events may be approaching extremes, with consistent moderate floods fewer and farther between (Figure 12a). Since 1983, the flood has reached Gobabeb during all but five summers. The dry years all occurred in the last thirteen years, with no flow in the summers of 2007, 2010, 2013, 2016, and 2019. Furthermore, precipitation at the headwaters of the Kuiseb also shows no consistent trend, with a shift from below average precipitation to a wet period at the turn of the century, followed by recurrent droughts in recent years (Figure 12b). It is worth noting, however, that these hydrological records do highlight a period of prolonged drought in the early 1980s, during which surface flow did not reach Gobabeb for four consecutive years from 1980 to 1983. Significant decline of the water table and widespread mortality of mature riparian trees-most notably F. albida-were documented during this severe dry period [51,95]. Thus, the observed 33% increase in vegetation cover since 1984 could represent a long-term recovery from drought. This theory could be tested using coarser resolution satellite data and/or archival aerial imagery acquired prior to this drought period, both of which would require substantial calibration for comparison to present-day imagery.
Version December 10, 2020 submitted to Remote Sens.
18 of 25 All of the Kuiseb species rely on consistent surface flow and associated aquifer recharge for survival, as 444 groundwater is their dominant water source [93]. While regular flooding is critical for the sustaining 445 the riparian forest, strong floods can increase mortality via uprooting or waterlogging vegetation [94], 446 and decreases in woody biomass have been observed following large flood events [95]. The increase in 447 vegetation cover in the Kuiseb is consistent with a recent analysis of repeat photography across the 448 Namib Desert that showed an increase in vegetation cover over the 20th century, potentially associated 449 with an increase in precipitation in the latter half of the century throughout the region [96]. However, at the headwaters of the Kuiseb also shows no consistent trend, with a shift from below average 457 precipitation to a wet period at the turn of the century, followed by recurrent droughts in recent 458 years (Fig. 12b). It is worth noting, however, that these hydrological records do highlight a period  The spatiotemporal analysis of fractional vegetation cover revealed that the most highly vegetated reach of the river, located around 95 to 145 km from the coast, was also the most dynamic over the last three-and-a-half decades showing greater increases than other parts of the study area (Figure 11b). This may reflect a longitudinal response of the vegetation to Remote Sens. 2021, 13, 51 20 of 24 hydrological and geomorphological trends. Like in other dryland river systems, transmission losses in the Kuiseb result in much of the flood load being deposited in the middle reaches, in this case between 90 and 190 km from the mouth of the river [27,29,31,33]. Ninety percent of Kuiseb floods between 1981 and 2006 ran dry in the middle reaches of the Kuiseb around Homeb (approximately 120 km from the coast) [6], roughly in the middle of the study area. Jacobson et al. [27] found that silt content and levels of organic carbon, nitrogen, and phosphorous peaked in the soils of this terminal alluvial zone. These soil nutrients have important effects on production and decomposition, and silt traps moisture in shallow soil, preventing downward absorption. The stream gradient decreases in this stretch, reaching a minimum in the most highly vegetated reach between 95 and 145 km, before becoming increasingly steep farther downstream [34], where vegetation cover is lower. That this stretch is the most vegetated and shows the highest increases in cover suggests that the terminal alluvial zone may be the most ecologically resilient river reach, sustaining the riparian forest throughout the recurrent droughts of the last thirteen years. The observed ecological response in this region is consistent with the geomorphological and hydrological activity in the terminal alluvial zone.

Conclusions
The combination of high spatial resolution UAV imagery and high temporal resolution satellite imagery facilitated the analysis of the structure of the Kuiseb riparian vegetation across longitudinal and temporal gradients, revealing an overall increase in fractional vegetation cover of 33% over the last three-and-a-half decades. Furthermore, fractional cover was highest in the stretch where the river runs dry, highlighting the importance of this terminal alluvial zone in controlling the structure of the riparian forest. That this area also exhibited the most change in vegetation cover between 1984 and 2019 raises important questions about how the Kuiseb woody species will respond to a potential increase in hydrologic variability in the future [97].
The observed spatiotemporal variability of the Kuiseb River highlights the importance of understanding the processes influencing riparian vegetation growth and establishment in ephemeral ecosystems, particularly in light of potential increases in hydrologic variability, expansion of invasive species, and anthropogenic alteration of the hydrologic regime. Riparian corridors act as thermal and moisture refugia for a range of plant and animal species, particularly in water-limited landscapes. Changes in riparian vegetation cover could therefore threaten the diversity of life these landscapes support. This, in addition to the downstream effects of changes in riparian forest structure on the geomorphology and hydrology of the riparian landscape in general, underscores the conservation significance of such ecosystems. Approximately one-third of the Earth's land surface is classified as arid or semi-arid [4], making ephemeral hydrology extremely important to a large portion of the world's biodiversity and human population. Furthermore, because desert rivers, like the Kuiseb, are adapted to such extremely variable conditions, understanding the behavior of and changes in such systems may provide valuable insights for the future of riparian forests. The changes observed in this study emphasize the utility of fusing highly spatially and highly temporally resolved data to elucidate the complex ecological dynamics in these highly variable arid river systems.