Regularized Dual-Channel Algorithm for the Retrieval of Soil Moisture and Vegetation Optical Depth From SMAP Measurements

In August 2020, soil moisture active passive (SMAP) released a new version of its soil moisture and vegetation optical depth (VOD) retrieval products. In this article, we review the methodology followed by the SMAP regularized dual-channel retrieval algorithm. We show that the new implementation generates SM retrievals that not only satisfy the SMAP accuracy requirements, but also show a performance comparable to the single-channel algorithm that uses the V polarized brightness temperature. Due to a lack of in situ measurements we cannot evaluate the accuracy of the VOD. In this article, we show analyses with the intention of providing an understanding of the VOD product. We compare the VOD results with those from SMOS. We also study the relation of the SMAP VOD with two vegetation parameters: tree height and biomass.


I. INTRODUCTION
T HE soil moisture active passive (SMAP) mission was designed to acquire and combine L-band radar and radiometer measurements for the estimation of soil moisture (SM) with an average unbiased root-mean square error (ubRMSE) of no more than 0.04 m 3 /m 3 volumetric accuracy in the top 5 cm of soil for vegetation with water content of less than 5 kg/m 2 [1], [2]. SMAP released new versions of its SM and vegetation optical depth (VOD) products in August 2020 (version 4 for the L2/3_SM_P_E product and version 7 for the L2/3_SM_P).
The new implementation of the dual channel algorithm (DCA), which uses the two polarized brightness temperature measurements (H and V), generates SM retrievals, not only satisfying the SMAP accuracy requirements, but also showing a performance comparable to the single-channel algorithm that uses the V polarized brightness temperature (SCA-V) [2], [3], [4].
While the accuracy of the DCA SM can be evaluated by comparison with in situ data, the lack of VOD in situ data raises concerns about the accuracy of the VOD product. Although the SMAP mission does not have a requirement for the accuracy of the retrieved VOD, it is of great value for the SMAP team and the science community in general since it provides critical information about the water content of the aboveground biomass and its seasonal variations. In order to understand the performance of the VOD product, it is common to compare it to similar products from other missions and to look at it in relation to other vegetation parameters such as tree height and biomass.
Recent efforts to retrieve SM and VOD from SMAP L-band brightness temperature data have resulted in significant progress [5]. Konings et al. [5] uses a multitemporal DCA (MT-DCA), which assumes that VOD changes more slowly than SM and can be assumed to be almost constant between every two consecutive overpasses. In addition, the MT-DCA approach allows for the retrieval of a single temporally constant value of the scattering albedo per pixel. The soil moisture and ocean salinity (SMOS) mission [6], produces simultaneous retrievals of SM and VOD based on angular information in its V-and H-pol brightness temperature products. Its L-band VOD retrievals have been analyzed by several authors [7]- [9]. SMOS-IC [10] presented an alternative approach to the retrieval of the SM and VOD but is still using angular brightness temperature information.
In this article, we detail the methodology for the DCA implementation in Section II. In Section III we present results of retrieved SM over core validation sites (CVSs) [4], [19]- [30], and the sparse network (SP) [17], [18]. We also present VOD results in Section III. In section IV we compare the SMAP VOD with the vegetation parameters, including tree height and biomass. In Section V, we compare the SMAP VOD results with those obtained by SMOS.

II. REGULARIZED DUAL-CHANNEL ALGORITHM
The newly implemented DCA simultaneously retrieves the SM and VOD (τ θ ) by minimizing the cost function where TB sim V is the V-polarized simulated brightness temperature and TB sim H is the H-polarized simulated brightness temperature, λ is a regularization parameter, and τ * an initial expected VOD value.
To simulate the L-band emission of the soil-vegetation system the SMAP team uses the zero-order approximation of the radiative transfer equations, known as the τ -ω emission model [10]. The brightness temperature equation, which includes emission components from the soil and the overlying vegetation canopy, is given by where the subscript p refers to polarization (V or H), T s is the soil effective temperature, T c is the vegetation temperature, τ p is the nadir vegetation opacity, ω p is the vegetation effective scattering albedo, r p is the rough soil reflectivity, e p is the rough soil emissivity and θ is the incidence angle. The surface roughness reflectivity is modeled by where Q (polarization decoupling factor which is related to h by the linear relation Q = 0.1771 h), h, and N are the roughness parameters and r * p (θ) is the Fresnel reflectivity of the smooth surface where the index p and q (q opposite to p) account for the polarization V or H. The baseline SMAP implementation of the retrieval algorithms assumes that in (2) T s = T c at the early morning descending overpass and that ω p = ω and τ p = τ are polarization independent to reduce the number of algorithm parameters [2]. Note that in (1) τ θ is related to the nadir vegetation opacity by τ θ = τ sec θ.
The implementation of (1)-(3) requires that several parameters need to be assumed: a prior value of the scattering albedo based on land cover, roughness parameters as detailed in [3], effective soil temperature, clay fraction to determine the soil dielectric constant, as well as the regularization parameters τ * and λ [2].

A. Selection of Parameter λ
The retrieval of VOD through the DCA algorithm with λ = 0 (MDCA, modified DCA [3]) produces VOD results with high variability in the temporal dimension as shown in Fig. 1(a)(blue curve) and also in the spatial dimension. Fig. 1(a) displays an example of the MDCA climatology at the Monte Buey CVS. It also shows the seven-day average and the VOD based on Normalized Difference Vegetation Index (NDVI) climatology from the moderate-resolution imaging spectroradiometer (MODIS).  The selection of λ determines how much the residual noise [see Fig. 1(b)] will be suppressed and how much freedom the optimization algorithm will have to converge to a solution of VOD apart from the expected value τ * . Indeed, the value of λ = 0 means no regularization and very high values of λ will force the retrieval of VOD to converge to τ * .
To select potential candidates of λ we considered five years (April 4, 2015-March 31, 2020) of MDCA VOD data for each 9km EASE-2 (equal area scalable earth grid version 2) grid cell and computed its daily climatology, Fig. 1(a). To remove the seasonal variation, we computed a seven-day average [red curve in Fig. 1(a)] which is subtracted from the MDCA climatology resulting in the residual shown in Fig. 1(b). We used the standard deviation (σ) of the residual to determine the amount of regularization needed and defined λ as 1/ σ. Fig. 2 displays a global map of σ (top) and the histogram of the corresponding  λ (bottom). Fig. 2 shows that σ varies across the globe and that the histogram peaks at λ∼ = 20.
We evaluated the performance of the SM retrievals compared to in situ SM data from 15 SMAP CVS by setting λ = 20 fixed over all the sites and also using the local values of lambda over each CVS. Fig. 3 displays the values of lambda for each CVS. In what follows we will refer to RDCA for the regularized DCA algorithm with regularization parameter λ = 20 and RDCA-λ for the regularized DCA algorithm with λ changing spatially. Table I gives the acronyms used to differentiate the different DCA implementations. Table II gives the algorithm performance metrics for SM retrievals. The metrics are the average of those at each CVS. The metrics show that even though a nonzero lambda is crucial to get good performance by the DCA algorithm, the selection of the values of λ is not relevant as long as λ is allowed to vary locally about the value λ = 20. For this reason and considering the ease of the operational implementation we decided to set λ = 20 globally.

B. Selection of the Initial Guess τ *
As an initial guess τ * , an estimate of VOD based on NDVI climatology was selected [2]. This selection aroused concerns regarding how the seasonal behavior of the RDCA VOD would be affected.
To address this, we compared the Pearson correlation of the daily climatology of RDCA VOD with daily climatology of NDVI VOD [see Fig. 4(a)] and daily climatology of MDCA VOD [see Fig. 4(b)]. The figure shows that while the correlation between RDCA VOD and the NDVI VOD is moderate (mean value of 0.51 and standard deviation of 0.51), there is a high correlation globally between RDCA VOD and MDCA VOD (mean value of 0.76 and standard deviation of 0.25) which is an indication that RDCA VOD follows the MDCA VOD seasonal variation more consistently than the seasonal variation of the NDVI VOD. We also looked into the time difference (in days) between the NDVI τ * peak location (see Fig. 1, black curve) and the peak location of the RDCA VOD climatology (see Fig. 1, red curve) and similarly for MDCA VOD and RDCA VOD. Fig. 5 displays the histograms of differences. We can see that while the differences are widespread for NDVI τ * minus RDCA VOD, for RDCA VOD -minus MDCA VOD most of the values concentrate around zero. This is another indication that both MDCA VOD and RDCA VOD reach their maximum value at approximately the same time. These two tests suggest that in general the use of the NDVI τ * does not affect the seasonal variation of the MDCA VOD significantly.

III. ASSESSMENT
The SMAP mission validates the accuracy of the retrieved SM using several sources of information [15]. Among them are CVS which provide the ground-based data in a timely manner to the SMAP project, and SPs, such as the USDA soil climate analysis network [16], the NOAA Climate Research Network [17] and the Oklahoma Mesonet [18]. Table III gives how the SCA-V, MDCA, and RDCA SM retrievals compare at the CVS. We can see that a significant improvement has been reached by the implementation of the regularization term in the DCA algorithm. The table also shows that the RDCA and SCA-V are statistically similar.    [4].

IV. SMAP VOD VERSUS TREE HEIGHT AND BIOMASS
In this section, we analyze the correlation between the SMAP RDCA VOD and two vegetation parameters: tree height in unit of meters (m) [13] and the aboveground biomass density of vegetation in units of Mg/ha [14]. Both sets of data were aggregated to the 9km EASE-2 to match the enhanced SMAP data, Fig. 6. Fig. 7(a) displays the density plots of VOD versus tree height and Fig. 7 (b) and (c) displays VOD versus biomass. Fig. 7 also displays the mean values for several bins and the fitting curve. To compute the mean values, the data were binned by intervals of tree height of 1 m and the biomass by intervals of 10 Mg/ha.
There is clear linearity between SMAP RDCA VOD and tree height spatially for values of tree height less than 20 m and the relationship remains fairly linear up to tree heights of ∼ 35 m/ the slope of the fitting curve is 0.033 and the offset is -0.06. The spatial correlation between SMAP RDCA VOD and the tree height map is R = 0.81 (strong).
The VOD versus biomass density plot [see Fig.7(b)] also shows linearity for values of biomass less than 90 Mg/ha. The VOD versus biomass (B) fitting curve for values of B between 0 and 90 Mg/ha is given by vod = 0.004854 B + 0.2541.
For values of biomass greater than 90 Mg/ha and less than 150 Mg/ha the VOD stays almost constant, this could be caused Heterogeneity of the grid pixel could be one of the factors causing the scatter in the density plots. The spatial correlation between the SMAP VOD and the biomass map is R = 0.53 (moderate).

V. SMAP RDCA VOD VERSUS SMOS VOD
The lack of VOD in situ data makes it difficult to evaluate the accuracy and performance of the SMAP RDCA VOD product. In order to understand the performance of the VOD product we compare the SMAP RDCA VOD product (L2_SM_P_E v4 [12]) with the SMOS Level 3 VOD. 1 For comparison, the SMAP RDCA VOD product was multiplied by cos(40 o ) to match the SMOS VOD product at nadir. Fig. 8 displays global maps of VOD. The SMAP data were aggregated to the 25 km EASE-2 grid to match the SMOS data. We observed that the SMOS VOD has higher values of VOD. In fact, the mean VOD and standard deviation for the SMOS are 0.42 and 0.28 respectively, while for SMAP, those values are 0.29 and 0.26, respectively. This difference in value may be caused by different levels of the roughness parameters used by the two missions. for SMOS, SMAP, and the NDVI VOD by land cover type following the MODIS-based IGBP (International Geosphere Biosphere) classification. For the computation of the statistics, we considered only pixels with Gini-Simpson-Index (GSI) less than 0.1. The GSI is commonly used in ecology as a measure of degree of homogeneity, where GSI = 0 means total homogeneity and is computed as where f i is the fraction of the area covered by the ith land use classification and n is the number of land cover types. Table V gives that in general there is very good agreement between the magnitude of the SMAP RDCA VOD and the magnitude of the NDVI VOD. This is somehow expected due to the nature of the roughness parameter h implemented by the SMAP algorithm [3]. In the τ -ω emission model τ and h cannot be seen as parameter independent of each other and the magnitude of the selected h will affect the magnitude of the retrieved τ . Since the values of h are obtained by a DCA-type algorithm involving NDVI τ as an input, we expect to have retrieved values of τ of magnitude similar to the NDVI τ . However, since the values of h are temporal invariant, the seasonality variation of τ should not be affected.
Table V also gives that there is very good agreement with SMOS VOD over forested areas although SMOS data seem to have more variability. In fact, SMOS data seem to have more variability for all the land cover types except for Urban and builtup settings. For land cover types other than forest we observed significant discrepancies in the statistics between SMOS and SMAP RDCA VOD. Fig. 9 displays the VOD differences for two consecutive years: 2015-2016 and 2016-2017. SMAP RDCA VOD tracks yearly changes in VOD unlike the NDVI VOD which is a ten-year climatology. The SMOS VOD product exhibits more variability than the SMAP RDCA VOD. SMAP RDCA VOD is smoother due to the use of NDVI VOD as regularization parameter τ * . There are some similarities between SMOS and SMAP but also some discrepancies. For example, over Australia, in Fig. 9 (top row), the trends seem to agree although SMOS shows greater differences. On the other hand, in Fig. 9(bottom row) over the same region the trends are distinctly different except for a portion in the east-central part of the country. We also see discrepancies in Fig. 9 (top row) for the east coast of the United States. Fig. 10 displays the Pearson correlation (R) between SMAP and SMOS VOD. The aggregated SMAP RDCA VOD and the SMOS VOD were averaged monthly, thus obtaining two data sets of dimensions (1388584,60) and then for each grid cell the temporal correlation was obtained. We observed that the correlation varies along the globe with a mean value of 0.344 (weak correlation) and standard deviation of 0.33. If we only consider correlation with p-values < 0.05 then the mean correlation value is 0.542 (moderate correlation) and the standard deviation is 0.26. It is noticeable that the correlation is mostly positive indicating a degree of agreement in trends.   Table VI gives the correlation by land cover types using only pixels with significant correlation, p-value < 0.05. We observed mostly moderate correlation, with the exception of permanent wetlands (PW) where the correlation is strong and evergreen needle leaf forest (ENF), evergreen broadleaf forest (EBF), closed shrublands (CS) and barren/sparse (BS) where the correlation is weak. Fig. 11 displays the monthly average of VOD at five different regions. From top to bottom the regions are as follows.      Deciduous broadleaf forest. The SMOS-SMAP correlation is strong R = 0.747 and p-value < 0.05. We see that in all the cases the SMOS and SMAP have consistent trends while NDVI VOD trends only agree with SMAP and SMOS over Chaco and South Fork. There is a big difference in the VOD magnitude over the Peruvian Amazonia and very weak correlation caused by the low seasonal change of VOD combined with the differences in short term variability. The correlation is not significant according to the p-value (we consider a correlation to be significant if the p-value is < 0.05). We also observed that the SMOS VOD has more variability which may be the cause of low correlation in some locations, as can be seen in the Zambia case, Fig. 11 second from the bottom.
The spatial correlation between SMOS VOD and SMAP RDCA VOD as shown in Fig. 8 is R = 0.83 (strong).

VI. CONCLUSION
In this article, we have shown that the regularized DCA algorithm (RDCA along this article) implemented in the new release (R17) allows for an accurate retrieval of SM and a reliable VOD (τ ). Indeed, we showed that the DCA SM not only satisfies the SMAP requirements but also showed accuracy levels comparable to the SMAP SCA-V baseline. We compared the SMAP RDCA VOD with the SMOS VOD. We showed that even though there are differences in magnitude, they have, in general, consistent temporal behavior tracking seasonal changes and strong spatial correlation.
Comparison of the SMAP RDCA VOD with tree height showed strong correlation and a linear relation especially for tree height less than 20 m.
Comparison of the SMAP RDCA VOD with vegetation biomass showed moderate correlation. We also observed linear correlation for biomass less than 90 Mg/ha.
The magnitude of the SMAP RDCA VOD is comparable to the magnitude of the NDVI VOD due to the nature of the selection of the roughness parameter h. However, since the values of h are temporal invariant, the seasonality variation of the retrieved VOD should not be affected.
The application of temporal variant roughness parameter h should be explored for further improvement of the retrieved VOD.

ACKNOWLEDGMENT
This article is carried out by the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration. The SMOS Data were obtained from the "Centre Aval de Traitement des Données SMOS" (CATDS), operated for the "Centre National d'Etudes Spatiales" (CNES, France) by IFREMER (Brest, France). USDA is an equal opportunity provider and employer. The research described in this publication also included measurements from the Long-Term Agroecosystem Research (LTAR) network. LTAR is supported by the U.S. Department of Agriculture. The University of Salamanca team involvement in this study was supported by the Spanish Ministry of Science and Innovation (projects ESP2017-89463-C3-3-R and PID2020-114623RB-C33), the Castilla y León Government (projects SA112P20 and CLU-2018-04) and the European Regional Development Fund (ERDF). The authors would like to thank Heather McNairn of the Agriculture and Agri-Food Canada (AAFC) for her work managing the core validation sites Carman. He is currently a Research Hydrologist with the Southeast Watershed Research Laboratory, Agricultural Research Service, Tifton, Georgia, USA. In 1986, he joined the Agricultural Research Service. He leads a watershed research program investigating the impacts of land use on water balance and quality. His research interests include watershed and landscape scale hydrology; agricultural impacts on water quality; hydrologic and solute transport modeling of watershed processes; riparian buffer hydrology and solute transport; and developing new methods for assessing the impact of agricultural chemicals on ground and surface water supplies.
Dr. Bosch has been active in the validation of remotely sensed soil moisture products since 2000. Todd Caldwell received the B.S. degree in geosciences from the University of New Mexico, Albuquerque, NM, USA, in 1997, and the M.S. and Ph.D. degrees in hydrogeology from the University of Nevada, Reno, NV, USA, in 1999 and 2011, respectively.
He is currently a Research Hydrologist with the U.S. Geological Survey in the Nevada Water Science Center. His research interests include in situ soil moisture measurements and scaling, near-surface geophysics, and numerical modeling and parameter estimation of soil hydraulic properties.