Characterizing natural variability in complex hydrological systems using passive microwave-based climate data records: a case study for the Okavango Delta

The Okavango River system in southern Africa is known for its strong interannual variability of hydrological conditions. Here, we present how this is exposed in surface soil moisture, land surface temperature, and vegetation optical depth as derived from the Land Parameter Retrieval Model, using an inter-calibrated, long-term, multi-sensor passive microwave satellite data record (1998–2020). We also investigate how these interannual variations relate to state-of-the-art climate reanalysis data from ERA5-Land. We analysed both the upstream river catchment and the Okavango delta, supported by independent data records of discharge measurements, precipitation, and vegetation dynamics observed by optical satellites. The seasonal vegetation optical depth anomalies have a strong correspondence with the MODIS leaf area index (correlation catchment: 0.74, delta: 0.88). Land surface temperature anomalies derived from passive microwave observations match best with those of ERA5-Land (catchment: 0.88, delta: 0.81) as compared to MODIS nighttime land surface temperature (LST) (catchment: 0.70, delta: 0.65). Although surface soil moisture anomalies from passive microwave observations and ERA5-Land correlate reasonably well (catchment: 0.72, delta: 0.69), an in-depth evaluation over the delta uncovered situations where passive microwave satellites record strong fluctuations, while ERA5-Land does not. This is further analysed using information on inundated area, river discharge, and precipitation. The passive microwave soil moisture signal demonstrates a response to both the inundated area and precipitation. ERA5-Land however, which, by default, does not account for any lateral influx from rivers, only shows a response to the precipitation information that is used as forcing. This also causes the reanalysis model to miss record low land surface temperature values as it underestimates the latent heat flux in certain years. These findings demonstrate the complexity of this hydrological system and suggest that future land surface model generations should also include lateral land surface exchange. Also, our study highlights the importance of maintaining and improving climate data records of soil moisture, vegetation, and land surface temperature from passive microwave observations and other observation systems.

1 Introduction

Long-term data records of key components of the climate system, known as essential climate variables (ECV), are important for improving our understanding and predictability of climate behaviour at different timescales (Hollmann et al., 2013; Bojinski et al., 2014). These records can help us to determine the root causes of observed climate change, e.g. natural or anthropogenic, assess its impacts and associated risks, and support mitigation and adaptation activities. In 2008, the European Space Agency (ESA) started the Climate Change Initiative programme (CCI) to develop these ECVs from satellite data records. This was done in response to the United Nations Framework Convention on Climate Change (UNFCCC) need for systematic monitoring of the climate system. Today, the CCI programme covers 21 satellite-based ECV records (projects https://www.esa.int/, last access: September 2021).

Surface soil moisture (SSM) is one of these ESA CCI ECVs. These records are based on a fusion of both passive (PMW) and active microwave satellite retrievals (Dorigo et al., 2017). The current version 6.1 spans from 1979 until 2020 (Scanlon et al., 2021), and contains three separate SSM products, which are derived from active, passive, and a combination of active and passive sensors. The methodology and evaluation of the harmonization and merging of the soil moisture retrievals from multiple satellites are described by Gruber et al. (2019). ESA CCI SSM data have been used for more than 10 years as the baseline for the annual evaluation and interpretation of global SSM conditions, as reported in the leading “State of the Climate” reports (Van der Schalie et al., 2021a) that are published as a Supplement to the Bulletin of the American Meteorological Society. Three datasets are produced as part of the passive input for the ESA CCI SM, which is SSM (SSM_MW), but also land surface temperature (LST_MW), and vegetation optical depth (VOD_MW).

SSM_MW datasets have been extensively evaluated with ground observations, models, other satellite products, and related ECVs, like precipitation (e.g. Hirschi et al., 2021; Beck et al., 2021; Dorigo et al., 2015; Al-Yaari et al., 2019; Albergel et al., 2013; Loew et al., 2013). VOD_MW has been used in multiple studies with a focus on seasonal and interannual vegetation dynamics (e.g. Liu et al., 2015; Moesinger et al., 2020; Teubner et al., 2019) or specifically on L-band VOD characteristics (e.g. Schwank et al., 2021; Bousquet et al., 2021; Rodríguez-Fernández et al., 2018). Research on the quality of LST_MW (e.g. Holmes et al., 2009, 2015) remains limited. The robustness of the interannual variability signals within these multi-decadal data records is still not always clear and a combined assessment of all three variables is necessary for understanding these datasets, as the current joint retrieval algorithm makes their values fundamentally intertwined (Owe et al., 2008). Such information provides unique opportunities for both monitoring and seasonal forecasting, e.g. over Africa (e.g. Cook et al., 2021).

The purpose of this paper is to improve insight into the interannual signals of the SSM_MW, LST_MW and VOD_MW by presenting a case study over a region with a complex hydrological system, i.e. the Okavango, and how their skill compares to state-of-the-art climate reanalysis data from ERA5-Land (Muñoz Sabater, 2019; Muñoz-Sabater et al., 2021). ERA5-Land aims to quantify the water and energy cycles over land in a consistent manner, therefore, allowing the characterization of trends and anomalies. Although ERA5-Land (E5) data are known to be of high quality in many regions around the globe, for use in any specific regions, these need to be properly evaluated. Therefore, this dataset does not only function as a benchmark in this study, but will also be analysed in more detail to evaluate its ability to properly detect the natural dynamics and variability in the Okavango and how this compares to the signal of the passive microwave-based datasets. Other datasets are used as support for determining which dataset (i.e. either PMW or E5L estimates of the same variable) is more likely to reflect true conditions. This research can help to improve the synergy between Earth Observation (EO) datasets and land surface models, and to identify both strengths and shortcomings of either one.

More specifically, the Okavango delta and Okavango River catchment in southern Africa were selected as the study areas. The Okavango delta (Republic of Botswana, 2013) consists of permanent marshlands and seasonally flooded plains, and is one of the few endorheic “delta” systems (geomorphologically, Okavango delta is an alluvial fan, Kgathi et al., 2006) that does not flow into the ocean. It is an exceptional example of the interaction between climatic, hydrological, and biological processes leading to a unique mix of flora and fauna, and has therefore been included in the UNESCO World Heritage List since 2014. Three features in the local hydrological system stand out, i.e. the strong interannual variability, the lateral water influx component of the Okavango River into the delta, and the seasonal characteristics with a lag between rainfall, river discharge, and flooding. Unfortunately, it is expected that global warming will affect this natural variability in the hydrological cycle over the Okavango delta (Wolski et al., 2012, 2014), for example, reducing high-water periods like in 2009–2011. These kinds of negative impacts increase the need for reliable monitoring capabilities.

The structure of the paper is as follows. Section 2 introduces the study area and includes the exact regions of interest (ROIs) that are used for the data extraction. Section 3.1 describes the passive microwave data and other data sources. Section 3.2 explains the methodology concerning the inter-calibration (Sect. 3.2.1), the Land Parameter Retrieval Model (LPRM, Sect. 3.2.2), the evaluation of the dataset anomalies (Sect. 3.2.3), and the evaluation of the river, flood, and precipitation contribution to SSM anomalies over the Okavango delta (Sect. 3.2.4). Sections 4–6 provide the Results, Discussion and Conclusions of these different steps.

2 Research area

With a length of approximately 1600 km, the Okavango River is one of the largest in southern Africa (Muzungaire et al., 2012). The river is known globally for its large terminating inland “delta”. The Okavango delta is a large seasonally pulsed inland wetland, a mixture of aquatic vegetation, open water, and dry land with an actively inundated area covering a part of the 28 000 km² alluvial cone (Ringrose et al., 1988).

Figure 1 The research area comprising ROI₁ (a part of the upstream area of the Cubango and the Cuito rivers) and ROI₂ (the surrounding of the Okavango delta) in relation to the Okavango drainage basin (grey). The black dot marks the location of the discharge station at Mohembo.

In line with both the interannual variation in local and upstream rainfall and the longer term effects of surface–groundwater interactions, substantial interannual variability in the delta's inundated area was recorded over the period 1932–2000 (Wolski and Murray-Hudson, 2008), with annual minima of about 3000 km² up to annual maxima of 12 000 km² (Wolski et al., 2017; Gumbricht et al., 2004; McCarthy et al., 2003). Whereas estimates for the total annual water budget stemming from direct rainfall in the Okavango delta ranges between 25 % to 50 %, the Okavango River inflow accounts for the other 50 % to 75 % (McCarthy et al., 1998, 2000; Ashton, 2000; Ashton and Neal, 2003; Wolski et al., 2006).

In this study, we focus on only two perennial rivers in the Okavango catchment – the Cubango River and the Cuito River (Ashton and Neal, 2003). Data were extracted from the catchment area within ROI₁ of Fig. 1. These rivers originate in Angola and are a vital lifeline to the Okavango delta, with an average yearly inflow at Mohembo of 9863 millions of m³ (approximately 300 m³ s⁻¹) in the period 1932–2001 and a 71.4 % contribution to the total water budget of the delta.

The Angolan part of the basin is characterized by a subtropical climate, while in Botswana and Namibia, parts are classified as semi-arid (Kgathi et al., 2006). During drought years in the 1980s and 1990s, the annual inflow at Mohembo reduced up to 45 % (McCarthy et al., 2000; Ashworth, 2002; Ashton, 2003; Ashton and Neal, 2003), which then coincided with proportional declines of the Okavango delta outflow to the Thamalakane and Boteti rivers (Ashton, 2000; Ashworth, 2002; Ashton and Neal, 2003). Throughout these periods, a growing water demand arose in Botswana and Namibia (MGDP, 1997; Ashton, 2003). Overall, the dry phase was caused by multi-decadal oscillations in rainfall and was likely related to processes of internal variability in the climate system (Wolski et al., 2012).

ROI₁ and ROI₂ were chosen to study how their significantly different water influxes affect the signal of the data sources used in the evaluation. The delta is of particular interest as it is mostly driven by a strong and highly variable lateral influx from the Okavango River that creates a pattern of seasonally varying wetness that is asynchronous or off-phase with the rainy season.

3 Material and methods

3.1 Data

3.1.1 Passive microwave observations

The three main variables that are used for the analysis are surface soil moisture (SSM_MW), vegetation optical depth (VOD_MW), and land surface temperature (LST_MW). These variables are derived from passive microwave observations from multiple satellite sensors that observe in similar frequencies and overlap in time.

Table 1 Overview and characteristics of passive microwave satellite sensors used in the study.

Table 2 Results of the multiple linear regression for estimating the relationships between the z-score anomalies of SSM, PR, ORD, and ODIA_MD over ROI₂.

The Advanced Microwave Scanning Radiometer for EOS (AMSR-E, Kawanishi et al., 2003) is a 12-channel, six-frequency, passive microwave radiometer developed by the Japan Aerospace Exploration Agency (JAXA) and was active between 2002 and 2011. AMSR-E is part of the payload carried onboard the Aqua (EOS PM-1) NASA scientific research satellite, which has a polar orbit with a 01:30/13:30 LT equatorial crossing time for ascending/descending swaths. AMSR-E was launched to obtain data to improve our understanding of global-scale water and energy cycles and played a key role in the development of soil moisture retrieval algorithms. For the technical specifics, see Table 1. Only descending brightness temperature data were used for this study, as due to the thermal equilibrium during nighttime, these are more stable and of higher quality (Owe et al., 2008; Van der Schalie et al., 2021b).

The Advanced Microwave Scanning Radiometer 2 (AMSR2, Imaoka et al., 2012) onboard the GCOM-W1 satellite is the follow-up of AMSR-E and was launched in 2012. Although incorporating improvements, both the general setup and orbital characteristics (e.g. overpass times) are similar to AMSR-E (see Table 2). However, there is a gap between AMSR-E and AMSR2 of about 9 months, making a direct intercalibration of time series complicated.

To overcome this gap and to extend the passive microwave observation record back to 1998, we make use of the Tropical Rainfall Measuring Mission's (TRMM, Kummerow et al., 1998) Microwave Imager (TMI). TMI observes in X-band and higher frequencies. TRMM is not in a polar orbit because of its focus on the tropical regions and, therefore, does not cover the entire globe. Data are only available between 40 ^∘N and 40 ^∘S, and due to its orbital characteristics, have a variable crossing time (see Table 1). To find a good balance between data availability and data stability (the more stable temperature distribution between the soil and vegetation close to thermal equilibrium), only brightness temperature data were used that had a local overpass time between 22:30 and 04:30 LT, to best match AMSR-E and AMSR2.

For this study, we use VOD_MW and SSM_MW derived from X-band brightness temperature data due to its availability on all three sensors, while Ka-band is the main frequency used for the LST_MW. All brightness temperatures were collected and gridded into a 0.25^∘ grid for the study area.

3.1.2 Ancillary datasets

In our analysis, we use several ancillary datasets to determine the ability of passive microwave-based satellite data records to correctly capture interannual variations. These ancillary datasets are split into two types.

Firstly, data were used from the ERA5-Land climate reanalysis model (Muñoz Sabater, 2019; Muñoz-Sabater et al., 2021), which is an enhanced resolution (9 km × 9 km) land-only offline rerun of the ECMWF ERA5 climate reanalysis (Hersbach et al., 2020). SSM_E5, LST_E5, and PR_E5 were extracted. For both SSM_E5 and LST_E5, the Layer 1 (0–7 cm) was used. LAI_E5 was excluded from the analysis as it only contained a climatology based on satellite EOs (no interannual variability). ERA5-Land data were extracted from the Copernicus Climate Change Service (C3S) Climate Data Store (CDS). As it has an hourly resolution, the values closest to the satellite overpasses were chosen. Data cover the complete period from 1998 to 2020.

Secondly, independent observational datasets are used, which have the sole purpose of functioning as a benchmark. These consist of the Okavango River discharge measurements (ORD, Okavango Research Institute, 2021), Okavango delta inundated area (ODIA_MD), leaf area index (LAI_MD, Yang et al., 2006) and nighttime LST (LST_MD, Wan, 2014) from the Moderate Resolution Imaging Spectroradiometer (MODIS), and precipitation from the NASA Global Precipitation mission's Integrated Multi-satellitE Retrievals for GPM (IMERG) product (PR_IM, Huffman et al., 2015).

A majority of the water entering the Okavango delta originates from the delta inlet at Mohembo. Therefore, we use ORD from the Mohembo station (see Fig. 1) to indicate the long-term variability of the lateral inflow into the delta. Measurements using E-type gauge plates are done on a regular (fortnightly) basis by the Botswana Department of Water Affairs and the data are shared by the Okavango Research Institute of the University of Botswana. The advantage of using this dataset is that it has a long historical record dating back to 1974. For this study, data were extracted for the 1998–2020 period.

ODIA_MD represents the inundated area in the Okavango delta and is derived from shortwave infrared (SWIR) observations from the MODIS sensor (Wolski et al., 2017). More specifically, data for band b7 from the MCD43A4 product were used. Reflectances of training areas are used to dynamically determine the threshold used for the derivations of the inundation. An automated and up-to-date monitoring tool for the flooding extent can also be found online (http://www.okavangodata.ub.bw/, last access: 22 May 2022).

The LAI_MD is defined as the one-sided green leaf area per unit of ground area (Chen and Black, 1992; Yang et al., 2006). The LAI_MD for the study area, including both the drainage catchment and the delta, was extracted from the MOD15A2H Version 6 MODIS dataset. This is an 8-daily product that uses the best available pixel within the 8 d period. The product has a spatial resolution of 500 m and the mean was extracted for the complete ROIs.

One kilometre nighttime, about 01:30 LT surface temperature from MODIS, was extracted from the MYD11A2.006 product, which is based on the average over 8 d of all available LST_MD observations. For this study, the mean values of the two areas were extracted. The temporal coverage is from February 2000 to the end of 2020 for the LAI_MD and July 2002 to the end of 2020 for the LST_MD.

For PR_IM, data were used from the Integrated Multi-satellitE Retrievals for GPM (IMERG, Huffman et al., 2015), which are produced at 0.1^∘ resolution. IMERG is a unified algorithm that provides rainfall estimates based on a combination of observations from multiple passive microwave sensors, infrared sensors, and precipitation gauges. Mean daily data were used from the GPM_3IMERGDF version 6, covering June 2000 to December 2020.

3.2 Methods

3.2.1 Intercalibration of PMW brightness temperatures

The intercalibration of AMSR-E, AMSR2, and TRMM is based on the methodology described in Van der Schalie et al. (2021b). In this approach, a two-step linear regression model is used, which first defines a global slope and afterwards a local intercept. Secondly, it uses a cost function that not only minimizes the differences between brightness temperatures of the individual polarizations, i.e. vertical (V) and horizontal (H), but also for the ratio between the two. This is because the Land Parameter Retrieval Model (LPRM, see next section) used for the SSM_MW, VOD_MW, and LST_MW retrievals is very sensitive to the polarization ratio. Inconsistencies in this ratio between different sensors can lead to an imbalance in how the radiative transfer model distinguishes between the emission from the soil and vegetation, respectively, leading to biases in the resulting retrievals. This intercalibration methodology was previously only applied only to the Ku-, K-, and Ka-band, but is here also used for application to the X-band data in the same way.

More specifically, the following cost function is minimized in the linear regression instead of a standard least squares approach:

$$\begin{array}{}\text{(1)}& \begin{array}{rl}\mathrm{Err}& =\sum \mathrm{RMSE}{\mathrm{TB}}_{\mathrm{H}}+\sum \mathrm{RMSE}{\mathrm{TB}}_{\mathrm{V}}\\ & +\sum \mathrm{RMSE}\mathrm{MPDI},\end{array}\end{array}$$

with:

$$\begin{array}{}\text{(2)}& {\displaystyle }\mathrm{RMSE}{\mathrm{TB}}_{\mathrm{H}/\mathrm{V}}=\sqrt{{\displaystyle \frac{\sum _{t=\mathrm{1}}^T\left({\mathrm{TB}}_{\mathrm{H}/\mathrm{V}}^{s_{\mathrm{1}}}-\left(\mathit{\alpha }\cdot {\mathrm{TB}}_{\mathrm{H}/\mathrm{V}}^{s_{\mathrm{2}}}+\mathit{\beta }\right)\right)}{T}}}\text{(3)}& {\displaystyle }\begin{array}{rl}& \mathrm{RMSE}\mathrm{MPDI}=\\ & \sqrt{{\displaystyle \frac{\sum _{t=\mathrm{1}}^T\left(\frac{{\mathrm{TB}}_{\mathrm{V}}^{s_{\mathrm{1}}}-{\mathrm{TB}}_{\mathrm{H}}^{s\mathrm{1}}}{{\mathrm{TB}}_{\mathrm{V}}^{s_{\mathrm{1}}}+{\mathrm{TB}}_{\mathrm{H}}^{s_{\mathrm{1}}}}-\frac{\left(\mathit{\alpha }\cdot {\mathrm{TB}}_{\mathrm{V}}^{s_{\mathrm{2}}}+\mathit{\beta }\right)-\left(\mathit{\alpha }\cdot {\mathrm{TB}}_{\mathrm{H}}^{s_{\mathrm{2}}}+\mathit{\beta }\right)}{\left(\mathit{\alpha }\cdot {\mathrm{TB}}_{\mathrm{V}}^{s_{\mathrm{2}}}+\mathit{\beta }\right)+\left(\mathit{\alpha }\cdot {\mathrm{TB}}_{\mathrm{H}}^{s_{\mathrm{2}}}+\mathit{\beta }\right)}\right)}{T}}},\end{array}\end{array}$$

where TB is the brightness temperature for the two polarizations and from the base (s₁) and calibrated (s₂) satellite. The α and β are the slope and intercept for the linear regression. The T refers to the time steps with overlapping observations for a single location.

After retrieving sensor-specific SSM_MW, VOD_MW and LST_MW from the inter-calibrated TB, a linear regression is applied between the different sensors using their respective overlap to remove any leftover inconsistencies. Improved inter-calibration between sensors can lead to a reduced need for break corrections (e.g. Preimesberger et al., 2020) and help to better address related issues at the source.

As this study focuses on anomalies at a seasonal timescale, the temporal coverage obtained by the current three sensors is sufficient. However, as was shown by Van der Schalie et al. (2021b) and as is done for the passive microwave-based data input for the ESA CCI SM, other sensors like GPM, FengYun-3B, and FengYun-3D can be included without issues, resulting in improved revisit times and coverage.

3.2.2 Land Parameter Retrieval Model

The Land Parameter Retrieval Model (LPRM, Owe et al., 2008) is a retrieval algorithm that simultaneously solves for SSM_MW, VOD_MW, and LST_MW without the use of any ancillary data sources on vegetation or temperature. The model is based on the tau omega (τ–ω) model (Mo et al., 1982), which simulates the top-of-atmosphere brightness temperatures by modelling the individual contribution of the soil, vegetation, and atmosphere. LPRM mainly distinguished itself from other algorithms through the analytical derivation of the VOD (Meesters et al., 2005) and the use of Ka-band observations for the LST_MW (Holmes et al., 2009). Here, we use version 6.0 of LPRM, as developed by Van der Schalie et al. (2017).

LPRM is currently the main algorithm used for all the passive microwave-based SSM retrievals in ESA CCI SM (Dorigo et al., 2017). Due to its unique analytical solution for the derivation of VOD_MW that uses no external source of information for the vegetation, LPRM has also been used in several studies of long-term vegetation dynamics (Liu et al., 2012, 2015), land degradation (Liu et al., 2013; Van Marle et al., 2017), and the development of climate data records of VOD_MW, like the VOD Climate Archive (VODCA, Moesinger et al., 2020).

3.2.3 Evaluation of anomalies

To have a better understanding of the quality of the different datasets in detecting interannual variability and anomalies, a two-step comparison analysis is done. First, the anomalies are visualized over time and their dynamics assessed. Second, the relations between related datasets are quantified using correlations and visualized using density plots. This is done separately for the catchment and the delta.

The SSM_MW is compared to the SSM_E5, both representative of the moisture conditions in the first few centimetres of the soil. As this is a direct comparison, in this step, the focus will be on their similarity and differences without analysing what causes them. Additionally, an extensive analysis is conducted (Sect. 3.2.4) to determine which of the datasets most likely reflects the ground conditions based on their relation to ORD, ODIA_MD, and PR.

For VOD_MW, there is a comparison with another regularly used satellite-based dataset, LAI_MD. Theoretically, the VOD_MW represents the attenuation of the microwave emission through the vegetation cover, which is related to both the structure and moisture content of the vegetation. The LAI_MD is representative of the projected single-sided green leaf area per unit ground area. Although VOD and LAI are fundamentally different, it is assumed that for dynamic and sparsely to moderately vegetated regions, i.e. excluding forests, the X-band also mostly measures the response of the leaves with the microwave signal via the vegetation water content (Jackson and Schmugge, 1991). Further defining the quality and ability of VOD_MW to detect interannual variability can be especially useful in improving the applicability and understanding of independent vegetation data records based on passive microwave observations like VODCA (Moesinger et al., 2020).

Here, the anomalies of LST from three different sources, e.g. passive microwave (LST_MW), model (LST_E5), and thermal infrared (LST_MD) are evaluated. These represent slightly different parts of the soil surface, being 0–7, 0–0.1, and 0 cm (skin), respectively. The mismatch in depth is also a reason why we choose nighttime comparisons as there is much more thermal stability. When looking over longer periods (e.g. weeks, months), we assume that the slightly different definitions of the soil temperature should still show a similarity in underlying anomalies, as in correlating well.

Because the focus is on the (seasonal) variability over a multi-decadal time span, a 91 d moving average (±45 d) is first applied to the datasets. The climatology for the anomaly calculation is based on the 2003–2020 period, as the LST_MD is only available from 2003 onward and the overall consistency for the baseline is preferred. As the window for the moving average is 91 d, little impact is assumed from data loss due to cloud cover in the MODIS datasets.

It is worth keeping in mind that none of these datasets provide the “truth” or measure exactly the same quantity, therefore, differences are to be expected. In the analysis component (see the following section), extra attention will be given to a specific case in the Okavango delta where a clear divergence is observed between the different SSM datasets.

3.2.4 Analysis of river flooding and precipitation contribution to soil moisture anomalies in the Okavango delta

As further in this study (Sect. 4.1) the signal of the two SSM data sources (SSM_MW and SSM_E5) is shown to diverge over the Okavango delta, an in-depth analysis is set up to explain the main drivers of their respective signals. This can help to better understand what the SSM datasets represent and give users insight into how to use them for their research activities and applications.

The first step is to directly compare the SSM datasets to the ORD, ODIA_MD, and both PR_E5 and PR_IM. These datasets can provide insight into what is the driver of the SSM anomalies in this region. As described in Sect. 2, about 50 %–75 % of the total influx of water into the Okavango delta comes from the ORD, while PR, on average, contributes 25 %–50 %, so we expect to see this reflected in the SSM either via the ORD or the ODIA_MD signals.

Following this, a multiple linear regression exercise is conducted. This is done to look into the influence of the ODIA_MD, ORD, and PR signals on the SSM_MW and SSM_E5 anomalies in the delta. This allows us to determine the drivers of the SSM anomalies and more importantly, how they differ between the two. Instead of using the absolute anomalies in this analysis, the z-score is preferred as this normalization removes issues with conversion of units and can be interpreted as standardized anomalies. A visualization will also be made of the climatologies from the different datasets, including their 10 and 90 percentiles, to further define the connection and time lag between the signals of the different parameters.

4 Results

4.1 Soil moisture

Figure 2 shows the anomalies of SSM_MW and SSM_E5 over the Okavango catchment and delta, with Fig. 2e and f comparing them directly to each other in a density plot. In both areas, the two datasets correlate moderately well, 0.72 and 0.70, respectively. In the delta however, a mismatch occurs on some occasions, leading to a visible flat line in the density plot where the anomalies of SSM_MW vary, while the anomalies of SSM_E5 are close to zero (Fig. 2e). The SSM_MW anomalies over the catchment and SSM_E5 anomalies over both the catchment and delta seem to have clear short-term variability, as can be seen from the peaks in the wet season, while the dry season remains mostly stable around zero. Only the anomalies of SSM_MW over the delta diverge from this and show a more multi-year variation, with highs in the years around 2011 and lows in the early and late periods of the time period. These cases will be further analysed in Sect. 4.4 in combination with the ORD and PR.

Figure 2 SSM_MW and SSM_E5 absolute anomalies over the Okavango catchment (a, b, e), with the intensity of the colouring based on the z-score of the positive (blue) and negative (red) anomalies and the Okavango delta (c, d, f) in time series and density plots. A daily time step is used from the moving average dataset.

The absolute range of the anomalies differs to some extent between the two products: SSM_MW anomalies range between −0.03 and 0.025 m³ m⁻³ in the catchment and −0.05 and 0.05 m³ m⁻³ in the delta, whereas SSM_E5 anomalies range between −0.10 and 0.06 m³ m⁻³ in the catchment and −0.08 and 0.10 m³ m⁻³ in the delta. However, the dynamics of the signal are very similar.

4.2 Vegetation optical depth

Figure 3 shows the anomalies of VOD_MW and LAI_MD over the Okavango catchment and delta, with Fig. 3e and f again showing a direct comparison in a density plot. The two datasets have a 0.74 correlation over the catchment and up to 0.88 in the delta. Generally, a similar pattern is visible for both regions. One exception can be seen during the 2008 to 2011 period in the catchment, where the VOD_MW anomaly remains high throughout multiple years, while the overall above average LAI_MD anomalies fluctuate to a greater extent. The lowest values in the delta were detected early in the study period, with VOD_MW recording an almost −0.08 anomaly during 1998 and 2003. This 2003 event is also seen in the LAI_MD dataset, while no data are available for 1998. In more recent years, no negative anomalies of that strength have been recorded.

Figure 3 VOD_MW and LAI_MD anomalies over the Okavango catchment (a, b, e), with the intensity of the colouring based on the z-score of the positive (green) and negative (brown) anomalies and the Okavango delta (c, d, f) in time series and density plots. A daily time step is used from the moving average dataset.

4.3 Land surface temperature

Figure 4a and b show the anomalies of LST_MW over the Okavango catchment and delta, with Fig. 4b–g showing a direct comparison in a density plot between LST_MW, LST_E5, and LST_MD (time series of LST_E5 and LST_MD anomalies can be found in the Appendix, Fig. A1). Because of the high correlation between LST_MW and LST_E5 of 0.88 in the catchment and 0.81 in the delta, the decision was made to only show the LST_MW time series to focus more on the density plots of the three different products. The correlation of LST_MW against LST_MD is much lower, with 0.64 and 0.26 for both regions, showing a low relation in the catchment. LST_E5 compares better to LST_MD with a correlation of 0.70 in the catchment and 0.65 in the delta, however, this is still significantly lower than the comparison with LST_MW. The absolute ranges in the anomalies as detected by the three products are very similar.

Figure 4 LST_MW time series over the Okavango catchment and the Okavango delta (a, b), with the intensity of the colouring based on the z-score of the positive (red) and negative (blue) anomalies. For the density plots; LST_MW compared to LST_E5 (c, f), LST_MW compared to LST_MD (d, g), and LST_E5 compared to LST_MD (e, h). A daily time step is used from the moving average dataset for the density plots.

The slightly lower correlation of LST_MW against LST_E5 in the delta is mostly caused by the period 2010 and 2011, when the LST_E5 anomaly (between −1 and 1 ^∘C) is smaller than that of LST_MW (between −3 and −1 ^∘C). This breakaway is clearly visible in the density plot of Fig. 4f on the lower left side. Below-average temperatures are recorded for a prolonged period between 2006 and 2014 in both regions. For the delta, the highest temperature anomalies are recorded in 2019 and 1998. In the catchment, this is seen in 2015 and 2019.

4.4 River and precipitation contribution to soil moisture anomalies in the Okavango delta

Figure 5a–c shows the anomalies of the ORD, ODIA_MD, and PR_E5 over the delta, which have visibly different signals. The ORD shows a strong multi-year signal with especially high values recorded from 2009 to 2012. Outside of that period, with the exception of 2004, values generally lie below the 2003 to 2020 climatology. The OIAD_MD shows a signal that is relatively similar to that of the ORD, however, smoother with less variability and lagging behind. The PR_E5 over the delta shows mostly values around 0 mm during this 2009 to 2012 period, and otherwise varies more dynamically from year to year with values above and below the climatology. Additional PR_E5 and PR_IM anomaly time series over both the catchment and the delta can be found in the Appendix (Fig. A2).

Figure 5 ORD, ODIA_MD, and PR_E5 time series over the Okavango delta (a–c), with the intensity of the colouring based on the z-score of the positive (blue) and negative (red) anomalies.

Although SSM_MW and SSM_E5 anomalies have an overall correlation of 0.69 in the delta, Fig. 2f shows many occasions where the SSM_MW had negative or positive anomalies, while the SSM_E5 did not diverge from the climatology. To better assess what causes this opposite signal, the climatology (using ±15 d moving average) of different parameters are provided in Fig. 6, including their 10 % and 90 % percentiles. Here, one can see the difference in the dynamics between SSM_MW (Fig. 6a) and SSM_E5 (Fig. 6b). The SSM_E5 shows a clear relation to the PR datasets (Fig. 6g and h), while the SSM_MW still picks up a moisture signal between April and September. When looking at the ORD and OIAD_MD, these are the moisture-related signals that still show strong variability in this time of the year, indicating that the SSM_MW could also contain information from other sources than PR. On a side note, Fig. 6 shows that besides matching well with long-term anomalies, LST_MW and VOD_MW also have a strong matching intraseasonal signal with LST_E5 and LAI_MD, respectively.

Figure 6 ROI₂ mean climatology (coloured thick line) and both 10 % and 90 % percentiles (black dashed lines) for SSM_MW (a), SSM_E5 (b), LST_MW (c), LST_E5 (d), VOD_MW (e), LAI_MD (f), PR_IM (g), PR_E5 (h), ODIA_MD (i), and ORD (j). Data using a ±15 d moving average were plotted to distinguish between intraseasonal signals.

Table 2 presents the results of a multiple linear regression to determine the drivers of the observed/modelled SSM anomaly signal in the delta, using ODIA_MD, ORD, and PR as inputs. The z-score anomalies are used to improve the comparability between the different datasets and their weight. The results show that the weighting for SSM_MW consists of a balance between the PR in the delta and the ODIA_MD, with an overall slightly higher weight for the ODIA_MD, and leading to a maximum correlation of 0.84 when using PR_E5 over PR_IM. This leads to a RMSE of about 0.44 for the z-score. The SSM_E5 anomalies are clearly driven by the PR_E5 anomalies, reaching a correlation of 0.87. The correlation strongly decreases to 0.64 when the PR_E5 is replaced with PR_IM, which reflects back in the RMSE of the z-score anomalies, which increases from about 0.37 to 0.57.

5 Discussion

Over both the delta and the catchment, a remarkably strong relationship between the LAI_MD and VOD_MW was observed, even though fundamentally, they measure two different things. The relationship is slightly weaker over the catchment where you see more of a buffer effect in the VOD_MW dataset as compared to the LAI_MD. This could be caused by a buildup of woody biomass as this would theoretically be better detected with the VOD_MW than with the LAI_MD. The period of sustained high VOD_MW in the catchment during the 2008 to 2012 period aligns well with the PR_E5, which recorded 5 years of above-average rainfall over the catchment. The ORD shows this increase above the climatology starting only the year afterwards (from 2009 to 2012), showing the lagged response of the system after a prolonged drier period.

The VOD_MW signal in the delta is more complex: the peaks in VOD_MW do not coincide with prolonged time spans of high water availability, but seem to peak during shorter periods of increased water availability during overall conditions with medium to low ODIA_MD. This can be explained by what the VOD represents: in this case, it is related to biomass that is above the surface. During prolonged periods of high water, a larger extent of these regions are flooded. Therefore, within the 0.25^∘ pixel, data that are not corrected for dynamic water bodies, the vegetation covered by these flooded areas might not be properly measured by the VOD_MW signal as it is also known that VOD_MW values can be underestimated during flooded conditions (Bousquet et al., 2021). Note that the negative SSM_MW and ORD anomalies in 2019 have not led to the same intensity of vegetation decline, while in 2019, the ODIA_MD was at a record low in the last 20 years. With the very strong relationship over the delta between the anomalies of both VOD_MW and LAI_MD – two independent satellite-observed datasets – these observations very likely reflect the conditions on the ground. These results show that future use of even longer VOD_MW records can help monitor complex regions like the Okavango catchment and delta. For example, following the progress on VODCA – which aims to build a data record similar to the ESA CCI SSM for VOD_MW – future releases will also include the latest calibrated datasets as used here.

Three different sources of LST were tested over the Okavango catchment and delta. The highest correlation can be found between the LST_E5 and LST_MW, which most likely best represent the actual ground conditions. Although LST_MD performs less well, the better correlation of LST_MD against LST_E5 than LST_MW might indicate that the overall best performing dataset is the LST_E5. However, in many cases, an observation-based long-term dataset (e.g. the LST_MW) is still preferred. For example, in 2010 and 2011, the LST_MW has the lowest temperature anomalies on record in the delta, going to −3 K, while the LST_E5 remains more neutral. This is most likely caused by the lack of lateral water influx modelling from the ORD and following ODIA_MD in ERA5-Land (Muñoz-Sabater et al., 2021), as shown in Sect. 4.4. The lack of moisture input into the model can lead to an underestimation of the latent heat flux and overestimation of the sensible heat flux, leading to an unrealistically high LST_E5.

In the delta, 2015, 2016, and 2019 have been warm compared to the years before. The LST_MW and LST_E5 both show that these are not unique occurrences as similar high values have been detected on multiple occasions before 2006. These seem to occur during periods of lower ODIA_MD, which shows dry anomalies of varying strength in these years. The catchment does see its highest and more prolonged peaks only in the last years, i.e. 2015 and 2019. These high peaks coincide with the strongest negative anomalies found for both SSM_MW and SSM_E5, linking the high temperature and reduced moisture availability.

The precipitation-driven SSM in the catchment aligns closely with SSM_MW and SSM_E5 datasets. Especially in the period after 2010, the signal in the anomalies is very similar. Before 2010, it appears that the SSM_MW shows slightly stronger dynamics than SSM_E5. In the delta, a mismatch is clearly seen between SSM_E5 and SSM_MW, especially with regard to the duration of the dry and wet peaks, but also in their intensity. With the knowledge that about 50 %–75 % of the water flux into the delta comes from the ORD, and about 25 %–50 % from the PR, an analysis using z-score anomalies was conducted to determine the driving signals behind the SSM anomalies, using the ORD, ODIA_MD, and PR as inputs. For SSM_E5, an almost one-to-one relationship was found with the PR, with little to no effects from the ORD or ODIA_MD. The SSM_MW anomalies, on the other hand, are almost equally driven by PR and ODIA_MD, which is much closer to the actual balance between the ORD and PR water fluxes for the delta, as expected from literature.

The almost one-to-one relationship between the SSM_E5 and PR_E5, and lack of signal related to the ORD due to the missing lateral water influx modelling, or alternatively dynamic open water bodies using the ODIA_MD, in ERA5-Land, indicates that in a complex region like the Okavango delta, important forcings are missing. This, for example, could also cause the difference in LST_MW and LST_E5 in 2010 and 2011 (not shown) as the model cannot correctly convert the incoming radiation into sensible and latent heat fluxes when the moisture conditions are inaccurate, leading to a false increase of LST. On the other hand, while the SSM_MW signal provides users with a better representation of total moisture conditions within the catchment, it can also not be interpreted as a pure SSM signal here as it includes moisture information driven by the ODIA_MD. In a dynamic environment such as the Okavango delta, users should therefore clearly define what they require of such datasets to avoid unwanted side effects.

6 Conclusion

The anomalies of three different parameters, i.e. SSM_MW, LST_MW, and VOD_MW, were evaluated against other satellite-observed datasets and data from the ERA5-Land climate reanalysis. Although SSM_MW and SSM_E5 correlate moderately well, structural differences were detected over the Okavango delta, where SSM_MW contains a clear multi-year signal that is not in the SSM_E5. To determine the cause of this mismatch, an analysis was conducted to determine the impact of three sources of water into the Okavango delta, i.e. the ORD, ODIA_MD, and the PR, on the SSM signal. The SSM_MW signal appears to be driven about equally by the ODIA_MD and the PR, while SSM_E5 is almost fully driven by the PR_E5. This indicates that ERA5-Land does not properly include the lateral influx of the Okavango River and, therefore, the use of SSM_MW is preferred in this region.

For the VOD_MW, a direct comparison against LAI_MD was made. Although the two parameters measure two different physical characteristics of the vegetation, their anomalies show a similar response, which were reflected back in their strong correlations (0.74 and 0.88). Over the catchment, a stronger multi-year signal was detected in the VOD_MW, which could be related to the buildup of biomass, to which VOD_MW is theoretically more sensitive. For the delta, both datasets are impacted by the increase in open water during long wet periods that can suppress the observed vegetation. This strong similarity, as observed between the two datasets, indicates that it is very likely they are both representative for the in situ conditions.

LST_MW was shown to be of good quality and correlated well with LST_E5 (>0.8). LST_MD still managed to reach a significant correlation with LST_E5, but not with LST_MW, indicating that in general LST_E5 could be of the highest quality of the three when looking at the temporal signal. However, at the record low values in LST_MW over the delta in 2010–2011, corresponding to the peak years of the ORD and ODIA_MD, it seems that LST_E5 cannot properly model the sensible and latent heat fluxes because it is missing the lateral water component. This can have a large impact on detecting extremes, which are especially important in the current changing climate.

The findings of this research show the importance of not only relying on climate reanalysis, but also the need for further development and maintenance of observational datasets like the ones derived from passive microwave observations. For example, within the ESA CCI soil moisture datasets, but also the development of new climate data records (CDRs) on VOD_MW like VODCA. Their ability to properly detect anomalies and extremes is very valuable in climate research and can especially help to improve our insight in complex regions where the current climate reanalysis datasets reach their limitations. With microwave data being available from 1978 onwards, the data can be used for long-term climate studies, near-real-time applications, e.g. monitoring complex natural systems like the Okavango delta, and to constrain climate reanalysis through data assimilation techniques to overcome known model weaknesses.

Appendix A

Figure A1 LST_E5 (a, b) and LST_IM (c, d) time series over the Okavango catchment and the Okavango delta, with the intensity of the colouring based on the z-score of the positive (red) and negative (blue) anomalies.

Figure A2 Time series of PR_E5 over the Okavango delta (a) and PR_IM over both the Okavango catchment and delta (b, c), with the intensity of the colouring based on the z-score of the positive (blue) and negative (red) anomalies.