remote sensing Assimilation of Satellite-Derived Soil Moisture and Brightness Temperature in Land Surface Models: A Review

: The correction of Soil Moisture (SM) estimates in Land Surface Models (LSMs) is considered essential for improving the performance of numerical weather forecasting and hydrologic models used in weather and climate studies. Along with surface screen-level variables, the satellite data, including Brightness Temperature (BT) from passive microwave sensors, and retrieved SM from active, passive, or combined active–passive sensor products have been used as two critical inputs in improvements of the LSM. The present study reviewed the current status in correcting LSM SM estimates, evaluating the results with in situ measurements. Based on ﬁndings from previous studies, a detailed analysis of related issues in the assimilation of SM in LSM, including bias correction of satellite data, applied LSMs and in situ observations, input data from various satellite sensors, sources of errors, calibration (both LSM and radiative transfer model), are discussed. Moreover, assimilation approaches are compared, and considerations for assimilation implementation are presented. A quantitative representation of results from the literature review, including ranges and variability of improvements in LSMs due to assimilation, are analyzed for both surface and root zone SM. A direction for future studies is then presented.


Introduction
Soil Moisture (SM) is a crucial variable in the partitioning of water (into infiltration and runoff) and energy (into sensible and latent heat flux) and has been considered in many atmospheric and hydrological studies [1,2]. SM influences the carbon cycling by photosynthesis and transpiration processes [3] and affects short and medium-range forecasts in meteorology [4]. The variability in SM is due to the spatial distribution of many factors such as rainfall, the spatial variation of the soil texture, vegetation, and topography. Volumetric SM measurements at different depths are limited to sparse ground measurements, which cannot provide required spatial and temporal resolution data for numerical weather initialization. To solve this problem, Land Surface Models (LSMs) are used to model the near-surface (~0-5 cm) and root zone (~5 to 100) SM (hereafter referred to as SSM and RSM) behavior through physical and hydrological laws. However, the estimation of SM by LSMs is limited by forcing data, built-in physical rules, and parameterizations.
In order to reduce the mentioned uncertainties, the assimilation of surface observations in LSMs was proposed. Generally, the assimilation procedure in LSMs consists of three main The review paper is organized as follows. In Section 2, the current LSM SM improvement procedures are presented, and the uncertainties in LSMs, satellite BT, and SM for assimilation are reported. The correction of bias between the LSM and satellite is also reviewed. This section also refers to the applied assimilation methods regarding their characteristics and limitations. Section 3 presents the outcomes of studies based on LSM classification, and Section 4 identifies some of the issues with current approaches assimilating SM and BT in LSMs and offers recommendations for improvement.

Procedures for the Improvement of SM in LSMs
Based on the reviewed papers, when using both pieces of satellite data inputs (BT and retrieved SM), the correction of LSM SM estimates can be divided into the following categories: (A) some papers followed the correction of LSM climatology by calibration (optimization of the fixed parameters) of the LSM, (B) others conducted the correction of forcing data (required data for initialization of LSMs) by assimilation. In Figure 2, the procedures for both calibration (approaches 1 and 2 in Figure 2) and assimilation (approaches 3 and 4 in Figure 2) in the LSM are presented. In approach 1, the LSM SM estimates are applied to simulate the BTs (through RTM), and then an optimization method is used to decrease the difference between the simulated and observed BTs, and hence to optimize the critical parameters in the LSM (e.g., [24]). In approach 2, the optimization of LSM parameters is conducted by comparing the estimated LSM SM to satellite retrieved SM products (e.g., [25]). In approach 3, RTM is calibrated by adjusting LSM SM to reduce the difference between the observed and simulated BTs (e.g., [26][27][28]). Systems such as the Coupled Land and Vegetation Data Assimilation System (CLVDAS) were designed (as a dual-pass assimilation system), which simultaneously calibrates the parameters and assimilates the BT in the LSM [29]. In approach 4, satellite SM is used through assimilation to correct the LSM SM estimates. In both approaches (calibration and assimilation), RTM calibration can be conducted.

Land Surface Models
Over the last decade, 2 m relative humidity, air temperature (from ground stations), and BT and SM products (from satellites) have been assimilated in different LSMs to improve SM estimates. Figure 3 shows different LSMs (total number: 73) applied in SM assimilations across the years in reviewed papers. According to Figure 3, the number of studies has increased since 2007, and most of these were conducted between 2013 and 2017. Most of the assimilation studies were implemented in Noah, CLSM, and the Community Land Model (CLM).  Table 1).
The applied LSMs along with their default soil layers are shown in Table 1. Table 1. The layering assumptions in LSMs; **: depends on land cover, tall vegetation fraction 0.15 and 1, bare soil 0.15. *: Discretized to 5.

LSM Number and Depth of Soil Layers (m) Reference
Noah 4 0.1, 0.3, 0.6, 1.0 [30,31] Australian Water Resources Assessment system (AWRA-L) 3 0.1, 0.2, and 8 [32] Catchment Land Surface Model (CLSM) 2 0.02, 1.0 [33] Community The sources of errors in LSMs include a large number of their parameters influenced by the vegetation and soil state estimation [48], model physics and structure, forcing data, and boundary conditions [49]. The role of precipitation as forcing data for assimilation was proven to be significant (discussed further). Figure 4 shows the frequency of applied forcing data for the assimilation of SM in LSMs in the reviewed papers, showing that Modern-Era Retrospective analysis for Research and Applications (MERRA) [50] was used the most. Blyverket et al. [51] compared the performance of MERRA-2 and NLDAS (with precipitation correction) over the contiguous United States (US). They suggested that the assimilation procedure has higher improvement skills with MERRA-2; the correction of precipitation in NLDAS may improve the results of SM estimation, and hence reduce the assimilation skill [52].

RTM
RTMs include a large number of parameters that are highly variable and sometimes unavailable [53] as well as showing dependency to land cover [28]. See Section 2.3 for more details.

Satellite Data
The inherent uncertainties in both BT and SM products that may influence assimilation are:

1.
SM retrieval approaches. Different strategies and simplifications applied for parameterization of RTMs and soil dielectric mixing models in SM retrieval. As stated in [54,55], most of the SM retrieval algorithms used T-ω zero-order RTM as the baseline. The major difference between these approaches originates from their parameterizations and estimation of vegetation optical depth. The modeling of roughness, soil, and canopy temperatures, vegetation structure and its optical depth, and atmospheric effects are the most significant factors in SM retrieval. Choosing between effective single or multiple scattering albedo, single or double polarization, static or dynamic roughness coefficients, and single-or multi-angle satellite observations, the applied frequency and ancillary data make the algorithms complex and define the quality of SM products. Moreover, as discussed in [54], the selection of mixing dielectric modeling is another critical issue that needs to be regarded when using SM products and evaluating assimilation results.

2.
The depth of satellite-derived SM. The depth of measurement is not definite and is a function of moisture, temperature, and soil texture. The introduced depths vary spatially and temporally [23,56,57].

3.
The time of satellite orbital pass. Some studies have investigated the accuracy of the SM products of a sensor at different acquisition times. For example, Martens et al. [58] and Yee et al. [59] showed that ascending SMOS data (about 6:00 a.m.) were more accurate than the descending pass (about 6:00 p.m.). However, Jing et al. [60] maintained that the descending pass of SMOS performed better than the ascending pass. In [17,61], it was stated that ASCAT SM data from the descending pass (9:30 a.m.) have more accuracy than the ascending one (9:30 p.m.). Jing et al. [60] noted that the AMSR-E ascending (1:30 p.m.) product was better than the AMSR-E descending (1:30 a.m.) product. Yee et al. [59] stated that AMSR2 X-band products have better performance in evening overpasses (1:30 pm) in contrast with morning overpasses (1:30 a.m.). As stated in [62], it is expected that night time or early morning SM products be more accurate due to: (a) assumed temperature equilibrium between soil and vegetation canopy in the SM retrieval modeling [63] and (b) expected minimum of Faraday rotations during the night [64,65]. However, as indicated above in some studies, other revisit times showed more accuracy, which can be attributed to other factors [62]. These are: (a) Radio Frequency Interference (RFI) contamination in some time periods and (b) diurnal variability of surface conditions such as water in vegetation. Therefore, it is expected that the accuracy of SM products in different acquisition times be checked in the study area before conducting an assimilation study. 4.
Land cover. The vegetation coverage and type affect the SM retrieval, especially at higher frequencies.
In [17,61], it was reported that X and C bands are opaque in dense forests and shrubs with green vegetation fractions of more than 0.5. Other land covers, such as frozen ground, snow cover, water body areas (e.g., flooded areas, rivers, wetlands, lakes, precipitation at the time of satellite passes), steep topography (e.g., rocks), urban area, and heavily forested areas [66] would decrease the quality flag of data recording and, as a result, reduce the available data for assimilation. 5.
Temperature limitation. Due to the freezing temperature, currently, no reliable satellitederived SM product is available for low-temperature latitudes and snow covers [67]. 6.
Temporal resolution. The higher the temporal resolution, the more accurate SM temporal dynamics (e.g., caused by irrigation) can be retrieved. Yin et al. [68] proposed that SMOPS products (as combined products) have a higher temporal resolution, which can capture SM dynamics with more detail. 7.
Spatial resolution. The low spatial resolution of passive microwave sensors, the mismatch between the LSM and satellite spatial resolutions, and the unmolded sub-pixel heterogeneity (mixed land covers mentioned in item 3, land cover, within a pixel) could increase the errors. Downscaling the SM products can be a solution, but the accuracy of downscaling depends on applied ancillary data, downscaling the method and accuracy of in situ measurements, which needs to be analyzed [69]. Therefore, in some cases, it may bring forth lower results. For example, Lievens et al. [70] compared the performance of downscaled SM (by MODIS thermal data) with coarse resolution and showed that downscaling could not provide better results. 8.
Polarization. Vertical and horizontal polarization might not retrieve the same BTs; vertical polarization BTs are less sensitive to vegetation heterogeneity and roughness than horizontal polarization while they are more sensitive to SM [4,70,71]. Tian et al. [26] showed that the horizontal polarization BTs from AMSR-E have lower temporal variability than the vertical, which can make them less suitable for assimilation experiments. The role of incidence angle could make the problem more complex; Lievens et al. [70] showed that horizontal polarization with incidence angle less than 42.5 degrees could provide better results than in combination with vertical polarization and multi-angle observations (with SMOS sensor). They stated that it could be attributed to the correlation between dual-polarization and multi-angle observations. The higher sensitivity of horizontal polarization to SM and the lower sensitivity of vertical polarization to vegetation and roughness are issues that require more research to find out the weight of each one of these characteristics. 9.
Dynamic drying rate. The drying rate of estimated SMs in LSMs could be different from satellite-derived SMs. Shellito et al. [72] compared the drying rates of SMAP and those simulated by Noah and showed that the drying SSM from SMAP is faster than Noah simulations. They also stated that when SM content (Noah and SMAP SM contents are linearly related) is high, potential evaporation is high, vegetation cover is low, and SMAP drying is the fastest. They reported that the effect of vegetation on the Noah model is simplified and is less than SMs retrieved by SMAP.
To find out the biases in satellite products, several studies compared them with ground networks. Although most of them stated that SMAP, SMOS-IC (V105 and V2), and CCI (based on LSM corrections) have slightly better performance [73][74][75], it should be considered that SM product biases differ based on soil and vegetation conditions. Figure 5 shows the applied satellite data in assimilation studies. As Figure 5 shows, several microwave data and an algorithm (ALEXI, Atmosphere-Land Exchange Inverse, based on the thermal infrared band and used for comparison with microwave data) were applied for the assimilation studies, and the AMSR-E was used the most. The increasing trend in the application of such data requires that the best products are identified based on the spatio-temporal conditions of the study area and the characteristics of applied evaluation methods [76].

Calibration of LSM and RTM Parameters
As LSM and RTM parameter values are assumed to be fixed, their resultant error in the assimilation process is inevitable. It is essential that the biases between RTM and LSM (through simulated BT) and satellite data decrease [77,78]. Therefore, these parameters in both models should be calibrated by comparison with satellite-observed BT. In parameter calibration (optimization), several issues need to be considered, including the applied methods for optimization, the selection of parameters to be optimized, uncertainties in calibration, time window setting, and the number of used microwave frequencies in the calibration process.
The applied methods for the optimization of parameters in LSM and RTM are: (1) Shuffled Complex Evolution (SCE-UA); this minimizes a cost function composed of observed and simulated BT (by LSM or RTM), wavelength, and polarization of observation and then optimizes the parameters in their value domain using complex evolution algorithm [26,[79][80][81]. (2) Ensemble Proper Orthogonal Decomposition-based Parameter (EnPOD_P); this constructs an ensemble of state variables (by the Monte Carlo method), including the perturbed states (simulated BTs) and perturbed parameters at each time of satellite BT, and finally optimizes a cost function and obtains the parameters [49,53,82].
(3) Particle Swarm Optimization; this uses multi-angular SMOS BTs and tries to minimize a cost function composed of the penalty terms for bias between simulated BTs and SMOS BTs, the variability (temporal standard deviation) and the parameters by Particle Swarm Optimization algorithm [83]. (4) Particle Filtering; this simultaneously estimates SM and some LSM parameters using state augmentation method [84]. (5) Markov Chain Monte Carlo (MCMC); this uses MCMC algorithm to estimate RTM parameters by minimizing the differences between simulated BTs and multi-angle SMOS BTs [21]. Tian et al. [82] showed that the EnPOD_P considerably outperforms the SCE-UA due to the simultaneous optimization of the model states and parameters. Furthermore, De Lannoy et al. [21] indicated that the parameter values estimated through Particle Swarm Optimization were in close agreement with those obtained from the MCMC simulation. Further research is essential to compare the performance of these optimization methods for both LSMs and RTMs in different conditions.
The selection of parameters for calibration is another crucial issue that is used in both RTM and LSM. Generally, the sensitive parameters in RTM to BT (microwave emissivity) are mainly from vegetation and surface roughness parameters, soil structure, and texture [77,85,86]. Based on the conducted sensitivity analysis in observation operators (Land Emissivity Model [87], Community Microwave Emission Model [88], and T-ω) [89], these parameters include the following. (1) For vegetation: leaf thickness, affecting vegetation optical thickness [53,90], vegetation structure parameter, affecting vegetation optical thickness [53,91], the radius of dense medium scatters, fraction volume of dense medium scatters [90], and b parameter, affecting vegetation optical depth in the T-ω model [91].
In LSM, the most sensitive parameters to SM also need to be considered. Sawada et al. [29] and Sawada et al. [95] conducted a sensitivity analysis and selected the most sensitive parameters (both hydrological and ecological). These parameters include saturated hydraulic conductivity, saturation water content (porosity), parameters in van Genuchten's water retention curve, the maximum rubisco capacity of top leaf, the normal turnover rate of leaves, the empirical parameter of root biomass fraction function, point of water-related stress loss, and point of temperature-related stress loss. After choosing the sensitive parameters, the works of [29,78,97] optimized the selected parameters by comparing simulated BTs and satellite-observed BTs. Koster et al. [27] applied a different method, where a time-invariant parameter was determined for calibration of LSM SM estimates, and through comparison with SMAP SM observations, the parameter value was determined.
Many factors should be considered during calibration. In RTM calibration, Zhao et al. [77] maintained that, due to the limited number of iterations in the optimization process, it is possible that calibrated parameters are not efficiently optimized. This can produce degraded results in contrast to non-calibrated parameters. Moreover, another issue is the consideration of land cover in calibration. For example, De Lannoy et al. [28] showed that the calibration of the RTM revealed higher errors in BT predictions over cropland. In LSM calibration, Yang et al. [97] recognized that the estimated parameter values might not be the truth values of soil parameters because, firstly, dense vegetation could reduce the model calibration effect; in areas with dense vegetation, microwave signals provide more information about vegetation rather than soil state. This may reduce the accuracy of the soil model calibration. Secondly, estimated parameter values may be biased due to inappropriate parameterizations or simplified physical assumptions in an LSM, such as the partitioning between runoff generation and SM content and the relationship between surface temperature and SM. Furthermore, current RTMs and LSMs are usually calibrated by passive microwave satellites with pixel sizes of more than 10 km. The issues of spatial matching between models, satellite data, and heterogeneity within each pixel could cause errors in the calibration process.
The time window for both RTM and LSM calibration is another issue that needs to be considered. For RTM, no study analyzed the effect of time window on the calibration process in the analyzed period of this study, but for LSM, [97] showed that a time window of 1 year for the model calibration might be reasonable for grasslands. More studies are needed to analyze the effect of a time window.
In most studies for RTM and LSM calibration, one satellite frequency was applied, but Sawada et al. [29] applied a cost function that used three frequencies to provide more stability as it could provide more robust results.

Bias Correction in SM Assimilation
The systematic error or bias removal between the LSM and satellite-derived SM should be considered for the assimilation procedure. The applied strategies to deal with these biases are shown in Figure 6. Several online and a priori methods have been introduced for bias correction between the LSM, in situ measurements, and satellite data [98][99][100][101][102][103][104][105][106]. In online (dynamic) bias correction, the biases are corrected in the assimilation process. This could be used when transient biases appear and they could not be identified in the climatology of both model and observation. In this case, the bias is attributed to either model [103,104] or observations [107].
The applied a priori bias correction methods (bias correction prior to assimilation) for satellite SM are Triple Collocation (TC), linear scaling (linear regression), and Cumulative Distribution Function (CDF). In TC, it is possible to use three independent datasets to remove the bias. Yilmaz et al. [101] evaluated the performance of least squares regression, variance matching, and TC. They showed that TC had better results compared to others.
Most of the studies reviewed here used the CDF technique for bias removal in SM assimilation. In CDF matching, the CDFs of both observation and model at each grid cell are extracted, and the observation values are bias-corrected to the model climatology through the mean, variance, and possibly higher-order moments; Lievens et al. [108] emphasized the mean as the most accurate rescaling method. A disadvantage of this method is that it needs a long data record to provide a suitable sample size. Reichle et al. [109] proposed the replacement of spatial sampling with temporal sampling to reduce the required time scale to one year. Regarding the required long time series for CDF matching, Mahfouf [110] proposed considering all points in the study area (global rescaling). Such a method was compared with localized recalling (point-to-point matching between model and satellite observations) in [111,112]. Kolassa et al. [111] showed that the global approach had higher skills in improving SM estimation by having a higher impact of observation on assimilation procedure. They stated that this approach is subjected to these conditions: (1) a good observation error characterization; (2) rigorous observation quality control; (3) potential component-wise assimilation to better isolate the reliable satellite information. Despite these results, in the point-wise approach, there is more sensitivity to locations where observation is uncertain or the range of change is high (mountainous areas). Schneider et al. [112] also showed the influential role of domain spatial diversity in models where global rescaling could remove important observation data. Blankenship et al. [113] and Blankenship et al. [114] used soil type-based CDF matching and showed that this method had better performance than localized matching. It needs to be considered that the spatial resolution (or subcategories) of soil types is an important variable. Blankenship et al. [114] also showed that vegetation could not have significant results in their applied CDF matching in the warm season of 2011 in central and southeastern United States.
According to some recent studies, from a temporal point of view, CDF matching can be conducted in two ways: lumped and seasonal [98,115]. In the lumped approach (one year CDF curve for all points), it is assumed that the difference between the model and observations is stationary in time. In the seasonal (or monthly) approach, the biases are corrected through temporal stratification. In the case of unmolded processes such as irrigation, the seasonal method would be more effective but it requires an appropriate sample size [116][117][118][119]. Yin et al. [120] compared the performance of lumped CDF matching, monthly CDF matching, and linear matching of monthly mean and standard deviation of SMAP and model SM simulations. They showed that monthly CDF matching had the best results in the prediction of precipitation.
As it is possible that in some situations (such as localized bias correction) the inaccurate model climatology is forced toward satellite one, some studies proposed a correction of the model toward observations. This would be applicable when the model biases are spatially known.
Some studies used bias-blind assimilation. In this case, it is assumed that observations and model estimates are not biased. Lin et al. [121] believe that all available data sets for assimilation should be kept intact to use the useful information within the data. Rescaling the observation to the model climatology could remove valuable information [113]. This is more critical when the sources of biases are unknown in both observation and the model.  Table 2 lists and describes the techniques used in these studies, while Figure 7 shows the specific methods used in the studies reviewed here. This figure also shows the increase in the number of studies using assimilation techniques (in particular, from 2014 to 2016), with most of the studies using EnKF (Ensemble Kalman Filter), as based on reviewed papers. The model states and error statistics are generated by perturbations in the model initial conditions (causing a set of ensembles) and are forwarded by a non-linear model without the need for linearization. [127]  11 Evolutionary (Evol) Relies on evolutionary algorithms to provide evolved members for analysis step in assimilation. [131] 12 Particle Filter (PF) The probability density function of model SM estimates is represented by particles; a set of weighted trajectories (based on observations) are defined between background and analysis steps. [132] 15 Newtonian Nudging (NNudg) The model is dynamically relaxed toward the observations, and therefore the model is considered as a weak constraint and observations are considered as perfect. [133,134] 16 Optimal Interpolation (OI) Based on the known errors, observations are weighted and the analysis is estimated based on the gain matrix to obtain the analysis. [135] 17 Variational (3D and 4D) (Var) Model state is estimated by minimization of cost function, which explains the misfit between the model simulations and observations. If the cost function is conducted in a time interval, it is converted to 4D Var type. [136]

Comparison of Assimilation Methods in Improvement of SM Estimates
Over the last decade, several studies have compared the performance of different assimilation methods on the improvement of SM in LSMs. De Rosnay et al. [5] replaced the Optimal Interpolation (OI) method with EKF using 2 m air temperature and relative humidity observations. They tested the possible replacement of ASCAT SM. They showed that EKF had superior performance over OI (as it accounts for the nonlinear relationship of meteorological forcing and SM). Dumedah et al. [137] used SMOS SM data and found that the surface estimation accuracy of SM for EnKF and Evolutionary Data Assimilation (EDA) assimilation methods increased, and EDA had better results than EnKF when evaluating with SMOS level 2 SMs. EDA could converge the model parameters and improve the SM estimation. Blyverket et al. [51] compared the performance of EnKF, EnOI, and climatological background error (EnOI-Clim) assimilation methods using Climate Change Initiative (ESA CCI) by assimilating SM retrievals from SMAP and SMOS over the contiguous US. They showed that, although EnKF had slightly better performance, EnOI is recommended because it is computationally cheaper. The disadvantage of EnOI is that it does not use flow-dependent ensemble spread. Comparisons of the EKF and EnKF showed that EnKF had slightly better results in the southeastern United States using synthetic observations [138] and over an experimental site in southwestern France [139]. Bi et al. [140] compared the performance of EnKF, the standard PF, and the Improved PF (IPF) using brightness temperatures from the AMSR-E into the VIC LSM in the NaQu network region at the Tibetan Plateau. The IPF improved the estimation along with requiring fewer particles (explaining the probability density function of the model by a set of random samples) than EnKF and the standard PF.
Yan et al. [141] analyzed the estimation of SM and hydrologic properties by PF (used AMSR-E) in the Oklahoma Little Washita Watershed. The assimilation methods were PF with commonly used sampling importance resampling (PF-SIR) and Markov chain Monte Carlo sampling (PFMCMC). The applied assimilation methods were overestimated in the dry season and underestimated in the wet season. PF-MCMC showed less error than PF-SIR simulations.

Optimal Ensemble Size in EnKF and PF
To conduct EnKF and its variants, it is necessary to determine the optimal ensemble members for calculating the covariance matrices. The determination of ensemble size has been studied in [142][143][144]. The number of ensemble members and particles that were considered in the analyzed papers are shown in Figure 8 (based on available data). For EnKF (and its variants), most studies used 12 ensembles. In the En4DVAR approach, the numbers vary between 30 and 100. For the PF method, the minimum number was 30 and the maximum was 1000. However, the optimal ensemble size may vary depending on the study and, more importantly, the computation cost.

The Results of Assimilating BT and SM in LSMs
In the following section, the results of assimilations (for both BT and SM) are categorized based on LSMs considering the papers reviewed in this work. Figures 9 and 10 compare the performance of non-assimilation cases with respect to assimilation cases based on bias ( Figure 9A), correlation ( Figure 9B), RMSE ( Figure 9C), and unbiased RMSE (uRMSE) ( Figure 9D) computed between the SM estimated in LSM and ground SM networks for SSM and RSM, respectively. Non-assimilation cases include those bias-corrected/non-bias corrected. In Figure 9A, which shows the bias for NOAH, Jules, and VIC, the SM bias range varies between −0.096 and 0.13 (m3 m-3) for non-assimilation studies and between −0.077 and 0.069 (m3 m-3) for assimilation studies. The NOAH has the greater variability in both cases, with assimilation cases having lower biases. In Figure 9B, the performances of NOAH, CLSM, CLM, Jules, ISBA, SiB2, and VIC are shown for a correlation index. The correlation varies between −0.243 and 0.95 in non-assimilations and between 0.343 and 0.93 for assimilations. The mean of assimilation (0.69) cases is higher than non-assimilation (0.64) cases. In Figure 9C, the RMSE in NOAH, CLM, Jules, SiB2, and VIC varies between 0.016 and 0.13 (m3 m-3) in non-assimilation cases and between 0.01 and 0.15 (m3 m-3) for assimilation cases. The mean of RMSE (0.058 m3 m-3) is lower in assimilation cases with respect to non-assimilations (0.069 m3 m-3). In Figure 9D, the uRMSE in NOAH and CLSM varies between 0.018 and 0.079 (m3 m-3) in non-assimilations and between 0.012 and 0.075 (m3 m-3) in assimilations. The uRMSE mean of assimilations (0.041 m3 m-3) is lower than non-assimilations (0.046 m3 m-3). Figure 10 is the same as Figure 9 but regarding RSM. In Figure 10A, the performances of NOAH and Jules in the bias index are shown. The range for non-assimilation cases varies between 0.025 and 0.1 (m3 m-3) and between −0.017 and 0.04 (m3 m-3) for assimilation cases. The mean for assimilation cases (0.005 m3 m-3) is lower compared with non-assimilation cases (0.054 m3 m-3). In Figure 10B, the correlations for NOAH, CLSM, Jules, ISBA, and SiB2 are shown. The ranges vary between 0.25 and 0.89 in non-assimilations and between 0.35 and 0.89 in assimilations. In Figure 10C, the RMSE is shown for NOAH, CLM, Jules, and SiB2, varying between 0.022 and 0.111 (m3 m-3) for non-assimilations and between 0.024 and 0.109 (m3 m-3) for assimilations. In Figure 10D, the uRMSE is shown for NOAH and CLSM. The ranges vary between 0.027 and 0.053 (m3 m-3) for non-assimilations and between 0.025 and 0.048 (m3 m-3) for assimilations. The uRMSE mean of assimilations (0.036 m3 m-3) is lower than that of non-assimilations (0.040 m3 m-3). Such results (both in SSM and RSM) are shown to demonstrate the quantitative range of error indices in the reviewed papers. These ranges are a function of many factors, such as the accuracy of in situ measurements, study area, the applied sensors, assimilation methods, etc. To show the effect of the assimilation procedures, the differences between assimilation and non-assimilation cases are reported (sorted based on correlation) for SSM ( Figure 11). Some authors conducted multiple approaches, which are shown as slots with different colors in each bar; for example, [29] conducted two assimilation scenarios with different configurations. As Figure 11 shows, most of assimilation approaches improved the SM estimation. In some papers [29,78,97], some configurations resulted in negative correlation differences where only one station was used, which cannot be considered a robust result. Cases providing correlation improvements of more than 0.12 are presented in Table 3.
Most of these scenarios used AMSR-E and SMOS in their assimilation studies. Overall, when BT and SM are used as input, no significant priority can be found between them. The highest improvement was observed in [53] with CLM LSM and the En4DVAR approach in China. Most configurations used EnKF for conducting the assimilation.
Based on Figure 11, RMSE improvements for assimilation cases vary between −0.002 and −0.07 (m3 m-3) (excluding positive values of retrogression). Table 4 shows RMSE improvements in SSM with RMSEs less than −0.024 (m3 m-3). These studies applied AMSR-E and SMOS BT and SM products. Out of ten configurations, five scenarios used the EnKF approach for assimilation.
From Figure 11, bias varies between −0.003 and −0.09 (m3 m-3) (excluding positive values of retrogression). Table 5 demonstrates the bias improvements in SSM with biases less than −0.023 (m3 m-3). These studies used AMSR-E and SMOS and EnKF, which is the most frequent assimilation approach.
In Figure 12, the improvements in correlation, RMSE, and bias are shown for RSM. In most cases, the correlation was improved. The correlation improvements vary between 0.005 and 0.184 (excluding negative values). In [29], a significant negative correlation was found, which is the subtraction of assimilation from open loop (the open loop is initialized with 35 days of assimilation results). It shows that when the Sib2 LSM was initialized with 35 days of assimilation in this study, better results were obtained, and the open-loop run in this configuration had better performance than assimilation in the whole study period. Regarding RMSE, most of the studies reported an improvement between −0.007 and −0.04 (m3 m-3) (excluding positive values). All available studies show bias improvements in RSM by assimilation. In Tables 6-8, the configuration of high score results are shown. Figure 11. The difference between metrics in assimilation and non-assimilation scenarios for SSM in terms of correlation, bias, and RMSE; the y axis shows the authors of papers, and the three columns show the metrics.  The "-" in the table shows that this item was not discussed in the reviewed paper or the paper did not apply this item. "S" character in column of "sensors" refers to different configurations that have been applied in a paper. The "-" in the table shows that this item was not discussed in the reviewed paper or the paper did not apply this item. "S" character in column of "sensors" refers to different configurations that have been applied in a paper. The "-" in the table shows that this item was not discussed in the reviewed paper or the paper did not apply this item. "S" character in column of "sensors" refers to different configurations that have been applied in a paper. The "-" in the table shows that this item was not discussed in the reviewed paper or the paper did not apply this item. "S" character in column of "sensors" refers to different configurations that have been applied in a paper. The "-" in the table shows that this item was not discussed in the reviewed paper or the paper did not apply this item. "S" character in column of "sensors" refers to different configurations that have been applied in a paper.  Table 6 displays the configurations of scenarios for RSM with correlation improvements greater than 0.09. These studies applied SMOS and AMSR-E SM data and 1D-VAR and PF assimilation methods. Table 7 shows RMSE improvements higher than −0.04 (m3 m-3). AMSR-E, SMOS, and TMI data were applied and EnKF, GPF, and evolutionary methods were used for conducting the assimilations. Table 8 shows the bias improvement in RSM for configurations with improvements of more than −0.07 (m3 m-3). They applied SMOS and AMSR-E SM and EnKF (most of them) and 1D-VAR as assimilation approaches.
Based on Figures 10 and 12 and experimental and physical conclusions in recent studies, the assimilation results show lower improvements in RSM. SM is highly dependent on the LSM structure (subsurface physics) and water balance modeling [151] along with the lag time (especially in wetter condition) that are required for the correction of RSM by assimilation [113].

Discussion
This study reviewed the studies that analyzed the improvement of LSM SM estimates by assimilation of satellite-derived SM and BTs. The details in assimilation and improvement of LSMs SM estimates through LSM and RTM calibration were presented and, finally, the results of reviewed studies were investigated. It should be noted that the assimilation results and the conclusions of the papers reviewed here could be spatially and temporally variable. It is necessary to analyze the outcome of high score assimilation configurations in analyzed studies in different spatial and temporal conditions. Based on the analyzed papers, along with the mentioned parameters (input data, bias correction, and assimilation details) in the last sections, the following considerations could be applied to gain better results in LSM SM correction in the future: a Joint Use of Calibration and Assimilation The calibration of LSM and RTM is a crucial approach, along with assimilation, for reducing biases. The calibration directly affects model climatology and the assimilation corrects the forcing data (precipitation) for driving the model. During calibration, model parameters (both fixed and time-variant) could be problematic, and the application of both calibration and assimilation may reduce such uncertainties. Zhao et al. [77] proposed pre-calibration of RTM parameters for gaining better results in both SSM and RSM estimation through assimilation. Tian et al. [26] applied a dual-pass assimilation system to simultaneously update model states and RTM parameters. Koster et al. [27] emphasized the effects of model calibration and assimilation as two separate methods for correcting model parameterization and errors induced by forcing data. They proposed that the joint use of these approaches could have complementary effects and result in more accurate surface SM estimates. However, it should be noted that such a simultaneous estimation is challenging when rapid changes are happening on vegetation opacity and soil surface roughness; this might reduce the compensatory role of parameters on each other [96]. b The Role of Precipitation, Vegetation, and Land Cover Recent studies show that assimilation results are dependent on precipitation forcing, vegetation, and land cover. The role of precipitation, wetness, temperature, and vegetation were investigated on assimilation results. Lin et al. [152] and Liu et al. [52] suggested that the assimilation of both SM and precipitation could provide better results. They stated that as both methods have independent corrections, their combination can have the largest impact on SM skills. Moreover, assimilation could be more beneficial in places where precipitation corrections are not possible [28,153].
As for vegetation cover, Draper et al. [154] showed that significant improvements were found over croplands and mixed cover for the root zone. Furthermore, Tian et al. [155] found better results over croplands, grasslands, and mixed cover for SSM. It was proposed that SM improvement in dry to moderate wetness conditions was higher in some studies [116,153]. Additionally, [156] reported that due to the presence of frozen soil (low temperature), high vegetation cover and noise in dry regions resulted in a reduction in the accuracy of satellite retrieval, and larger improvements could be observed in areas with moderate vegetation cover. c The Joint Assimilation of Satellite Data It is vital to use the best data; Gevaert et al. [153] stated that the joint assimilation of L-and C-band (or X-band) retrievals had higher results in contrast with individually assimilating C-band (or X-band) retrievals in AWRA-L LSM. Draper et al. [154] used ASCAT and AMSR-E SM products individually and in combination. They stated that when these two active and passive data are used together, better results are obtained. Yin et al. [68] used SMOPS SM products and showed that they had superior results compared to ASCATA, ASCATB, SMAP, SMOS, and AMSR2. However, Kolassa et al. [157] conducted assimilation experiments from active and passive SM products in CLSM LSM. They used the retrieved SM from both active (ASCAT) and passive (AMSR-E) platforms and found similarly improved results in the surface from both experiments in contrast with individual products. As shown in Tables 3-8 (based on available data), the highest improvements may be gained by the SMOS and AMSR-E sensors.
The application of other satellite products can also be a beneficial solution. Gavahi et al. [158] applied MODIS evaporation (MODIS16 ET) along with SMOPS SM product in an assimilation study and found improved results compared to a univariate assimilation procedure. According to the abovementioned studies, the joint use of satellite products may provide better results in the future. d

Selection of Assimilation Method
The selection of the assimilation method is highly dependent on the nature of the problem and observation data [159]. According to Moradkhani et al. [159], the following factors should be considered. (1) The linearity and non-linearity of the LSM and updating procedure. As the SM probability distribution in LSMs is realized as a non-linear process, it is expected that a non-linear assimilation assumption could provide more plausible results.
(2) Filtering and smoothing updating. It can be more beneficial in SM estimation to include backward observations simultaneously (4D-VAR). (3) The deterministic (KF, variation EKF) and ensemble-based assimilation procedures (EnKF and PF). The deterministic approaches could result in poor results in complex LSMs. (4) The computation cost. The required runtime in frequently used assimilation methods (Variational, EnKF, and PF) could be a challenge, especially for operational applications.
As shown in Section 3, most of the reviewed studies with the highest improvements used EnKF and variational approaches, respectively. It should be noted that in EnKF, the assumption of the non-Gaussian error structure in state variables and observations and the linear updating procedure can have negative impacts on assimilation performance [160][161][162][163]. Although variational methods are not bound to solve state-space modeling, they need the derivation of an adjoin model (linearized form of the model) to optimize the error minimization formulas. Moreover, the assumption of time-invariant model covariance could result in inherent errors in SM improvements [163]. The PF can be a solution when skewed model ensembles and non-Gaussian error structures are observed. In contrary to EnKF, where model state variables are updated at each time step, in PF, the probability distributions of state variables and observations evolve through time. Although PF was proposed to overcome the suboptimal problems of EnKF, it suffers from computational burden and an inherent problem called "particle degeneracy" when the weights of particles are not considerable, which causes inaccurate posteriori distribution estimates. To solve this problem, several methods were proposed (the details can be tracked in [162]).
As shown in Section 2.5.2, some studies compared the performance of assimilation methods regarding the improvement of SM in LSMs. The results from these studies may support the theoretical conclusions above that EKF (dealing with non-linearity) had better performance than OI [5] and that EnKF had better performance than EKF [138,139,164] and EnOI (it is cheaper but does not use the flow-dependent ensembles) [51]. EDA and PF (proposed IPF) had better results than EnKF [137,165].
Although the internal setting in the reviewed papers in this study are different and the number of studies using PF is low, in Table 3 (correlation difference in SSM), PF had better results than EnKF, but in Table 4 (RMSE difference in SSM), the PF had lower results than EnKF. This can be an issue for future studies in which evaluation metrics for PF have better performance than EnKF. e

Selection of Assimilation Period
Draper et al. [165] stated that periods longer than 31 days are essential for achieving better results from assimilation in the root zone; SM anomalies can take between 1 and 3 months [166,167]. Draper et al. [168] showed that the sub-seasonal scale had the largest impact on SM estimation. It is necessary that future studies investigate the role of the assimilation period while considering diverse spatial and temporal conditions. f Satellite BTs and SMs Kolassa, et al. [111] (2017) showed that, on average, using both BT and SM retrieval in assimilation has similar results (excluding local effects). De Lannoy et al. [169] showed that BT assimilation has slightly better results than SM products. They stated that this may show that BT could contain more information than the filtered information within SM retrievals. However, the results from our analysis in Section 3 showed no significantly greater improvements in BT with respect to the SM products. g The Spatial Resolution of Satellite Data According to Peng et al. [170], the application of high spatial resolution SM products (less than 1 km) is highly demanding in the hydrometeorological community. Promising projects such as the Sat'elite Argentino de Observacion' COn Microondas (SAOCOM) mission, the NASAISRO Synthetic Aperture Radar (NISAR), the Radar Observing System for Europe L (ROSE-L), and the Tandem-L satellites could provide such required resolutions, but at the present time, ESA Sentinel-1 European Radar Observatory can provide the highest resolution. Moreover, recently proposed downscaling algorithms could be applied to acquire better results in mid to low spatial scales using atmospheric and geophysical data from high spatial resolution remote-sensing data from other platforms [69,171].

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The full form of applied acronyms in this study are shown in the following: