Evaluation of an Adaptive Soil Moisture Bias Correction Approach in the ECMWF Land Data Assimilation System

1. Introduction

Inherent biases exist when assimilating satellite-derived soil moisture (SM) observations in land surface models. Observation biases may result from instrument errors, vegetation effects and errors of representativeness whilst model biases can originate from the model physics, parameterizations, initial conditions and atmospheric forcing [1,2]. Biases can be addressed by data assimilation when it is not possible to correct the source directly.

Commonly in atmospheric and ocean DA, a variational bias-correction (BC) approach is employed for satellite radiances, whereby the corrections for the computed forward radiances are updated within the cost function minimization [3,4]. Additionally, a network of anchor observations are employed as a reference (e.g. radiosondes) to prevent the analysis from drifting to the model bias [5]. Land surface models are much more heterogeneous than atmospheric or ocean models and pointwise reference soil moisture measurements are not generally representative of the large-scale footprint associated with model or satellite-derived products [6]. Without an accurate ground reference, Numerical Weather Prediction (NWP) centres typically assume that the land surface model is perfect and all the biases belong to the observations. Traditionally a cumulative distribution function (CDF) matching approach has been used to rescale the mean and variance of the observations to match the model climatology [7,8,9,10]. Machine learning methods have been advocated in recent years, including neural networks that convert level 1 satellite data to level 2 soil moisture [11,12]. Whilst these methods work well in many practical applications, large calibration/training datasets are required and the assumption of stationary biases is often suboptimal. Biases in the observations can change over time due to issues with the instruments and modifications in the retrieval algorithms, whilst land surface models are also updated periodically. Both these factors change the observation-model climatology and necessitate a recalibration of the climatological bias-correction.

At the European Centre for Medium Range Weather Forecasts (ECMWF), Advanced Scatterometer (ASCAT) and Soil Moisture Ocean Salinity (SMOS) derived soil moisture observations are both assimilated operationally in the SM analysis using a simplified extended Kalman filter (SEKF) [13,14]. The ASCAT SM observations are bias-corrected via the aforementioned rescaling approach. On the other hand, SMOS SM is derived from the level 1 product using a neural network (NN) trained on the ECMWF operational analysis [11]. In this study an adaptive SM bias-correction is tested for ASCAT and SMOS in order to capture non-stationary biases. The two-stage adaptive BC was introduced by Dee and Silva [15], whereby the filter and the bias-correction updates are performed independently. Applications can be found for the atmosphere [1], ocean [16,17] and land [18,19,20,21]. In Draper et al. [22] the bias-correction is applied to skin-temperature assimilation in the catchment and surface modelling component of the Goddard Earth Observing System. A similar approach is employed here, in the context of satellite-derived SM assimilation in the ECMWF integrated forecasting system (IFS) for NWP applications.

The manuscript is structured as follows. Firstly, the materials and methods are presented in Section 2, including the current SEKF formulation and the two-stage bias filter. The results are given in Section 3. Internal DA diagnostics are evaluated in terms of the first guess (observation-model) departures and the SM increments. Then a validation of the SM and soil temperature (ST) analysis performance is carried out using in situ data from the International Soil Moisture Network (ISMN). The validation of the subsequent NWP forecasts is performed using the ECMWF high resolution analysis as a reference. Finally, Section 4 discusses the key findings and recommends future work.

2. Materials and Methods

2.1. The bias-free SEKF soil moisture analysis

A detailed documentation of the ECMWF land (ECland) data assimilation for IFS cycle 47r3 can be found in ECMWF [23]. The operational SM analysis is based on a point-wise SEKF with 12-hour assimilation windows, which was implemented in 2010 to assimilate proxy screen-level observations of 2m temperature (T2m) and relative humidity (RH2m) [13]. These proxy observations come from the analysed T2m and RH2m states at the synoptic times (every 6 hours), which are derived by assimilating SYNOP screen-level observations with a 2D optimal interpolation scheme. ASCAT-derived SM observations have been assimilated in operations since 2015 followed by SMOS NN assimilation since 2019 [11]. Up to 8 observations can be assimilated for each grid point per assimilation window, with a maximum of 2 each for ASCAT, SMOS, T2m and RH2m. The top 3 layers (top metre) of the ECland model are analysed with depths of 0-7 cm, 7-28 cm and 28-100 cm from top to bottom. Using the notation of Ide et al. [24], the pointwise SM analysis state update at time ($t_i$) is expressed as:

$${\mathbf{x}}^a\left(t_i\right)={\mathbf{x}}^b\left(t_i\right)+{\mathbf{K}}_i[{\mathbf{y}}^o\left(t_i\right)-H_i\left({\mathbf{x}}^b\right)],$$

(1)

where the superscripts a, b, o and subscript i denote the analysis (of dimension n), background (of dimension n), observations (of dimension p) and timestep respectively. The observation operator H maps the model state to observation space. As the preprocessed ASCAT and SMOS NN observations are already in the modelled volumetric SM units, $H_i\left({\mathbf{x}}^b\right)$ is simply the model value in the top layer at the nearest gridpoint and timestep to the observation. The weights of the observations and background state in the analysis update are determined by the Kalman gain matrix K (of dimension $n\times p$):

$${\mathbf{K}}_i={[{\mathbf{B}}^{-1}+{\mathbf{H}}_i^T{\mathbf{R}}^{-1}{\mathbf{H}}_i]}^{-1}{\mathbf{H}}_i^T{\mathbf{R}}^{-1},$$

(2)

where $\mathbf{R}$ (of dimension $p\times p$) and $\mathbf{B}$ (of dimension $n\times n$) represent the observation and background error-covariance matrices respectively, which are assumed to be static and diagonal (uncorrelated). The SM background errors are prescribed with values ${\sigma }_b=0.01{\mathrm{m}}^3/{\mathrm{m}}^3$. ASCAT and SMOS observations are prescribed values of ${\sigma }_{ascat}=0.05{\mathrm{m}}^3/{\mathrm{m}}^3$ and ${\sigma }_{s\mathrm{m}os}=0.02+ϵ$, where $ϵ$ is a situation dependent uncertainty output provided by the NN itself. The screen-level variable (SLV) observations are assigned values of ${\sigma }_{T2\mathrm{m}}=1K$ and ${\sigma }_{RH2\mathrm{m}}=4\%$ [23]. Flow-dependent uncertainty information from an ensemble of data assimilations (EDA, Bonavita et al. [25]) is implicitly propagated from the observations to the analysis state via the Jacobian matrix $\mathbf{H}$ (of dimension $p\times n$). The Jacobian linking the $k^{th}$ observation to modelled soil moisture layer j is given by:

$${\mathbf{H}}_{k,j}=\frac{cov(H_k\left({\mathbf{x}}^{eda}\right),{\mathbf{x}}_j^{eda})}{var\left({\mathbf{x}}_j^{eda}\right)}.c_j,$$

(3)

where “cov" stands for “covariance", “var" stands for "variance" and “eda" implies that the EDA ensemble is used. Tapering coefficients $c_j=1/(1+(j-1).{\alpha }_{sekf})$ were empirically derived to optimise the impact for the different SM layers ($j=1,2,3$), with ${\alpha }_{sekf}=0.6$. In the following assimilation window, the background (prior) state is a model simulation initialised by the current analysis state:

$${\mathbf{x}}^b\left(t_{i+1}\right)=\mathrm{m}\left[{\mathbf{x}}^a\left(t_i\right)\right],$$

(4)

where M is the nonlinear coupled forecast model.

2.2. The two-stage bias filter

In the bias-aware SM analysis, the biased observations $\left({\tilde{\mathbf{y}}}^o\right)$ are partitioned into the analysed bias ( ${\mathbf{z}}^a\left(t_i\right)$) and the non-biased term (${\mathbf{y}}^o\left(t_i\right)$):

$${\tilde{\mathbf{y}}}^o\left(t_i\right)={\mathbf{z}}^a\left(t_i\right)+{\mathbf{y}}^o\left(t_i\right).$$

(5)

The state update then includes the bias-corrected observations:

$${\mathbf{x}}^a\left(t_i\right)={\mathbf{x}}^b\left(t_i\right)+{\mathbf{K}}_i[{\tilde{\mathbf{y}}}^o\left(t_i\right)-{\mathbf{z}}^a\left(t_i\right)-H_i\left({\mathbf{x}}^b\right)].$$

(6)

In this study we only investigate biases in the ASCAT and SMOS NN observations and therefore we assume that the proxy SLV observations are bias-free. For simplicity, the ASCAT and SMOS NN bias updates are performed separately over each gridpoint i.e. they are assumed to be uncorrelated. The bias update is calculated following the approach of Draper et al. [22], as:

$${\mathbf{z}}_l^a\left(t_i\right)={\mathbf{z}}_l^b\left(t_i\right)+{\mathbf{L}}_{i,l}[{\tilde{\mathbf{y}}}_l^o\left(t_i\right)-{\mathbf{z}}_l^b-H_{i,l}\left({\mathbf{x}}^b\right)],$$

(7)

$${\mathbf{L}}_{i,l}=\left[{\mathbf{B}}_l^{\mathbf{z}}\right]{[{\mathbf{R}}_l+{\mathbf{B}}_l^{\mathbf{z}}+{\mathbf{H}}_{i,l}\mathbf{B}{\mathbf{H}}_{i,l}^T]}^{-1},$$

(8)

where the subscript l indicates the observation type and ${\mathbf{B}}_l^{\mathbf{z}}$ is the prior observation bias covariance matrix. The matrix ${\mathbf{R}}_l$ represents the random part of the observation-error covariance matrix using the bias-free values of $\mathbf{R}$ defined for ASCAT and SMOS NN in equation (2). The Jacobian matrix ${\mathbf{H}}_{i,l}$ is the subset of the Jacobian matrix ${\mathbf{H}}_i$. Following [15,19], ${\mathbf{B}}_l^{\mathbf{z}}$ is a diagonal matrix, which is proportional to the SM background-error covariance matrix:

$${\mathbf{B}}_l^{\mathbf{z}}=\frac{\gamma }{1-\gamma }{\mathbf{H}}_{i,l}\mathbf{B}{\mathbf{H}}_{i,l}^T,$$

(9)

where $\gamma$ is a scalar parameter. After putting (9) in (8) and re-arranging, the Kalman gain can be expressed as:

$${\mathbf{L}}_{i,l}=\gamma {\mathbf{H}}_{i,l}\mathbf{B}{\mathbf{H}}_{i,l}^T{[{\mathbf{H}}_{i,l}\mathbf{B}{\mathbf{H}}_{i,l}^T+(1-\gamma ){\mathbf{R}}_l]}^{-1}.$$

(10)

Increasing the value of $\gamma$ between 0 and 1 effectively reduces the memory of the bias. In this study the value of $\gamma =0.25$ was chosen empirically in order to have a long enough memory to capture seasonal-scale changes in the bias. After concatenating the updated bias correction for each observation type from equation (7) into a single vector ${\mathbf{z}}^a\left(t_i\right)$ (of dimension p), the SM analysis is calculated with equation (6). Recall that the screen-level observations are assumed to be bias-free and therefore the elements of ${\mathbf{z}}^a\left(t_i\right)$ are set to zero for these observations. Over gridpoints where no observations are available per type l, the bias state is not updated. If only one observation is available, equation (7) is a scalar update. In cases where 2 observations of the same type are available for assimilation, ${\mathbf{z}}_l^a\left(t_i\right)$ is a vector of size 2. In the SM analysis the observations are assumed to be representative of the model gridpoint, hence a weighted average of the bias update is stored and used for the following cycle:

$${\tilde{z}}_l^a\left(t_i\right)=\frac{1}{K}\sum _{k=1}^K{z}_{l,k}^a\left(t_i\right),$$

(11)

where k is the observation index and K is the total number of observations of type l in the gridpoint (maximum of 2 for ASCAT or SMOS NN). The persistence model is employed for the bias forecast, which is a reasonable assumption if the bias evolves slowly between cycles [15,19]. Hence the background bias state for the following cycle is defined as:

$${\tilde{z}}_l^b\left(t_{i+1}\right)={\tilde{z}}_l^a\left(t_i\right).$$

(12)

The background bias states for each gridpoint are then used for all the observations of type l in the gridpoint for the following cycle:

$${\mathbf{z}}_l^b\left(t_{i+1}\right)={\left({\mathbf{1}}_K\right)}^T{\tilde{z}}_l^b\left(t_{i+1}\right),$$

(13)

where ${\left({\mathbf{1}}_K\right)}^T$ is a vector of ones of dimension K.

2.3. Satellite-derived soil moisture observations and pre-processing

An introduction to the ASCAT and SMOS derived SM data is given below along with the current SM bias correction techniques.

Active C-band scatterometer data is provided by the ASCAT sensors on board the Metop satellites. These observations are converted into a liquid soil wetness percentage using the change detection approach [26,27] with an approximate spatial resolution of 25 km and a revisit time of about 1-3 days. During the study period (2019-2022), ASCAT from Metop-B was available throughout. ASCAT Metop-C assimilation was introduced in October 2019 and Metop-A assimilation was retired in November 2021. ASCAT SM is measured as a soil wetness percentage (between 0 and $100\%$) and the modelled SM is expressed in volumetric soil water content. Cumulative Distribution Function (CDF) matching bias correction rescales the ASCAT level 2 SM such that the observed CDF matches the modelled SM CDF. Systematic differences between the observations and the model are subsequently removed and the observations are effectively converted into volumetric units. At ECMWF and various other centres, only the first two moments of the observation CDF, the mean and variance, are rescaled to match the model mean and variance [9]. A CDF matching which accounts for seasonal variability was implemented by [28,29] and has been adopted by ECMWF in the Simplified Extended Kalman filter (SEKF) SM analysis [13] as well as for SMOS brightness temperature bias correction [10]. Hereafter this method will be called “seasonal rescaling". The seasonal rescaling parameters are derived for each calendar month using a 3-month moving average. For example, for the month of May the parameters are derived based on the data from April, May and June over a multi-year calibration period. The operational seasonal rescaling parameters for IFS cycle 47r3 are based on rescaling ASCAT-A/B derived SM to the ERA5 surface soil moisture (SSM) over the period 2009-2018 [30]. The ASCAT observations undergo a quality control check prior to assimilation, which screens observations during frozen conditions, over mountainous regions and where the estimated noise exceeds $15\%$.

The SMOS SM observations are derived from L-band brightness temperature (Tb) using a neural network approach, which is summarised below but full details can be found in Rodríguez-Fernández et al. [11,31]. They are trained and validated using the high resolution ECMWF operational SM analysis over the period 2015-2020. Firstly the SMOS Tb observations are colocated with the model gridpoints in space and time. In order to maximise the data availability, Tb observations are extracted from incidence angles ranging from 30 to ${45}^{\circ }$. Twelve predictors are used, consisting of 6 Tb and 6 SSM linear expectations and three angular bins. Furthermore, the ECMWF operational snow depth and temperature analysis are used to filter out data during frozen conditions. After filtering, 60% of the samples are used for training, 20% are used for an evaluation of the training dataset and a further $20\%$ for the validation of the SMOS NN derived SM. A gradient back-propagation approach is employed for the training using the Levenberg-Marquardt algorithm [32]. The architecture consists of a 2 layer NN with one hidden layer and 5 neurons. The NN is trained globally which means that the global biases between the SMOS SM observations and the model SM are small but significant regional biases remain.

2.4. ECland surface model

The ECLand surface model is used in the IFS [33]. It is inherited from the Hydrology Tiled ECMWF Scheme for Surface Exchanges over Land (HTESSEL; [34]). ECland models land surface processes including soil moisture, snow and vegetation. Each gridpoint is divided into 8 tiles which represent the fraction of different land cover types (including vegetation types, soil and snow cover) using data from the US Department of Agriculture Global Land Cover Climatology (GLCC) map [35]. The seasonal vegetation cycle is constructed from the monthly leaf area index climatology [36]. Soil moisture is represented by 4 vertical layers with thicknesses of 7 cm, 21 cm, 72 cm and 189 cm from top to bottom. The vertical SM exchanges are based on the equation of Richards [37] and a multi-layer snow scheme controls the snowpack evolution [38].

2.5. The stand-alone surface analysis

This study makes use of the stand-alone surface analysis (SSA) of Fairbairn et al. [30]. In each cycle, the atmospheric initial conditions for the coupled model are forced by an external atmospheric analysis, which avoids the computational burden of running the 4D-Var atmospheric analysis. Atmospheric forcing consists of the the ERA5 reanalysis fields [39]. However, land-atmosphere feedback is provided by the coupled forecast model between cycles. Fairbairn et al. [30] demonstrated that SSA can approximate the impact of changes to the soil moisture analysis on NWP forecasts in the IFS with a small fraction ($<30\%$) of the computing time compared with the operational weakly coupled DA approach. SSA is particularly useful for running multi-year experiments in a reasonable time frame, which is necessary to monitor the long timescales associated with SM bias correction.

Other than the atmospheric analysis, SSA benefits from most of the functionalities in the IFS. The land surface analysis, as described by ECMWF [23], includes a 2D OI for 2 metre temperature and relative humidity and the snow depth analysis. Additionally, soil temperature is analysed with a 1D optimal interpolation (OI). As explained by Haseler [40], two 12-h assimilation windows run from 2100 to 0900 and 0900 to 2100 UTC, which provide the initial conditions for the 10-day coupled NWP forecasts at 00 UTC and 12 UTC respectively.

2.6. Experiments

Five SSA experiments were validated globally over a 3-yr period (1 January 2020-31 December 2022) after 1-yr spinup (1 January 2019-31 December 2019). All the experiments assimilate both ASCAT and SMOS derived SM observations. Table 1 summarises the differences between the experiments: (i) C is the control with the operational configuration of ASCAT and SMOS assimilation with no adaptive BC, (ii) ${\mathrm{E}}_{\mathrm{A}}$ includes ASCAT adaptive BC, (iii) ${\mathrm{E}}_{\mathrm{S}}$ includes SMOS adaptive BC and (iv) ${\mathrm{E}}_{\mathrm{A},\mathrm{S}}$ includes both ASCAT and SMOS adaptive BC. In experiments (i) to (iv) ASCAT seasonal rescaling is applied before the adaptive bias-correction. Experiment (v) (${\mathrm{E}}_{\mathrm{A},\mathrm{S}}*$) is designed to test whether the adaptive BC alone can improve on the seasonal rescaling for ASCAT. In experiment (v) the level 2 ASCAT SM for each grid point was converted to volumetric units by rescaling the percentage to SWI (between 0 and 1) and multiplying the SWI by the modelled saturation value. Each experiment was initialized by the ERA5 reanalysis at 00 UTC on 1 January 2019. In subsequent cycles, the atmospheric analysis was initialized by ERA5 (as described in Section 2.5). The land surface analysis and the coupled model were performed with cycle 47r3 of the ECMWF IFS. The experiments were run on a cubic octahedral reduced Gaussian grid at Tco319 (approximately 31 km ). This resolution is very close to the ERA5 native resolution, avoiding issues with spatial interpolation of the atmospheric fields.

2.7. SM and ST validation approach

The surface (0-7 cm) and root-zone (0-100 cm) layers of SM and ST were validated using sparse in situ data from the international soil moisture network [41]. Observations were extracted at the nearest hour to the analysis times (00/12 UTC ±30 min) from 7 networks over the United States, Europe and Australia. The North American networks consist of the U.S. Climate Reference Network (USCRN, [42]) and the U.S. Department of Agriculture’s Soil Climate Analysis Network (SCAN)/Snowpack Telemetry (SNOTEL) networks (Schaefer et al. 2007). In Europe, they consist of the Soil Moisture Observing System–Meteorological Automatic Network Integrated Application (SMOSMANIA, [43] in southwestern France, the Soil Moisture Measurement Stations Network of the University of Salamanca (REMEDHUS) in the central sector of the River Duero basin in Spain [44], and the Network of Terrestrial Environmental Observatories in Germany (TERENO, [45]). In Australia, the Murrumbidgee soil moisture monitoring network data set was used (OZNET, [46]). The depths, locations and number of stations in each network are listed in Table 2. Although the validation of model grid point data using in situ observations is affected by representativeness and instrument errors, the locations of these networks represent a wide range of vegetation and soil types. Similar validations have been performed at ECMWF [30,47] and elsewhere [48,49] for validating global SM analyses and satellite products. The in situ SM observations are compared with the root-zone soil moisture (RZSM) over the top metre of soil by taking a vertical average, with weights that are proportional to the spacings between the sensor depths (as in [49]). For SNOTEL, SMOSMANIA, OZNET and TERENO, the deepest measurement is assumed to represent the depth of the observation down to 1 m. In each network the top measurement at 5 cm depth is compared to the surface layer of the soil moisture analysis (0-7 cm depth). The observations undergo a rigorous quality control check by the ISMN facility [50] with flags during frozen conditions and implausible SM values (e.g. spikes). Additional quality control screening is implemented in this study where the analysed soil temperature is below 4°C, which reduces the risk of frozen conditions. Further, a minimum of 50 observations is required for each station during the validation period to reduce sampling errors.

The validation metrics consist of the Pearson R anomaly correlation coefficient (${\mathrm{R}}_{\mathrm{ano}}$), unbiased RMSE (UbRMSE) and the bias, as described by Fairbairn et al. [30]. These metrics offer a complementary validation of the SM performance [51]. The correlation coefficient R is affected by the seasonal SM cycle. For ${\mathrm{R}}_{\mathrm{ano}}$, these seasonal-scale correlations are removed by subtracting the centred 35-day moving average from the time series. As in Albergel et al. [47], the p-value provides a measure of the significance of the correlations and the scores for each station are only retained where the p-value if less than 0.05. Further, a Fisher-z transform with a lag-1 auto-correlation was implemented to find the 95% confidence interval (CI) while mitigating the autocorrelations associated with the seasonal cycle [48].

2.8. Atmospheric validation approach

In line with the ECMWF NWP setup, 10-day coupled forecasts are initialized at 0000 and 1200 UTC for all the experiments. The upper air verification is performed by comparing the forecasts against the ECMWF operational analysis. Following Geer [52], the normalized RMSE differences (dRMSE) are used to validate the relative performance of the experiment forecasts compared to the control forecasts:

$$\mathrm{d}\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}=\frac{∥{\mathbf{z}}_{t:t+T}^e-{\mathbf{z}}_{t+T}^r∥-∥{\mathbf{z}}_{t:t+T}^c-{\mathbf{z}}_{t+T}^r∥}{∥{\mathbf{z}}_{t:t+T}^c-{\mathbf{z}}_{t+T}^r∥},$$

(14)

where ${\mathbf{z}}_{t:t+T}^c$ (${\mathbf{z}}_{t:t+T}^e$) are the control (experiment) forecasts of length T from analysis time t and ${\mathbf{z}}_{t+T}^r$ is the reference at time $t+T$. The $∥.∥$ in this case is the RMSE of the time series, which is calculated independently for each grid point. The reference is provided by the ECMWF high resolution deterministic analysis. Further, the difference in mean error between the experiments and the control is computed:

$$\mathrm{d}\mathrm{M}\mathrm{E}=\overline{{\mathbf{z}}_{t:t+T}^e-{\mathbf{z}}_{t+T}^r}-\overline{{\mathbf{z}}_{t:t+T}^c-{\mathbf{z}}_{t+T}^r}.$$

(15)

3. Results

3.1. Internal DA diagnostics

Firstly the impact of the adaptive BC on the ASCAT SM first guess departures is evaluated. Recall that the departures are calculated as $({\mathbf{y}}^o-H\left({\mathbf{x}}^b\right))$ in the control and $({\tilde{\mathbf{y}}}^o-{\mathbf{z}}^a-H\left({\mathbf{x}}^b\right))$ in the adaptive BC experiments. Figure 1 (a) shows the global mean monthly ASCAT SM departures over 2019-2021 for the control ($\mathrm{C}$) and experiment $\left(E_{A,S}\right)$. Whilst the mean departure oscillates close to zero in the first few months, it then trends upwards in $\mathrm{C}$. At the same time the BC $(-{\mathbf{z}}^a)$ applied to ${\mathrm{E}}_{\mathrm{A},\mathrm{S}}$ trends downwards, which results in less biased departures. The spin-up for the bias-correction takes about a year, after which the absolute departures (Figure 1 (b)) are about $10-20\%$ larger in $\mathrm{C}$ compared to ${\mathrm{E}}_{\mathrm{A},\mathrm{S}}$, thus confirming that the bias is a substantial part of the departure magnitude. Figure 1 (c) displays the used ASCAT SM observation count. Up to $6\%$ more observations are assimilated in ${\mathrm{E}}_{\mathrm{A},\mathrm{S}}$ compared to $\mathrm{C}$, with the biggest gap seen during the northern hemisphere summer/autumn, when ASCAT assimilation is most active in the northern hemisphere. Further investigation revealed that the quality control was rejecting a greater number of observations in $\mathrm{C}$ compared to ${\mathrm{E}}_{\mathrm{A},\mathrm{S}}$ when the departure size exceeded the maximum threshold of 0.1m³/m³ (not shown). Note that the bias-correction update is still applied when the departure size exceeds this threshold, which enables systematic large departures to be corrected. Figure 1(c) also demonstrates a substantial increase in the peak observation count after the introduction of ASCAT-C and a similarly large reduction in the peak count following the retirement of ASCAT-A, which highlights the non-stationary nature of satellite observations.

Figure 1 (d), (e) and (f) demonstrate the impact of the adaptive BC on the SMOS SM departures. In Figure 1(d) the mean global bias in $\mathrm{C}$ is highly variable, ranging from slightly negative to highly positive. The experiment ${\mathrm{E}}_{\mathrm{A},\mathrm{S}}$ responds with a negative bias-correction, although it sometimes over-corrects in the negative direction. Nevertheless, the absolute departures are about $20-30\%$ smaller in ${\mathrm{E}}_{\mathrm{A},\mathrm{S}}$ than $\mathrm{C}$ (Figure 1(e)), which indicates that the BC is mainly effective. Overall, up to $12\%$ more SMOS SM observations are assimilated when adaptive BC is turned on (Figure 1(e)), with the biggest gap again seen in the northern hemisphere summer.

Next the average ASCAT and SMOS SM departures are evaluated over August 2022 in Figure 2. In the control, ASCAT has mainly positive departures over eastern Europe. On the other hand, SMOS SM departures are generally positive over South America and South Africa. Large departures are seen over some areas that are known to have quality control issues, such as the Andes mountain range in South America for SMOS SM and across the high latitudes for both SMOS and ASCAT. Evidently the bias-correction reduces the departure magnitudes in these areas. Under closer inspection it is evident that the bias-correction is more active over some regions than others. For example, the strong negative SMOS SM departures over India and the Sahel region of Africa are barely corrected by the adaptive BC. Further analysis has shown that the negative departure patterns in these regions were not lasting sufficiently long (a few months, not shown) for the adaptive BC to spin-up.

In Figure 3 the mean soil moisture increments are shown for the surface layer (0-7 cm depth) for August 2022. The patterns are generally similar for $\mathrm{C}$ and ${\mathrm{E}}_{\mathrm{A},\mathrm{S}}$. However, in a latitude band from about ${55}^{\circ }$ to ${80}^{\circ }$ North there are smaller magnitude increments in ${\mathrm{E}}_{\mathrm{A},\mathrm{S}}$ compared to $\mathrm{C}$. In particular, northeast Asia shows a reduced positive signal when the adaptive BC is active. However, the impact of the adaptive BC on the increments in the lower latitudes ($<{55}^{\circ }$ North) is generally small, even in regions where the BC is large e.g. South America. In some cases for SMOS SM this is related to larger observation errors being prescribed in regions known to have quality control issues, such as the Andes mountain range. For the RZSM (top metre of soil) the depth-integrated increments are also impacted in the higher latitudes by the BC, but to a lesser extent than the surface layer (not shown).

Focusing on a point in central-eastern Australia with latitude ${25}^{\circ }\mathrm{S}$ and longitude ${140}^{\circ }\mathrm{E}$, Figure 4 shows the time series of the SMOS SM departures for (a) $\mathrm{C}$ and (c) ${\mathrm{E}}_{\mathrm{A},\mathrm{S}}$. Over the 4-year period the departures are almost entirely positive for $\mathrm{C}$, which demonstrates a local positive bias in the SMOS NN observations with respect to the model. The time series is effectively shifted in the negative direction by the BC, while preserving the temporal variability in the departures. Figure 4 (b) and (d) present the distribution of the innovations normalized by $(\mathbf{H}\mathbf{B}{\mathbf{H}}^T+\mathbf{R})$ according to Equation (1) of Desroziers et al. [53]. The positive bias in $\mathrm{C}$ and the negative shift in the distribution for ${\mathrm{E}}_{\mathrm{A},\mathrm{S}}$ are clearly evident. Neither distribution follows a Gaussian curve, which might be related to the nonlinear nature of the land surface model and the associated background and observation uncertainties.

Table 3 summarises the global average internal diagnostics for the five experiments. The results suggest that the ASCAT and SMOS biases are mainly independent of each other as the BC of SMOS alone (${\mathrm{E}}_{\mathrm{S}}$) does little to reduce the ASCAT SM departures and vice-versa when assimilating ASCAT alone (${\mathrm{E}}_{\mathrm{A}}$). Experiment ${\mathrm{E}}_{\mathrm{A},\mathrm{S}}*$ has larger absolute ASCAT SM departures than $\mathrm{C}$, which suggests that the ASCAT adaptive BC and ASCAT seasonal rescaling complement each other. The mean SSM increments are reduced by the combination of adaptive BC and ASCAT seasonal rescaling. Furthermore, the absolute values are slightly reduced. However, the impact of the adaptive BC on the RZSM increments is relatively small compared to the SSM increments.

3.2. SM and ST validation

The analysed SM and ST are validated using observations from the ISMN over 2020-2022. Figure 5 shows the locations of the stations for each network over the US, Australia, France, Spain and Germany. Also shown is the relative performance of the Pearson R anomaly for ${\mathrm{E}}_{\mathrm{A},\mathrm{S}}$ compared to C for the SSM. Overall, $1.7\%$ of all the stations in ${\mathrm{E}}_{\mathrm{A},\mathrm{S}}$ significantly improved (at the $95\%$ confidence level) on $\mathrm{C}$ and $0.5\%$ of stations were significantly degraded relative to $\mathrm{C}$, with all these sites located in the US (see Figure 5(a)). Table 4 gives the global average SSM and RZSM scores for the different experiments. There is a slight improvement in the SSM and RZSM anomaly correlations for the adaptive BC experiments compared to the control. The other scores are largely equivalent when averaged over all the networks. Whilst the impact of the BC is relatively small for most stations, there are some strong local impacts. For example, Figure 6 shows the time series of the SSM analysis for the different experiments at a point located in the USCRN network in the US, which corresponds to the station circled in Figure 5. In this case the adaptive BC reduces the positive bias in the experiments and subsequently improves the fit to the observations, especially when both ASCAT and SMOS are bias-corrected (${\mathrm{E}}_{\mathrm{A},\mathrm{S}}$).

3.3. Atmospheric validation

Figure 7 presents the global difference in normalized relative humidity RMSE (dRMSE) for ${\mathrm{E}}_{\mathrm{S}}$ and ${\mathrm{E}}_{\mathrm{A}}$ with respect to $\mathrm{C}$, at 1000 hPa and for forecast lead times of 1-10 days. Evidently the impact of the adaptive BC is small in both cases and over the different regions, with the SMOS adaptive BC giving a slightly positive impact over the southern extratropics and the tropics and the ASCAT adaptive BC giving a slightly negative impact over the northern hemisphere extratropics. Figure 8 shows a latitude-pressure plot of the difference in mean errors (dME) for ${\mathrm{E}}_{\mathrm{S}}$ (top), ${\mathrm{E}}_{\mathrm{A},\mathrm{S}}$ (middle) and ${\mathrm{E}}_{\mathrm{A},\mathrm{S}}*$ (bottom) for a lead times of 72 hours. The mean errors for ${\mathrm{E}}_{\mathrm{S}}$ are significantly improved in the lower troposphere (above 850 hPa) between 60 and 70 degrees North. The pattern is similar for ${\mathrm{E}}_{\mathrm{A},\mathrm{S}}$ but the signal is spread over a larger latitude band between 30 and 70 degrees north. On the other hand, the signal is reduced going from ${\mathrm{E}}_{\mathrm{A},\mathrm{S}}$ to ${\mathrm{E}}_{\mathrm{A},\mathrm{S}}*$. Together these results demonstrate the complementary impacts of the adaptive BC and the ASCAT seasonal rescaling on the mean relative humidity in the boundary layer for the northern hemisphere midlatitudes. The dME improvements persisted to about day 5 (not shown). Figure 9 shows the mean differences in RH between ${\mathrm{E}}_{\mathrm{A},\mathrm{S}}$ and $\mathrm{C}$. There are reductions in relative humidity over eastern China and South America, which roughly correspond with the BC seen in Figure 2 (c) and (f) respectively. Likewise, the positive difference over the southeastern US agrees with the positive BC in Figure 2 (f). In the case of atmospheric temperature, neither the dME nor the dRMSE were significantly impacted by the adaptive BC (not shown).

4. Summary and discussion

In this study, a two-stage observation bias-correction and state update filter was introduced in the ECMWF soil moisture analysis and evaluated over the period 2020-2022. Compared with the existing soil moisture analysis, the two-stage filter is computationally efficient and only requires the additional storage of the ASCAT and SMOS biases for each gridpoint. As in Draper et al. [22], the observation-model bias is effectively treated as the observation bias.

The impact of the two-stage bias filter on the internal DA diagnostics was assessed for ASCAT and SMOS SM by evaluating the first guess (observation-model) departures and the SM increments. It should be recognised that the current seasonal rescaling for ASCAT SM accounts for stationary biases in the mean and variance. On the other hand, the adaptive BC accounts for non-stationary biases but it only corrects the mean soil moisture state. Therefore the adaptive BC alone could not outperform the seasonal rescaling approach. However, when the adaptive BC was employed on top of the seasonal rescaling the magnitude of the first guess departures was reduced by up to $20\%$ for ASCAT SM, with the greatest impact when the assimilation was most active during the northern hemisphere summer. Furthermore, up to $6\%$ more ASCAT observations were assimilated as excessive first guess departures were removed less often by the quality control.

Over the four-year experiment period, a slight positive global trend was found for the ASCAT SM departures. Further analysis has demonstrated that the positive trend can be largely attributed to moistening ASCAT observations (not shown). Whilst 4 years is a relatively short timescale, these results appear to be consistent with Hahn et al. [54], who demonstrated that the backscatter signal from ASCAT SM observations has been increasing in some regions over the last 15 years. This is partly related to land-use change, including urbanisation and deforestation. They found that a periodic recalibration of the dry and wet backscatter reference can help to mitigate the spurious trend, and future ASCAT soil moisture products will be recalibrated accordingly. In this study the adaptive SM BC effectively removed a global positive trend in ASCAT SM.

In the case of SMOS, large local soil moisture biases in the first guess departures were found over several regions, including a strong positive bias over the Andes in South America and over parts of Australia, and a negative bias over the eastern US and high latitude regions in Asia. The SMOS neural network is trained globally, which may partly explain these results. The biases could be alleviated by introducing local predictors in the training dataset, such as latitude and longitude, which will be investigated in future versions of the neural network. Nevertheless, the adaptive BC corrected many of these local biases, which was illustrated for a persistent positive bias over a point in eastern Australia. Globally the adaptive BC reduced the magnitude of the first guess departures by up to $30\%$ with typically $10\%$ more observations assimilated in the northern hemisphere summer months.

The impact of the adaptive BC on soil moisture and soil temperature performance was evaluated using observations from the international soil moisture network. On average, the adaptive BC slightly improved the anomaly correlation of the surface and root-zone SM with the observations. There were also some locally significant impacts for individual stations. However, many regions could not be validated due to the lack of in situ data, including northeast Asia, where some of the largest biases were corrected. Nevertheless, the adaptive BC experiments demonstrated statistically significant improvements over the control in the 1- to 5-day forecasts of mean relative humidity in the lower boundary layer, mainly over mid to high latitudes in the northern hemisphere.

Whilst it is acknowledged that the perfect model assumption in this study is likely to be incorrect, the heterogeneity of soil moisture together with the lack of an accurate ground reference makes it very challenging to differentiate between model and observation errors in global land data assimilation (DA) systems. Furthermore, land surface models rely on calibrated parameterisations to relate soil moisture prognostic variables to model fluxes, such as hydraulic conductivity and latent/sensible heat exchanges with the boundary layer. Arguably, accounting for modelled soil moisture biases can only be achieved by correcting the model parameters in conjunction with the state estimates [55,56]. Joint initiatives between ECMWF and other organisations are working to improve land surface benchmarking datasets, such as the PLUMBER2 model inter-comparison project [57].

Finally, it is planned that the adaptive BC will be implemented operationally at ECMWF. A simplified version of the IFS was employed for the experiments in this study in order to minimise the computational cost, which is called the “stand-alone surface analysis" (SSA). Although Fairbairn et al. [30] demonstrated that SSA is a useful research tool, a further evaluation of the adaptive SM BC will be carried out using the full weakly coupled land-atmosphere DA system.