Parameterization of Vegetation Scattering Albedo in the Tau-Omega Model for Soil Moisture Retrieval on Croplands

: An accurate radiative transfer model (RTM) is essential for the retrieval of soil moisture (SM) from microwave remote sensing data, such as the passive microwave measurements from the Soil Moisture Active Passive (SMAP) mission. This mission delivers soil moisture products based upon L-band brightness temperature data, via retrieval algorithms for surface and root-zone soil moisture, the latter is retrieved using data assimilation and model support. We found that the RTM based on the tau-omega ( 𝜏 - ω) model, can suffer from significant errors over croplands (in average between -9.4K and + 12.0K for Single Channel Algorithm SCA; -8K and + 9.7K for Dual-Channel Algorithm DCA) if the vegetation scattering albedo (omega) is treated as a constant and the temporal variations are not accounted. In order to reduce this uncertainty, we propose a time-varying parameterization of omega for the widely established zeroth order radiative transfer 𝜏 - ω model. The main assumption is that omega can be expressed by a functional relationship between vegetation optical depth (tau) and the Green Vegetation Fraction (GVF). The validation was performed from 14 May to 13 December 2015 over 61 Climate Reference Network sites (SCRN) classified as croplands. The application of the proposed time-varying vegetation scattering albedo results in a consistent improvement for the unbiased root mean square error of 16% for SCA and 15% for DCA. The reduction for positive and negative biases was 45% and 5% for SCA and 26% and 12% for DCA, respectively. This indicates that vegetation dynamics on croplands are better represented by a time-dynamic single scattering albedo.


Introduction
The prediction of weather extreme events, such as heat waves and cold surges, is important in time spans from one week to several months (S2S: sub-seasonal to seasonal) [1]. However, existing weather and climate models still perform very poorly in the prediction for this time scale. This issue is known as the weather and climate prediction gap [4]. At this time scale, the initial conditions of atmosphere, land and ocean components affect the prediction skill. One of the important missing pieces in S2S predictions is the role of the land surface; in particular, soil moisture, which is the main variable transferring water and energy to atmosphere. Furthermore, soil moisture plays an important role in cloud and precipitation formation emphasized in the recent modeling and land-atmosphere feedback studies [5][6][7]. Estimation of global soil moisture variability, particularly within the root zone, can be realized in a land surface model using data assimilation of remote sensing measurements [8]. Assimilation systems opens new possibilities to improve the accuracy and robustness of land surface models with microwave brightness temperature assimilated from satellite such as the SMAP mission [9][10][11] and SMOS mission [12][13][14]. For this purpose, accurate and realistic microwave radiative transfer modelling (RTM) is essential as an operator for simulating microwave brightness temperature (Tb). One of the uncertainty sources in microwave RTM is modelling of wave-canopy interaction, which is commonly represented with a zeroth-order RTM using vegetation optical depth (VOD) and single scattering albedo (omega) [15][16][17].
Currently, the SMAP baseline soil moisture algorithm (SCA, single channel algorithm) use an NDVI climatology-based VOD in its RTM [18]. The wavelength, or frequency, limits the penetration of electromagnetic waves through vegetation. At shorter the wavelengths there is less capacity it has to penetrate through the vegetation saturating the VOD at lower vegetation density. Longer wavelengths have the ability to capture VOD over wider range of the vegetation growth stages [19]. Therefore, low-frequency microwave measurements from L-band Radiometers such as SMAP and using algorithms such as Dual-Channel Algorithm, DCA [20] and Multi-Temporal DCA, MT-DCA [21] allows full penetration of wide variety of vegetation types.
In this study, we will investigate the improvement of the operational SMAP SCA and DCA algorithms by proposing a time-varying parameterization of omega for the two algorithms. Currently, both of the SMAP operational algorithms consider the scattering albedo as constant, e.g. value of 0.05, for cropland type, while in experimental algorithms such as the MT-DCA omega is varying in space but fixed over time domain of the retrieval period for SM and tau. An important difference between DCA and MT-DCA is whether omega is estimated by the cost function minimization along with tau and soil moisture. The main assumption in the minimization of the cost function of MT-DCA is that the temporal variability of scattering albedo is much larger than soil moisture and tau. However, the assumption of one fixed omega for each vegetation type domain may be invalid over heterogeneous surfaces and for fast growing crops. This heterogeneity issue ultimately adds to the uncertainty of soil moisture estimation using SCA, DCA and MT-DCA algorithms (e.g. [22]). In this study, we apply time-and space-varying omega that is synchronized with VOD and investigate whether a newly parameterized (time-varying vegetation scattering albedo) tau-omega ( -ω) radiative transfer model based on SCA and DCA is able to simulate Tb more accurately over spatially and temporally heterogeneous croplands.

The -ω model of vegetated soil emission
The -ω model represents a zeroth order solution of the radiative transfer equation [17] and is a common basis of current passive microwave electromagnetic interaction modeling with vegetated soils at L-band. This model is also applied in the SMAP soil moisture retrieval algorithms [18]. It expresses the aggregated brightness temperature in the resolution cell over of view as follows [23]: where, = exp(− /cos ) is the brightness temperature, emitted from land surface; is the soil emissivity; indicates the transmissivity of canopy which is determined by vegetation optical thickness at nadir incidence ; T is the physical surface temperature, and is the single-scattering albedo, set to a constant of 0.05 ( 0.05 ) for croplands in the SMAP SCA. In this study, this approach is called the fixed-omega approach. The basis to estimate the value of (or VOD) has arguably improved from the NDVI-based used in SCA. In DCA and MT-DCA, or and are directly determined from the polarimetric microwave L-band Tb, respectively. In this study, we focus on improving the scattering parameter , which is a constant in space and time for SCA and DCA and a constant in time for MT-DCA. In contrast, we by establish a spatially heterogeneous and temporally varying to account for the heterogeneity of vegetation scattering albedo in croplands and their dynamics. Owing to the varying omega, we differentiate this approach from the MT-DCA (multitemporal dual channel algorithm) where omega is a time-constant value over the optimization period. The latter is retrieved from a model selection during multi-temporal optimization of the estimation of and permittivity [24].

New Parameterization of vegetation scattering albedo with GVF
In order to derive the temporally varying vegetation scattering albedo (ω) within the -ω model, we assume that omega can be derived based on a proportionality to the sub-grid scale vegetation fraction, Green Vegetation Fraction (GVF) [25].
Based on this assumption, the temporal variability of is determined by the temporal variability of the vegetation fraction . With no vegetation scattering condition for the bare soil fraction (1-GVF), 0 becomes 0, which leads Eq. (2) to:

Combining tau and omega via GVF
In this study, our hypothesis is that we can parametrize the 2-D (spatial) vegetation cover fraction (GVF) with the measured VOD via a power-law function. Firstly, VOD (or ) can be expressed with a parameter b and the vegetation water content VWC [24], where b is a parameter related to the wavelength and vegetation structural characteristics. Now, we define the vegetation cover fraction with the vegetated area, A, per unit ground area.
Studies in the past have established empirical relationships between above ground biomass (AGB) and tree height, H. The allometric relationship has been derived as AGB ~ H 2 for forest by [26,27]. As vegetation grows, it typically increases in height (H) and covers a larger area (A). The height and area of vegetation can be related using allometric functions. Using allometric functions we express the 1-D height in terms of the 2-D area of the vegetation calculated with tree diameter D [28][29][30]. However, instead of using the ln(H)-ln(D) 2 non-linear approach, we apply an H-D linear approach without violation of their physical units as shown in Eq (6), where c is a non-dimensional factor related with environmental variables and model uncertainty from the proposed function [28]. In a recent study, total VWC in a SMAP grid was scrutinized in terms of volume and height of canopy by [31], 4 of 21 where physical density of plant elements ( ), density of canopy in plant elements ( ), volume of vegetation ( ), height of vegetation layer ( ). In this study, we express V as a function of a (area of a plant element) and h (the unique thickness the unique thickness of the plant element): If all plant elements are homogeneous in a measured resolution cell, we can compute the vegetation area as shown in Eq. (9) where is the number density of the plant elements. Then we can get the volume of a plan element from A and h.
This results in a new GVF [cm 2 /cm 2 or %] -[-] relationship. The new relationship is differentiated from the exponential function of LAI which can be estimated as tau or VWC via the approximated relation ( = 0.5*LAI) [32] to estimate the vegetation fraction proposed by [33]. Chaubell et. al and Fernandez-Moran et. al ( [34][35][36]) suggested that VOD is proportional to grass or crop height linearly. But the non-linearity between VOD -vegetation fraction turned out to be the power-law function with 2/3 exponent as shown in Eq. (13). Finally, without ancillary input, can be derived as power-law function of tau based on Eq. (3) and (10) as following.

Experimental Results and Validation of Parameterization
In order to confirm the developed time-dynamic vegetation scattering albedo approach, we performed a validation process. The control cases (SCA1 and DCA1) are used for Tb simulation with in-situ SM which is the reference input as shown in Fig.1. In this step, the difference between the simulated and observed Tb is considered as the modelling mismatch (mainly ω in this study).
The standard -ω is used for Tb simulations with in-situ SM as reference input. In this simulation, the difference between Tb simulated and the observed is considered as an error. With the same in-situ SM input, we simulate Tb but this time by applying the new parameterization of vegetation scattering albedo, . We evaluate the differences between the newly parameterized, time-varying -ω model (SCA2 & DCA2) with the results obtained using the control runs (SCA1 & DCA1). The amount of reduction (SCA2 -SCA1 and DCA2 -DCA1) represents the RTM improvement due to the timevariation of the vegetation scattering albedo, .

Data
In-situ soil moisture from the U.S. Surface Climate Observing Reference Networks (USCRN) soil moisture network [41] was used as the input for Tb simulations from May to November 2015, which are used as the reference for the comparisons. The USCRN sites and soil moisture networks selected for the investigation are presented in Fig. 9 in Appendix: A are located on croplands (with information of crop type) according to MODIS IGBP land cover classification. The detailed description of the study sites is provided in Table 1.
In the SCA ( ) case, the -ω model uses a value estimated from the MODIS NDVI data. In the DCA ( ) case, is retrieved simultaneously in addition to the SM. In both cases, omega is constant 0.05 for the crop surface type following [42].
For the newly parameterized approach, the Tb simulations consider the canopy interaction heterogeneity in the -ω model by applying time-and space-variable , which is a function of the estimated in SCA or DCA. The heterogeneity inclusion in the DCA and SCA will be investigated by comparing the SMAPL2 soil moisture product [43][44][45] in the specific crop sites over USCRN. Furthermore, the validation SMAP Level 2 Enhanced Passive Soil Moisture Product [18,46] will be performed from 2015 to 2019 presented in Table 2. The detailed description of the validation data with SMAP level 2 products at USCRN validation sites and SMAP Level2 Enhanced Products in core validation sites are provided in Table 6 and Table 7.

New parameterization of in the -ω model
The parameter required in Eq. (11) is determined from temporal average of and VIIRS GVF measurements over the calibration sites (TERENO, HOBE, REMEDHUS, RISMA) as shown in Fig. 2. The determined in this study is 1.12 for the GVF simulation (P-value from Wilcoxon rank sum test is 0.6817, which means our hypothesis is reliable enough). The computation of the time-varying based on Eq. (11) requires also the maximum . For the new parameterization of forward model parameters, the time-varying was tuned via the optimal gamma (1.12) and . The results of the calibration and validation are presented in Table  2.

Quality assessment of the new parameterization in the -ω model
We investigate if the -ω model error in the simulation of Tb is reduced by replacing the timeconstant omega ( 0.05 ) with time-varying omega ( var ) that depends on the value of . Eq. (11) indicates that a higher measured in a SMAP resolution cell is likely to have a higher effective value of omega; Higher → larger vegetation fraction (less bare soil) in a grid → higher effective Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 30 July 2020 doi:10.20944/preprints202007.0717.v1 Figure 3a shows the significant overestimation of SMAP that SM RTMs can produce. Particularly, the SCA based SM reached the limit value, up to 0.6 m 3 /m 3 during half of the time-series. These errors (Fig. 3b) were estimated by deducting in-situ SM and are temporally correlated with the varying (Fig. 3c). The SM estimation is affected by the required ancillary parameters of vegetation, and . If one of the ancillaries is not realistic -in this study the time-constant -it will affect the SM estimation. In other words, one of the error sources in SM are the vegetation properties within RTM and this error is at least to a certain extent addressed with time-and space-varying . This result confirms the validity of the hypothesis that this can be approximated with . The -derived was more than 0.1 and two times larger than the constant applied in SMAP baseline algorithm. We investigated the improvement by applying the time-and space-varying . The time series of Tb in (Fig. 3d) shows that the overestimated Tb (blue) decreases in Tb simulations with DMM (red curve). The effect of the new parameterization of in the -ω model is displayed in Fig. 3e. The application of varying significantly reduces the Tb bias form SCA1 to SCA2. Over cropland as presented in Fig.3 a, this type of SM bias seems to be more of a serious issue in SCA than DCA-based soil moisture retrievals. We can expect the unrealistically overestimated SM from both approaches will have a positive effect by applying varying during the SM estimation process from the measured Tb. The SMAP DCA SM estimates in the right panel of (3a) are close to 0.6 cm 3 /cm 3 missing the seasonal SM evolution observed in-situ. The time-varying ranges between 0.08 and 0.12, which is a much larger value than the default value, 0.05. In addition to the large difference in the absolute value, a temporally varying pattern that exhibits a similar pattern of the SM error due to the constant . In this case, the application of the time-varying also significantly reduced the overestimation of Tb.In On the other hand, in the Fig. 4, the overestimation of SM by using the constant in SMAP RTM is less severe than Fig. 3 showing the limited SM value in all-time series in SCA approach. But in this case SM estimated by SCA is much closer to the in-situ SM than the one by DCA. The SCA used in the computation of (c) in Fig. 4 are lower than the one in Fig. 3. Still, the DCA of Fig. 4 range from 0.8 to 0.12, which is similar to the Fig. 3. With given TB and higher , the SM is higher in the simultaneous optimal estimation. It means that the DCA in Fig. 4c should be lower. Particularly, DCA SM error becomes larger when was high in DOY 220-230 and 240 & 280, which leads to DCA-based SM much higher than the in-situ in this period. Probably, the further improvement of DCA approach for simultaneous estimation of and SM can be expected in the minimization process finding optimal with the temporally varying than the constant . As a result, SM estimated by SMAPL2 SCA and DCA was overestimated as shown in Fig 3 and 4 (a). The difference of the SMAP SM to the in-situ in (b) shows the temporal correlation with the changes of the omega in DCA of Fig. 3 and 4. It means that both SCA and DCA approaches suffer from a low value of ; DCA can detect the temporal changes of vegetation better, which is revealed in its SM error. The improvement in Tb simulation is mostly originated from the overall larger value of the new in both SCA and DCA and less because of the temporal variance. This uncertainty is attributed to the scattering properties of the -ω model which was the reason we replaced the constant with the varying . The results in Figs. 3 and 4 suggest that the soil moisture estimations using the -ω model based on the fixed were mostly underestimated and the new is on a higher level than the constant oneshowing reduced bias compared to the measured Tb. As a result, the vegetation variability in the newly parameterized -ω model improves the Tb simulation. Bias and ubRMSE tend to decrease. Owing to this, the SM estimation from the SMAP Tb will be closer to the in-situ SM. Furthermore, the newly parameterized -ω model provides a more accurate observation operator for data assimilation, which would result in more accurate soil moisture update to NWP. The validation over the crop sites matched with 9km Tb products, showed a very little improvement by varying (SCA1→SCA2 and DCA1→DCA2). Further case studies have been performed and the results are summarized in Table 4.  Overall, the biases were reduced (SCA1→SCA2 and DCA1→DCA2) and ubRMSE becomes closer to zero for croplands as shown in Table 5 and Fig. 7.  More details on the validation statistics for the sites used in Fig.8 can be found in Table 8 of the supplement results in the end of the manuscript.

Summary and discussions
In this study, we found that the soil moisture estimated with SCA and DCA from the SMAP mission suffers from over-and under-estimations for cropland sites. In order to tackle this bias, we derived a varying omega (ωvar) based on the assumption of a power law relationship between GVF and VOD instead of a time-constant omega (ω0.05) used in the SMAP baseline algorithms (SCA and DCA). The formulation allows us to express a time-varying omega ωvar based on the dynamics of . Hence, ωvar is able to account at least partly for the temporal variation of the vegetation properties in cropland. In this study, we focused on linking the measured VOD and the effective value of omega (effective single scattering albedo) mainly via vegetation volumetric traits such as the height and area fraction within the measured resolution cell.
The assessment was performed with the SMAPL2 brightness temperature (Tb). It is matched with the forward modelled brightness temperature using in-situ stations of the USCRN (27 croplands (11 corn, 7 soybean, 2 cotton, grapes, alfalfa, 1 wheat, citrus, unknown sites)) in 2015 and in the SMAPL2 Enhanced H-pol brightness core validation sites (3 cropland) from 2015 to 2019. As a result, we were able to reduce the positive Tb bias for several reference sites over cropland (C1, 3-11, 13-17, 20, 22, 25, 27 d in Table 8) including Gadsden-19-N (Atlanta) and Durham-2-N (Boston) presented in the figures 3 and 4. This bias reduction mitigates the overestimation of Tb (K) by 80% and 35% in the SCA and DCA approach, respectively. These results demonstrate that owing to a different phenology of the VOD time series over cropland, the time-varying omega parametrized with the VOD can implement more realistic -ω model than the one applied with the constant omega. In a future study, further experiments will be performed including organic matter (OM) to the applied dielectric mixing model [46]. The soil moisture of USCRN used for validation are also estimated from dielectric constant measurement [47], where the soil organic matter is not considered. We assume that this missing consideration in the reference soil moisture values might affect the validation of Tb simulations.
As satellite remote sensing is the only operational way to determine global soil moisture, an accurate radiative transfer model is essential. We propose that the presented parameterization for a time-varying vegetation scattering albedo from VOD dynamics implemented within the -ω model provides more realistic retrievals of soil moisture dynamics. The key feature of the approach is that no more variables are added with this new parameterization of the -ω model contributing to a more accurate but less complex global soil moisture estimation. This is equally important for retrieval and data assimilation approaches based on microwave brightness temperature measurements from SMOS and AMSR-2.    The crop type information was extracted from the 30m resolution Cropland Data Layer database [54] within SMAP's 36km grid boundaries for the year 2015 (*: both corn and soybean are dominated within the SMAP grid). Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 30 July 2020 doi:10.20944/preprints202007.0717.v1