Multiscale Assimilation of Sentinel and Landsat Data for Soil Moisture and Leaf Area Index Predictions Using an Ensemble-Kalman-Filter-Based Assimilation Approach in a Heterogeneous Ecosystem

: Data assimilation techniques allow researchers to optimally merge remote sensing observations in ecohydrological models, guiding them for improving land surface ﬂuxes predictions. Presently, freely available remote sensing products, such as those of Sentinel 1 radar, Landsat 8 sensors, and Sentinel 2 sensors, allow the monitoring of land surface variables (e.g., radar backscatter for soil moisture and the normalized difference vegetation index (NDVI) and for leaf area index (LAI)) at unprecedentedly high spatial and time resolutions, appropriate for heterogeneous ecosystems, typical of semiarid ecosystems characterized by contrasting vegetation components (grass and trees) competing for water use. A multiscale assimilation approach that assimilates radar backscatter and grass and tree NDVI in a coupled vegetation dynamic–land surface model is proposed. It is based on the ensemble Kalman ﬁlter (EnKF), and it is not limited to assimilating remote sensing data for model predictions, but it uses assimilated data for dynamically updating key model parameters (the ENKFdc approach), including saturated hydraulic conductivity and grass and tree maintenance respiration coefﬁcients, which are highly sensitive parameters of soil–water balance and biomass budget models, respectively. The proposed EnKFdc assimilation approach facilitated good predictions of soil moisture, grass, and tree LAI in a heterogeneous ecosystem in Sardinia for a 3-year period with contrasting hydrometeorological (dry vs. wet) conditions. Contrary to the EnKF-based approach, the proposed EnKFdc approach performed well for the full range of hydrometeorological conditions and parameters, even assuming extremely biased model conditions with very high or low parameter values compared with the calibrated (“true”) values. The EnKFdc approach is crucial for soil moisture and LAI predictions in winter and spring, key seasons for water resources management in Mediterranean water-limited ecosystems. The use of ENKFdc also enabled us to predict evapotranspiration and carbon ﬂux well, with errors of less than 4% and 15%, respectively; such results were obtained even with extremely biased initial model conditions.


Introduction
Recent improvements in satellite remote sensing techniques obtain land information for ecohydrological studies at unprecedented fine spatial and time resolutions [1][2][3], which are also suitable for heterogeneous ecosystems, typical of semiarid and arid climates [4,5]. In these ecosystems, trees, grass, and bare soil components coexist, and grass and tree covers vary in time and almost randomly in space, at scales of the size of tree clumps (usually less than 50 m [6,7]). For these ecosystems, there is a need for satellite data at fine spatial and temporal resolutions, such as those provided nowadays by remote sensors such as Sentinel and Landsat 8 [8,9], which are also freely available. 2.1.1. Field Data A 10 m micrometeorological station was operating at the site to measure land-atmosphere fluxes of energy, water, and carbon in addition to key state variables. The apparatus included a Campbell Scientific CSAT-3 sonic anemometer and a Licor-7500 CO2/H2O infrared gas analyzer positioned adjacent to each other at the top of the tower. These two instruments measured velocity, temperature, and gas concentrations for the estimation of

. Field Data
A 10 m micrometeorological station was operating at the site to measure land-atmosphere fluxes of energy, water, and carbon in addition to key state variables. The apparatus included a Campbell Scientific CSAT-3 sonic anemometer and a Licor-7500 CO 2 /H 2 O infrared gas analyzer positioned adjacent to each other at the top of the tower. These two instruments measured velocity, temperature, and gas concentrations for the estimation of sensible heat flux, evapotranspiration (ET), and CO 2 exchanges (F c ) with the standard eddy covariance method (e.g., [59]). Half-hourly statistics were computed.
The two-dimensional footprint model of Detto et al. [6], previously tested for this site, was used for interpreting eddy correlation measurements in the context of the contributing land cover area. The combined use of the footprint model and the satellite images allowed us to interpret the eddy-covariance-observed surface flux and distinguish the source area of each vegetation component and bare soil from the measured flux, using the methodology in Detto et al. [6]. Considering the main eddy covariance tower footprint at the Sardinian site [6], an area of~500 m × 500 m around the tower (~90% of the footprint) was considered for model predictions ( Figure 1).
The surface temperature of the tree canopy and grass/bare soil patches (using IRTS-P by Apogee Instruments), incoming and outgoing shortwave and longwave radiation (from which to derive net radiation) (using CNR-1 by Kipp and Zonen), and the soil heat flux and temperature (using HFT3 REBS) at two locations close to the eddy covariance tower were monitored and half-hourly means were recorded. Precipitation was measured using a PMB2 CAE rain gauge. Seven frequency domain reflectometer probes (FDR, Campbell Scientific Model CS-616) were inserted in the soil close to the tower (3.3-5.5 m away) to estimate moisture at half-hourly intervals in the thin soil layer. Complete details of these measurements and data processing are available from Detto et al. [6], Montaldo et al. [14], and Montaldo et al. [58].
The LAI was measured indirectly through a ceptometer (Accupar model PAR-80, Decagon Devices Inc., Washington, DC, USA), which measures the PAR in the 400-700nm waveband, and estimates the LAI from these readings (details are given in the instruction manual edited by Decagon Devices Inc.). LAI measurements were performed mainly during the grass growth season [14]. Finally, specific leaf areas (LAI divided by dry biomass) of predominant grass (0.01 m 2 gDM −1 ) and woody vegetation (0.005 m 2 gDM −1 ) species were measured directly by weighing the dry biomass.

Remote Sensing Data
The Sentinel 1 radar data originated from the S1A and S1B satellites, and the level-1 ground-range-detected (GRD) data were used. The images were calibrated, noise corrected with a Lee filter (7 × 7), and resampled from 10 to 30 m spatial resolution [33]. S1A images were available from January 2015. From September 2016, S1B satellite images were also available.
Images from the Sentinel 2 radiometer were acquired at the L1C level and atmospherically corrected with the Sen2Cor tool of the Sentinel Application Platform (SNAP), or directly at the L2A level (already corrected). For Landsat 8, the L1TP product was used (it is radiometrically calibrated and orthorectified using ground control points and a digital elevation model), and the dark object subtraction (DOS) method was used for the atmospheric correction.

The Proposed Assimilation Approach
The assimilation approach includes: (1) remote sensing data, (2) the ecohydrological model, (3) the EnKF, and (4) the updating procedure of key LSM and VDM parameters ( Figure 2). Below, each component of the proposed assimilation approach is described.

The Proposed Assimilation Approach
The assimilation approach includes: (1) remote sensing data, (2) the ecohydrological model, (3) the EnKF, and (4) the updating procedure of key LSM and VDM parameters ( Figure 2). Below, each component of the proposed assimilation approach is described. Figure 2. The scheme of the multiscale assimilation approach: soil moisture (θ) and leaf area index (LAI) are the two assimilated and predicted variables by the coupled land surface model (LSM) and vegetation dynamic model (VDM); the saturated hydraulic conductivity (ks) and the maintenance respiration coefficient (ma) are the dynamically updated parameters through the assimilation approach; EnKF is the ensemble Kalman filter applied to both LSM and VDM. The timescales of the Figure 2. The scheme of the multiscale assimilation approach: soil moisture (θ) and leaf area index (LAI) are the two assimilated and predicted variables by the coupled land surface model (LSM) and vegetation dynamic model (VDM); the saturated hydraulic conductivity (k s ) and the maintenance respiration coefficient (m a ) are the dynamically updated parameters through the assimilation approach; EnKF is the ensemble Kalman filter applied to both LSM and VDM. The timescales of the models, ENKF, and parameter updating are reported. The equations for the θ, LAI, k s , and m a adjustments are referenced within parenthesis.

Optical Remote Sensing Data for LAI Estimate
NDVI is estimated from red and near-infrared spectral reflectance measurements of satellite remote sensors. We used mainly the Operational Land Imager (OLI) Landsat 8 data, and, secondly, the Sentinel 2 radiometer data for increasing the database. Landsat 8 has a temporal resolution of 16 days and a spatial resolution of 30 m for the optical bands, which is coarser than the resolution of Sentinel 2 data (up to 10 m; temporal resolution: up to 5 days), and both are freely available.
From Landsat 8 and Sentinel 2 data, NDVI is estimated at a 30 m spatial resolution. LAI is related to NDVI through the Γ N operator using an empirical approach (e.g., [18,60,61]): where β 1 , β 2 , and β 3 are coefficients for vegetation species. Note that analogous solutions should be derivable from different Γ L (NDVI) relationships. The NDVI map was also used for identifying the f v,t fraction of tree cover in the field, which was estimated as NDVI/NDVI max following Detto et al. [6], with NDVI max as the maximum value of NDVI in the investigated field.

Radar Images for Soil Moisture Retrieval
The dielectric constant of the surface soil, which is related to soil moisture, can be detected from radar data. We used the Sentinel 1 radar in VV polarization because VV polarization is more sensitive to soil moisture [28,62] and less sensitive to vegetation compared with VH polarization [33,63]. The images were calibrated, noise-corrected, and Remote Sens. 2022, 14, 3458 6 of 27 resampled at a 30 m spatial resolution. The temporal resolution was approximately six days. Sentinel 1 images provide backscatter data in VV polarization (σ 0 vv ). We used the semi-empirical Dubois et al. [34] model for relating the ε dielectric constant to the backscatter signal. The Dubois et al. [34] model accounts for co-polarized backscatter only and was formulated using scatterometer data collected at six frequencies between 2.5 and 11 GHz. The validity range is ks < 2.5 (ks is the normalized RMS surface roughness), and the incidence angles were greater than 30 • and were restricted to co-polarization (VV or HH). Here, only the VV empirical relationship was used: where λ is the wavelength; k is the wave number equal to 2 π/λ; σ is the surface roughness; and β is the local incident angle related to the radar beam angle and the latitude, exposition, and slope of the site. The inversion of Equation (2) estimates the dielectric constant from the VV polarized backscatter coefficient, knowing the soil surface roughness and specific radar configuration parameters (wavelength and incidence angle). While the radar configuration parameters are known, σ is undetermined. We used the approach of Montaldo et al. [33], which relates σ with NDVI through an empirical σ (NDVI) relationship, for accounting for the grass growth effect on radar signal, and which was just estimated at the Orroli field site ( Figure 4b of Montaldo et al. [33]): With the objective of simplifying the parameterization of the retrieval models, the use of the σ (NDVI) of Montaldo et al. [33] allowed us to integrate vegetation effects in the roughness parameter, as previously suggested by Capodici et al. [64]. The use of Equation (3) in Equation (2) estimates ε, which is finally related to θ soil moisture through the Γ θ operator [65]:

The Ecohydrological Model
The ecohydrological model is a three-component coupled land surface-vegetation dynamic model (LSM-VDM). The LSM predicts the soil water and energy balances. The VDM estimates the LAI evolution through time for two vegetation components (grass and trees), which are used by the LSM for computations of the energy exchanges between soil and vegetation. The details are given in Montaldo et al. [14] and Montaldo et al. [57]. Here, a summary of the main components is provided.

The Land Surface Model
The LSM predicts the dynamics of water and energy fluxes at the land surface on a half-hour time step (Figure 2). It includes three components of the land surface: bare soil, grass, and trees, representing two vegetated components. It is derived from the LSM of Montaldo and Albertson [66] and the surface temperature states are estimated through the force restore method [67]. The root zone supplies the bare soil and vegetation with soil moisture for evapotranspiration and controls the infiltration and runoff mechanisms. The base of the root zone represents the lower boundary of the LSM. Equations for surface temperature and the components of the energy balance (sensible heat flux, soil heat flux, and net radiation) are applied separately for each land cover component, so that the model predicts the energy balance distinctly for each land cover component [14] ( Table 1).
The soil-water balance equation of the root zone is computed by ∂θ ∂t where θ is the soil moisture; d rz is the root zone depth; I bs is the infiltration rate on bare soil; I t and I gr are the throughfall rates infiltrating into the soil covered by trees and grass, respectively; q D is the rate of drainage out of the bottom of the root zone; E bs is the rate of bare soil evaporation; E t and E g are the rates of transpiration of trees and grass, respectively; f v,g is the fraction of grass cover; and f bs is the fraction of bare soil [14]. The throughfall rate is modeled through a balance equation of the intercepted water by the canopy reservoir (a function of the LAI), which produces throughfall when the reservoir is saturated [66,67]. An infiltration excess mechanism, based on the Philip's infiltration equation [68], was used for the infiltration. In unsaturated soil, the Clapp and Hornberger [69] relationships were used to describe the nonlinear dependencies of volumetric soil moisture and hydraulic conductivity on the matric potential. The q D rate was estimated using the unit head gradient assumption (Table 1; [11,14]). E t and E g were estimated distinctly using the Penman-Monteith equation (e.g., [70], p. 224) for each vegetation component. Canopy resistances, accounting for environmental stresses, were estimated using a typical Jarvis [71] approach ( Table 1). The actual rate of bare soil evaporation was determined as α(θ)PE, where α(θ) is a rate-limiting function, estimated by the polynomial function of Parlange et al. [72], and PE is the potential evaporation estimated by the Penman equation (e.g., [70], Equations 10.15,10.16,and 10.19). Hence, the total evapotranspiration was estimated as: Paralleling the approach for ET estimation, a three-component approach was implemented for estimating the total net CO 2 flux [57]: where F c,t and F c,g are the carbon exchange of trees and grass, respectively, and R bs is the soil respiration. Carbon exchange rates for each PFT (i.e., F c,t , F c,g ) were computed as the difference between photosynthesis and growth respiration ( Table 2). Soil respiration was estimated as a function of the temperature (Table 2, [57,[73][74][75]). The model parameters are presented in Table 3. Table 1. Equations of drainage (q D ), canopy resistance (r c ) with stress functions of soil moisture (θ), air temperature (T a ) and vapor pressure deficit (VPD), sensible heat flux (H), net radiation (R n ), soil heat flux (G), and surface temperature (T s ) in the LSM. Parameters are defined in Table 3.

Equations Source
Drainage Canopy resistance r c = for T a ≤ T a,min and T a > T a,max 1 − T a,opt −T a T a,opt −T a,min for T a,min < T a < T a,opt 1 for T a,opt ≤ T a ≤ T a,max f 3 = 1 − ω log(VPD) [74] Sensible heat flux H = ρ a c p C H u(T s − T a ), where C H the heat transfer coefficient [14] Remote Sens. 2022, 14, 3458 8 of 27

Ecophysiological Term Equations Source
Photosynthesis Allocation For the tree cover: ξ a + ξ s + ξ r = 1; λ = e −k e LAI [14] For grass cover: ξ a + ξ r = 1 [14] Respiration Maintenance and growth respirations of biomass components: R l,µ = m a f 4 (T)B l ; R l,γ = g a a a P g R s,µ = m s f 4 (T)B s ; R s,γ = g s a s P g R r,µ = m r f 4 (T)B r ; R r,γ = g r a r P g [77] f 4 (T) = Q Tm 10 10 with T m = mean daily temperature [76] Soil respirationR bs = R 10 Q Tm 10 N where R 10 is the reference respiration rate at 10 • C and Q N is the soil respiration sensitivity to temperature [57] Senescence S l = d a B l S s = d s B s S r = d r B r [14,77] Litterfall L a = k a B d [14,77]

The Vegetation Dynamic Model
The VDM computes the change in biomass over time from the difference between the rates of biomass production (photosynthesis) and loss, mainly respiration (e.g., [10,75]). The VDM distinguishes tree and grass components and was adapted from Montaldo et al. [76], who derived a VDM for grass species starting from the Nouvellon et al. [77] model. In the VDM of trees, four separate biomass states were considered: green leaves (B l ), stem (B s ), living root (B r ), and standing dead (B d ). The VDM of grass only distinguishes three biomass states (green leaves, roots, and standing dead). The biomass [g DM m −2 ] components were simulated using the approach of Montaldo et al. [76], which consists of a balance between biomass production (related to photosynthesis for green leaves, stem, and roots biomass) and biomass destruction (respiration and senescence for green leaves, stem, and roots biomass), through ordinary differential equations, integrated numerically at a daily time step.
where Ph is the gross photosynthesis; a a , a s and a r are the allocation (partitioning) coefficients for leaves, stem, and root states, respectively; R l,µ , R s,µ , and R r,µ are the maintenance respiration rates from leaves, stem, and root biomass, respectively; R l,γ , R s,γ , and R r,γ are the growth respiration rates from leaves, stem, and root biomass, respectively; S g , S s , and S r are the senescence rates of leaves, stem, and root biomass, respectively; L a is the litterfall. The model equations are given in Table 2 and the parameters are presented in Table 3. The leaf area index was estimated from the biomass by a linear relationship [13,[76][77][78]: where c l is the specific leaf area of the green biomass. The VDM provides estimates of daily values of leaf biomass and, thus, the LAI of the tree and grass was used by the LSM to estimate evapotranspiration, energy flux, rainfall interception, carbon assimilation, and the soil water content at a half-hour time step [14]. The LSM provides soil moisture and aerodynamic resistances to the VDM. The coupled model was calibrated and validated in Montaldo et al. [14] and Montaldo et al. [57]. The details are given in Montaldo et al. [76], Montaldo et al. [14], and Montaldo et al. [57].

The Ensemble Kalman Filter
Using the EnKF, we assimilated observations of 1) NDVI, which is related to LAI through Equation (1) in the VDM that describes the evolution of the leaf biomass through (8) that is related to LAI through (12); and 2) the radar-derived dielectric constant of the ground, which is related to θ through (4), in the LSM that describes the evolution of θ using (5). Hence, in the proposed approach, the EnKF is applied distinctly to the LSM and the VDM of both grass and tree components ( Figure 2). In general, in the Kalman filters, → ϕ is a vector of surface state variables (i.e., θ or LAI). The equation describing the evolution of these states (Equation (5) for θ and Equations (8) and (12) for LAI), as determined by a nonlinear model ( → f ), can be written in general as (e.g., [79]): where → ω relates errors in model physics, parameterization, and/or forcing data, and is taken to be with mean zero and covariance → Ω.
→ H is the operator that represents the observation process that relates → ϕ to the → δ t j actual measurements available at the time, t j , as follows: where → represents the vector of measurement errors, assuming a probabilistic distribution with a zero mean and covariance → R. In the EnKF [39,40], an ensemble of ϕ ζ (ζ= 1, . . . , N e , with N e the size of the ensemble) is predicted in parallel, using (13). The EnKF updates each ensemble member separately, using the → δ t j observation and the diagnosed state error covariance → P − t j (e.g., [39], Equation (6)). The superscripts "−" and "+" refer to the state estimates before and after the update at time t j , respectively. Ensemble members are updated using a combination of forecast model states and the observations [39], as follows: where → K is the Kalman gain, which depends on ζ is a random realization of the measurement error, which should have the same statistical properties as the error included in (14) [80]. The mean of the ensemble members, ϕ + t j , is the state estimate of the variables (i.e., θ + or LAI + ).
Model errors in the EnKF are included through errors in the model initial conditions, physical parameters, and forcing data. In the assimilation of radar backscatter data in the LSM we included errors in (1) soil moisture initial conditions, (2) precipitation (whose uncertainty is expected to have significant impacts on the distribution of soil moisture, e.g., [80]), and (3) a key parameter-the saturated hydraulic conductivity (k s )-following Montaldo et al. [38]. The ensemble of soil moisture initial values is generated by altering a particular value of soil moisture through the addition of a normally distributed perturbation with a zero mean and SD θ (standard deviation). At each time step, the ensemble of precipitation is generated by multiplying the recorded precipitation value by a normally distributed random variable. An ensemble of saturated hydraulic conductivity values (k ζ s ) is generated as being log 10 normally distributed with a mean of log(k s ) (indicating withk s the base (i.e., best guess) value of the k ζ s ensemble) and the standard deviation of SD logks . In this way, an ensemble of θ ζ , which includes model errors, is generated and evolves in time according to (5).
We also assimilated grass and tree NDVI data in the VDM, including errors of (i) LAI initial conditions; (ii) incoming short-wave solar radiation (R swin ), and the photosynthetically active radiation (PAR); and (iii) model parameters-the maintenance respiration coefficients for the aboveground biomass (m a ) of grass and trees ( Table 2). We chose m a as the VDM parameter for data assimilation after a sensitivity analysis of LAI to VDM parameters, which proved the high sensitivity of grass and trees LAI to m a [56]. The ensemble of LAI initial values was generated by altering a particular value of LAI through the addition of a normally distributed perturbation with a zero mean and SD LAI standard deviation. At each time step, the ensembles of R swin and PAR were generated by multiplying the recorded R swin and PAR values by normally distributed random variables. The ensembles of grass and tree maintenance respiration coefficients (m ζ a,g for grass and m ζ a,t for trees), were generated as being normally distributed with means ofm a,g andm a,t and standard deviations of SD mag and SD mat , respectively. In this way, ensembles of LAI ζ of grass and trees, which include model errors, were generated and evolved in time according to (8) and (12). The time steps of models and observations are shown in Figure 2. The radar → δ observations available at time t j,θ were obtained including the → l random error in the ε observations derived from Sentinel 1 according to (14), where the operator → H is the inverse of Γ θ in (4). Similarly, the NDVI observations available at time t j,L derived from Landsat 8 (or Sentinel 2) were altered randomly according to (14), where the operator → H is the inverse of Γ L in (1).
When observations from Sentinel 1 are available, the ensemble of θ ζ (i.e., θ ζ− t j,θ ) is replaced by (e.g., updated to) the ensemble θ ζ+ t j,θ that is optimally estimated by (15) using the radar backscatter observations. When observations from Landsat 8 (or Sentinel 2) are available, the ensemble of LAI ζ (i.e., LAI ζ− t j,L ) is replaced by (or updated to) the ensemble LAI ζ+ t j,L , which is optimally estimated by (15) using the NDVI observations.

The Updating of Model Parameters through the Assimilation
The EnKF approach compensates for both inaccurate initial conditions and moderate model parameter errors. In presence of high inaccuracy of model parameters, they can be adjusted dynamically through the assimilation process. The assimilation procedure includes an update of the k s parameter of the LSM, and the grass and tree maintenance respiration coefficients of the VDM.
The k s parameter is updated using the approach of Montaldo et al. [38], based on an expression derived by Montaldo and Albertson [81], that estimates the biased error in k s from analysis of the persistent-state variable bias (as defined by a longer time average). Each component of the k ζ s ensemble is updated over an appropriate averaging time interval (∆t 5 ; Figure 2), which coincides with time steps t ξ,θ , through where ∆t 3 is the radar observation time steps (Figure 2). The overbar in (17) and (18) provides an averaging in the ∆t 5 time steps (≥∆t 3 ) for capturing an estimate of the "persistent" moisture bias estimating the required change in the saturated conductivity. In this way, the biased model error can be removed after a learning (calibration) period, and the Kalman filter assumption, the zero mean model error, can be recovered. The maintenance respiration coefficient for the aboveground biomass of the VDM is updated using the approach of Montaldo and Gaspa [56] (Appendix A), which updates (i.e., dynamically adjusts) the m a based on observations of persistent bias in the modeled biomass (i.e., LAI). The proposed procedure derives the required m a adjustment from the conservation equation of the biomass (i.e., LAI), and we applied it for both grass and tree LAI. Each component of the m ζ a,g and m ζ a,t ensembles was updated over the ∆t 6 time interval, which coincides with time steps t ξ,L (Appendix A), as follows: where ∆t 4 is the NDVI observation time steps (Figure 2), and the overbar in (20) and (21) provides an averaging in the ∆t 6 time steps (≥∆t 4 ). In this way, an estimate of the "persistent" LAI bias is used for evaluating the necessary change in m a . Thereby, after a learning (calibration) period, the error of the model can be eliminated. We used the same solution for grass (m ζ+ a,g ) and tree (m ζ+ a,t ) maintenance respiration coefficients.

The Multiscale Assimilation Approach
In summary, the multiscale assimilation scheme includes the following elements ( Figure 2):

1.
A land surface model that predicts the ensemble of soil moisture states through (5) at the half-hourly timescale (∆t 1 ); 2.
A vegetation dynamic model that predicts the ensembles of grass and tree LAI through (8) and (12) at a daily timescale (∆t 2 ); 3.
EnKF filters of the ε observations (4), which are available every 6 days on average (∆t 3 ); these account for moderate LSM errors and provide optimal updates of the ensemble of θ ζ− t j through (15) to arrive at θ ζ+ ; 4.
EnKF filters of the NDVI remote data (1) of grass and trees, available over the weekly timescale on average (∆t 4 ), which optimally update the ensembles of LAI ζ− of grass and trees through (15) to arrive at LAI ζ+ ; 5.
An ensemble of the key LSM parameter, k ζ s , which is updated through (16) over the weekly timescale (∆t 5 ); 6.
Finally, the m ζ a ensembles of grass and trees that are updated through (19) at > weekly (e.g., 3 weeks) timescale (∆t 6 ).
Note that the time step of k ζ s updating was lower than m ζ a updating because the dynamics of the water in the soil layer, especially in the case of thin soil layers, are faster than the slow vegetation change dynamics. Note also that the assimilation procedure of radar backscatter data (step 3) is independent to the assimilation of optical image data (step 4), so that the timing of ε observations and NDVI remote data can be different.
Hereafter, we indicate the ensemble open loop without assimilation (i.e., only steps 1 and 2) as "EnOL". "EnKF" indicates the assimilation approach that includes the ensemble Kalman filter only (i.e., steps 1, 2, 3, and 4), and "EnKFdc" indicates the assimilation approach that includes the six steps described above.
We proved the performance of the assimilation approach for increasing the uncertainty of the model (given the observation errors) so that the proposed assimilation approach will be tested for increasing prescribed errors in k s and the grass and tree m a model parameters (compared with the calibrated values), comparing EnOL, EnKF, and EnKFdc performances.

Application of the Assimilation Approach to the Case Study
In total, Sentinel 1 data were collected for 153 days from January 2016 to August 2018, for which σ 0 vv backscattering coefficients were available. Data were collected, analyzed, and corrected to include the vegetation growth effect on the radar backscatter in Montaldo et al. [33], so that we used the ε time series produced by Montaldo et al. [33].
A total of 124 images of optical sensors were acquired (75 images from Landsat 8 and 49 images from Sentinel 2, see Figure 1c for the 2016-2018 period), from which the NDVI was derived at a 30 m spatial resolution. The coefficients of (1) were estimated using simultaneous NDVI data from remote sensors and LAI observations in the field (a total of 24 simultaneous days) distinguishing grass and trees (β 1 = −0.435, β 2 = 1.014 and β 3 = 0.4029 for grass in the autumn-winter period; β 1 = −0.141, β 2 = 1.720, and β 3 = 1.674 for grass in the spring-summer period; β 1 = 0, β 2 =5.392 and β 3 = 0.486 for trees).
In this case study, the LSM time step was half an hour (∆t 1 ), the VDM time step was one day (∆t 2 ), the assimilation time step of grass and tree NDVI data (∆t 4 ) was variable according to data availability, ranging from 2 days to 20 days with an average of 6 days, and the time step of the m ζ a,g and m ζ a,t updating was 3 weeks for both grass and trees (∆t 6 , Figure 2). Note that, because the VDM was applied distinctly for grass and tree in the Sardinian heterogeneous ecosystem, the grass and tree cells need to be selected in the field as representatives of the two main vegetation components for distinctly assimilating grass and tree NDVI, which were estimated in those cells in the VDM. The assimilation time step of the Sentinel 1 dielectric constant was ≤6 days (∆t 3 ), and k ζ s was updated every 6 days (∆t 5 ). In the EnKF, N e was 100, which is a sufficiently large number for accurate predictions [38,39,79].
We assumed the measurement errors to be with zero mean with a standard deviation of 0.025 for both grass and tree NDVI. The measurement error ε was assumed to be a zero mean with a standard deviation of 0.1, which corresponded to an error of about 5% in the θ observations. In the VDM, we generated the ensembles of the initial LAI values for grass (LAI g ) and tree (LAI t ) from a Gaussian distribution with means of 0.5 and 5.5, respectively, intentionally different from the observations, and a standard deviation σ LAI of 0.2 for both grass and tree LAI. In the LSM, the ensembles of initial θ ζ were generated from a Gaussian distribution with a mean of 0.2; this was intentionally lower (20%) than the observed value, with a standard deviation of 0.05. At each time step, we generated the ensembles of the following: (i) the precipitation, by multiplying the recorded precipitation value by a normally distributed random variable with a zero mean and a standard deviation equal to 20%; (ii) the incoming solar radiation; (iii) the PAR by multiplying the measured values by a normally distributed random variable with mean zero and a standard deviation equal to 10%. It should be noted that the errors of the initial model states and parameters were uncorrelated.

Results
Aerial photography on a dry summer day (July 2016; Figure 1b) successfully depicts the heterogeneity of the field site-a typical Mediterranean landscape, with open areas covered in grass or bare soil depending on the season, and surrounded by trees (in this case, wild olive trees) (Figure 1a). The 30 m spatial resolution of the NDVI map estimated by the Landsat 8 image of a close day (26 July 2016; Figure 1c) was enough to identify the spatial variability of the vegetation cover (differences of 5% of the fraction of tree cover in the footprint). We identified the two representative cells of grass and tree from the NDVI map: the tree cell was the cell with the highest NDVI, while the grass cell, in which the vegetation contribution was absent in the dry summer because bare soil was predominant, was the cell with the lowest NDVI (Figure 1c). During the year, tree NDVI values were always high ranging between~0.6 (in summer) and~0.8 (in spring), while grass NDVI changed widely with the seasons from~0.2 (in summer) to~0.8 (in spring) (Figure 1c). The NDVI data of these representative cells were assimilated in the coupled LSM-VDM model as representative of the two main vegetation components.
We tested the proposed assimilation approach for the case of extremely biased model parameters. In the following, the θ, grass, and tree LAI observations are those estimated from remote sensing data (i.e., using (4) and (1), respectively). We generated initial k ζ s , m ζ a,g , and m ζ a,t ensembles with initialk s ,m a,g , andm a,t values (=5 × 10 −4 m/s, 0.12 d −1 , and 0.01 d −1 , respectively), which were greatly higher than the corresponding calibrated values (Figure 3). The use of the EnKF was not enough for guiding the models, because grass LAI was still underpredicted during the growing seasons, similar to the results with the EnOL configuration; RMSE values were still high, up to 0.8; and observed grass LAI was higher Remote Sens. 2022, 14, 3458 15 of 27 than 1 (Figure 3d,g). The predicted tree LAI using the EnKF was even worse, with values lower than 0.3, while the observed tree LAI from remote data was around 4, and RMSE became close to 4 after 8 months of simulation when the initial model conditions were lost (Figure 3e,h). The results for θ predictions using the EnKF configuration were slightly better when compared with those using the EnOL configuration (Figure 3f,i); however, the model was still underpredicting the soil moisture during the wet months (up to 50% in autumn 2017) when the hydraulic conductivity greatly affected the soil moisture budget predictions due to the prescribed errors in k s . The EnKFdc approach dynamically calibrated the three key model parameters, which converged to values close to the calibrated values after 8-13 months (Figure 3), becoming the coupled model recalibrated for the 2017 and 2018 predictions. In 2016, the model was not yet recalibrated and the RMSE was still high, because the approach requires time for capturing and correcting the persistent model errors. Thanks to the model parameter updating, after almost one year, the RMSE of tree LAI became almost negligible (Figure 3h), and the RMSE of grass LAI decreased to values lower than 0.3 (Figure 3g). Soil moisture was better predicted using the EnKFdc configuration ( Figure 3f), with values close to the radar-observed soil moisture during the wet months due to the correction of the hydraulic conductivity.
Remote Sens. 2022, 13, x FOR PEER REVIEW 17 of 29 moisture budget predictions due to the prescribed errors in ks. The EnKFdc approach dynamically calibrated the three key model parameters, which converged to values close to the calibrated values after 8-13 months (Figure 3), becoming the coupled model recalibrated for the 2017 and 2018 predictions. In 2016, the model was not yet recalibrated and the RMSE was still high, because the approach requires time for capturing and correcting the persistent model errors. Thanks to the model parameter updating, after almost one year, the RMSE of tree LAI became almost negligible (Figure 3h), and the RMSE of grass LAI decreased to values lower than 0.3 (Figure 3g). Soil moisture was better predicted using the EnKFdc configuration (Figure 3f), with values close to the radar-observed soil moisture during the wet months due to the correction of the hydraulic conductivity. Similar supportive results were obtained generating , , , and , ensembles with initial , , , and , values, respectively (= 5 × 10 −8 m/s, 0.0032 d −1 , and 0.0001 d −1 , respectively). These were greatly lower than the corresponding calibrated values (Figure 4). Again, while the EnKF configuration was not enough to guide the model, the use of the EnKFdc approach updated the three model parameters, which reached values close to the corresponding calibrated values after almost one year (Figure 4a-c). Using the EnKFdc approach, the errors of grass and tree LAI predictions when compared to satellite observations were almost negligible (RMSE < 0.15 for grass LAI and RMSE < 0.1 for tree LAI in the 2017 and 2018 years; Figure 4g,h). These also decreased for the soil moisture predictions, especially during the unusually wet 2018 (Figure 4i). , , and , ensembles, respectively, using the EnKFdc approach (the means of the ensembles are in black thick lines; for reference, the calibrated values of ks, ma,g, and ma,t are reported in dotted horizontal lines); (d), (e), and (f) show the comparison between LAI and soil moisture observations derived from assimilated remote data (obs.) and the ensemble means predicted using the EnOL, EnKF, and EnKFdc approaches, respectively; (g), (h), and (i) show the evolutions of the RMSE of the ensemble mean of soil moisture and LAI predicted using the EnOL, EnKF, and EnKFdc approaches concerning the observed soil moisture and LAI (derived from remote sensing data) using a 60-day window, translated in 10-day increments. (d-f) show the comparison between LAI and soil moisture observations derived from assimilated remote data (obs.) and the ensemble means predicted using the EnOL, EnKF, and EnKFdc approaches, respectively; (g-i) show the evolutions of the RMSE of the ensemble mean of soil moisture and LAI predicted using the EnOL, EnKF, and EnKFdc approaches concerning the observed soil moisture and LAI (derived from remote sensing data) using a 60-day window, translated in 10-day increments.
Similar supportive results were obtained generating k ζ s , m ζ a,g , and m ζ a,t ensembles with initialk s ,m a,g , andm a,t values, respectively (=5 × 10 −8 m/s, 0.0032 d −1 , and 0.0001 d −1 , respectively). These were greatly lower than the corresponding calibrated values (Figure 4). Again, while the EnKF configuration was not enough to guide the model, the use of the EnKFdc approach updated the three model parameters, which reached values close to the corresponding calibrated values after almost one year (Figure 4a-c). Using the EnKFdc approach, the errors of grass and tree LAI predictions when compared to satellite observations were almost negligible (RMSE < 0.15 for grass LAI and RMSE < 0.1 for tree LAI in the 2017 and 2018 years; Figure 4g,h). These also decreased for the soil moisture predictions, especially during the unusually wet 2018 (Figure 4i). For fully evaluating the proposed assimilation approach, we compared the performance of the EnKFdc, EnKF, and EnOL approaches on soil moisture and grass and tree LAI predictions for all the ranges of initial , , , , and values, which vary independently ( Figures 5 and 6). For assessing the contribution of radar backscatter and NDVI assimilations separately, we compared the assimilation results in the cases of (i) assimilation of radar backscatter only; (ii) assimilation of grass and tree NDVI only; and (iii) assimilation of both radar backscatter and grass and tree NDVI (Figures 5 and 6). Compared to the EnOL-based results, when radar backscatter was assimilated using the EnKF approach, the soil moisture predictions improved for the whole range of parameters (RMSE < 0.06 for most of the parameter ranges; Figure 5b). These results were improved compared with those assimilating NDVI data only, while still in the EnKF configuration (RMSE was still high, up to 0.12, for high , and , ; Figure 5c). However, the use of the EnKFdc approach further improved the model performance, which was again better when radar backscatter was assimilated. The assimilation of both radar backscatter and grass and tree NDVI enabled us to guide and correct the model for the full range of parameters, with the RMSE of soil moisture always lower than 0.045 (Figure 5g), which is a sort of minimum RMSE value due to the intrinsic errors of the model and observations themselves in soil moisture estimates. For fully evaluating the proposed assimilation approach, we compared the performance of the EnKFdc, EnKF, and EnOL approaches on soil moisture and grass and tree LAI predictions for all the ranges of initialm a,g ,m a,t , andk s values, which vary independently (Figures 5 and 6). For assessing the contribution of radar backscatter and NDVI assimilations separately, we compared the assimilation results in the cases of (i) assimilation of radar backscatter only; (ii) assimilation of grass and tree NDVI only; and (iii) assimilation of both radar backscatter and grass and tree NDVI (Figures 5 and 6). Compared to the EnOL-based results, when radar backscatter was assimilated using the EnKF approach, the soil moisture predictions improved for the whole range of parameters (RMSE < 0.06 for most of the parameter ranges; Figure 5b). These results were improved compared with those assimilating NDVI data only, while still in the EnKF configuration (RMSE was still high, up to 0.12, for highm a,g andm a,t ; Figure 5c). However, the use of the EnKFdc approach further improved the model performance, which was again better when radar backscatter was assimilated. The assimilation of both radar backscatter and grass and tree NDVI enabled us to guide and correct the model for the full range of parameters, with the RMSE of soil moisture always lower than 0.045 (Figure 5g), which is a sort of minimum RMSE value due to the intrinsic errors of the model and observations themselves in soil moisture estimates.
Although the EnKF-based assimilation of NDVI only was already able to sufficiently predict the grass LAI for most parameter ranges with RMSE lower than 1.1 (Figure 6c), the use of the proposed EnKFdc approach-by assimilating both radar backscatter and NDVI-enabled us to make very good predictions of the grass LAI with RMSEs lower than 0.14 ( Figure 6g). Instead, the assimilation of only radar backscatter was not enough for successfully predicting grass LAI using both EnKF and EnKFdc form a,g andm a,t values lower than the calibrated values (Figure 6b,e). Similarly, the use of the proposed EnKFdc assimilation approach predicted the tree LAI well (Figure 6n), with RMSEs lower than 0.1 for the full ranges of parameter values. Instead, again, the use of only the EnKF led to a large misestimate of the tree LAI form a,g andm a,t values higher than the calibrated values ( Figure 6k).
Remote Sens. 2022, 13, x FOR PEER REVIEW 19 of 29 Figure 5. The root mean square error (rmse) of soil moisture (θ) predictions using the (a) EnOL, EnKF, and EnKFdc approaches with the assimilation of the radar backscatter (related to θ) only (b and e), the assimilation of grass and tree NDVI (related to LAI) only (c and f), and the assimilation of both radar backscatter and grass and tree NDVI (d and g), while varying the initial , , , , and values (low ma,g and ma,t values correspond to 0.0032 d −1 and 0.0001 d −1 , respectively, while high ma,g and ma,t values correspond to 0.12 d −1 and 0.01 d −1 , respectively; cal-calibrated values of Table 3).
Although the EnKF-based assimilation of NDVI only was already able to sufficiently predict the grass LAI for most parameter ranges with RMSE lower than 1.1 (Figure 6c), the use of the proposed EnKFdc approach-by assimilating both radar backscatter and NDVI-enabled us to make very good predictions of the grass LAI with RMSEs lower than 0.14 ( Figure 6g). Instead, the assimilation of only radar backscatter was not enough for successfully predicting grass LAI using both EnKF and EnKFdc for , and , values lower than the calibrated values (Figure 6b,e). Similarly, the use of the proposed EnKFdc assimilation approach predicted the tree LAI well (Figure 6n), with RMSEs lower than 0.1 for the full ranges of parameter values. Instead, again, the use of only the EnKF led to a large misestimate of the tree LAI for , and , values higher than the calibrated values (Figure 6k). Figure 5. The root mean square error (rmse) of soil moisture (θ) predictions using the (a) EnOL, EnKF, and EnKFdc approaches with the assimilation of the radar backscatter (related to θ) only (b,e), the assimilation of grass and tree NDVI (related to LAI) only (c,f), and the assimilation of both radar backscatter and grass and tree NDVI (d,g), while varying the initialm a,g ,m a,t , andk s values (low m a,g and m a,t values correspond to 0.0032 d −1 and 0.0001 d −1 , respectively, while high m a,g and m a,t values correspond to 0.12 d −1 and 0.01 d −1 , respectively; cal-calibrated values of Table 3).
The accuracy of the proposed assimilation approach when both radar backscatter and grass and tree NDVI were assimilated was tested seasonally, comparing model predictions of soil moisture and grass and tree LAI using EnOL, EnKF, and EnKFdc approaches in the four seasons for the last two investigated years and using all the combinations of initialk s , m a,g , andm a,t values ( Figure 7). The errors in the soil moisture predictions were removed partially using the EnKF in fall and summer, with a low variability range of RMSE-nearby to the mean RMSE value of~0.06 ( Figure 7a); meanwhile, errors were still high and highly variable in winter and spring (with RMSE values ranging between 0.03 and 0.06). The use of the EnKFdc removed the model bias for the soil moisture predictions in all the seasons, where the RMSE values were always coincidental with the seasonal mean RMSE values (~0.04 in winter and spring and~0.05 in fall and summer; Figure 7a). In a similar way, in all seasons, the bias in the tree LAI predictions was only corrected using the EnKFdc, with RMSE becoming close to 0 and with negligible variability in RMSE using all parameter combinations (Figure 7b). Instead, using the EnKF approach, the bias in the grass LAI predictions was large in the spring-the key season for grass growth-when the EnKFbased approach was not able to guide the model sufficiently well and when the RMSE still reached high values of up to 2 (Figure 8c). Again, only the use of the proposed EnKFdc approach enabled us to completely remove the model bias, which decreased the RMSE values to 0.1 in spring, and showed negligible variability of RMSE values for all parameter combinations (Figure 7c). Note that the statistics of the model performance were computed for the October 2016-September 2018 period-that is, a period after the calibration period of k s , m ag , and m at (e.g., Figure 3). The accuracy of the proposed assimilation approach when both radar backscatter and grass and tree NDVI were assimilated was tested seasonally, comparing model predictions of soil moisture and grass and tree LAI using EnOL, EnKF, and EnKFdc approaches in the four seasons for the last two investigated years and using all the combinations of initial , , , and , values (Figure 7). The errors in the soil moisture predictions were removed partially using the EnKF in fall and summer, with a low variability range of RMSE-nearby to the mean RMSE value of ~0.06 ( Figure 7a); meanwhile, errors were still high and highly variable in winter and spring (with RMSE values ranging between 0.03 and 0.06). The use of the EnKFdc removed the model bias for the soil moisture predictions in all the seasons, where the RMSE values were always coincidental with the seasonal mean RMSE values (~0.04 in winter and spring and ~0.05 in fall and summer; Figure  7a). In a similar way, in all seasons, the bias in the tree LAI predictions was only corrected using the EnKFdc, with RMSE becoming close to 0 and with negligible variability in RMSE using all parameter combinations (Figure 7b). Instead, using the EnKF approach, the bias in the grass LAI predictions was large in the spring-the key season for grass growth- Finally, we evaluated the impact of the use of the proposed assimilation approach on the predictions of two main land surface fluxes-evapotranspiration (ET) and carbon exchange (F c )-which are strictly related to soil moisture and vegetation growth. We used the eddy-covariance-based ET and F c observations to evaluate the model. Although it improved the ET predictions when compared with the predictions using the EnOL configuration, the use of the EnKF was not enough to guide the model and predict the ET well for the full range of parameter combinations: total ET was underpredicted by up to 70% and was overpredicted by up to 10% for high and lowm a,g andm a,t , respectively (Figure 8b). The use of the proposed EnKFdc allowed us to remove the model bias, and the ET was accurately predicted for the full range of the parameters (misprediction of the total ET was lower than 4%, see Figure 8c). ues for all parameter combinations (Figure 7c). Note that the statistics o formance were computed for the October 2016-September 2018 periodafter the calibration period of ks, mag, and mat (e.g., Figure 3). Finally, we evaluated the impact of the use of the proposed assimilat the predictions of two main land surface fluxes-evapotranspiration (ET change (Fc)-which are strictly related to soil moisture and vegetation g the eddy-covariance-based ET and Fc observations to evaluate the model. proved the ET predictions when compared with the predictions using the ration, the use of the EnKF was not enough to guide the model and predic the full range of parameter combinations: total ET was underpredicted b was overpredicted by up to 10% for high and low , and , , respecti The use of the proposed EnKFdc allowed us to remove the model bias, Similar results were obtained for F c predictions. Indeed, the use of the EnKF approach did not remove the model bias in F c predictions when extremely high and lowm a,g andm a,t values were initially assumed: underprediction and overprediction of the total observed F c occurred by up to 80% and 70%, respectively (Figure 8e). The dynamic calibration of the key model parameters using the EnKFdc approach allowed us to remove the model bias and predict the F c well for the full ranges of parameters (misprediction of the total F c was lower than 15%, see Figure 8f). accurately predicted for the full range of the parameters (misprediction of the total ET was lower than 4%, see Figure 8c). Similar results were obtained for Fc predictions. Indeed, the use of the EnKF approach did not remove the model bias in Fc predictions when extremely high and low , and , values were initially assumed: underprediction and overprediction of the total observed Fc occurred by up to 80% and 70%, respectively (Figure 8e). The dynamic calibration of the key model parameters using the EnKFdc approach allowed us to remove the model bias and predict the Fc well for the full ranges of parameters (misprediction of the total Fc was lower than 15%, see Figure 8f).

Discussion
The hydrologic database for the Sardinian site, covering almost 3 years, was very useful for testing the proposed assimilation approach, due to the wide range of hydrometeorological conditions [33], involving an extremely dry 2017 (the soil dried earlier than normal in March, with total precipitation of just 6.9 mm from March to August, see Figure  3) and a wet 2018 (soil dried shortly in July only, with total precipitation of 167.9 mm from March to August, see Figure 3).
The Sardinian field is heterogeneous, with the common Mediterranean tree species wild olives [82][83][84] randomly arranged in space (Figure 1). Although the tree cover size of these tree species was reduced, the freely available Landsat 8 and Sentinel 2 products captured the spatial variability of the NDVI and tree cover in the field (Figure 1). Montaldo et al. [33] demonstrated that the Sentinel 1 satellite can be used for reliable estimates of the soil moisture in the Sardinian site thanks to its fine spatial scale (up to 10 m), as has been shown in other Mediterranean ecosystems [85][86][87]. We demonstrated that the

Discussion
The hydrologic database for the Sardinian site, covering almost 3 years, was very useful for testing the proposed assimilation approach, due to the wide range of hydrometeorological conditions [33], involving an extremely dry 2017 (the soil dried earlier than normal in March, with total precipitation of just 6.9 mm from March to August, see Figure 3) and a wet 2018 (soil dried shortly in July only, with total precipitation of 167.9 mm from March to August, see Figure 3).
The Sardinian field is heterogeneous, with the common Mediterranean tree species wild olives [82][83][84] randomly arranged in space (Figure 1). Although the tree cover size of these tree species was reduced, the freely available Landsat 8 and Sentinel 2 products captured the spatial variability of the NDVI and tree cover in the field (Figure 1). Montaldo et al. [33] demonstrated that the Sentinel 1 satellite can be used for reliable estimates of the soil moisture in the Sardinian site thanks to its fine spatial scale (up to 10 m), as has been shown in other Mediterranean ecosystems [85][86][87]. We demonstrated that the availability of such satellite observations at adequate high spatial and time resolutions allows to monitor the main land surface variables in heterogeneous fields. Furthermore, we have showed that they can be used for guiding ecohydrological modeling, enabling researcher to reach the final objective of an operative data assimilation approach. The proposed data assimilation approach includes the assimilation of both radar backscatter data and NDVI optical data. In this context, Pan et al. [47] and Zhuo et al. [48] assimilated radar Sentinel 1 and optical Sentinel 2 data in the WOFOST model using the EnKF for agricultural land at fine spatial resolution (10-30 m). We demonstrated that both tree and grass NDVI can be assimilated for tree and grass LAI predictions in a heterogeneous field.
The proposed multiscale assimilation approach also includes the dynamic updating through assimilation of three key model parameters-saturated hydraulic conductivity, which is mainly related to soil moisture; and the tree and grass maintenance respiration coefficients, which are mainly related to tree and grass LAI. Montaldo et al. [38], using field measurements as a proxy for remote sensor observations, highlighted the limits of the EnKF model for soil moisture prediction, especially when the key model parameters largely differ from the calibrated values. This is a typical problem of operational data assimilation approaches, which are often used in case studies with limited knowledge of model parameters and field properties [88,89]. The dynamic updating of the model parameters occurs at timescales which are longer than the LSM-VDM timescales, as proposed by Montaldo et al. [38], to capture the persistent bias in soil and biomass budget predictions. Our analysis showed that EnKF-based assimilation performed sufficiently well when the parameters were moderately different from the calibrated values both for soil moisture and for LAI predictions; meanwhile, its performance decreased when the parameters largely diverged from the calibrated values, especially for soil moisture and tree LAI predictions (Figures 5 and 6). EnKF performed sufficiently well for the grass LAI predictions; although, the use of the EnKFdc approach still improved the model's performance for the low initial values of the maintenance respiration coefficients ( Figure 6).
The proposed EnKFdc approach performed well for the full range of parameters and hydrometeorological conditions. The EnKF approach was not enough to guide the models when the key parameters were largely biased, because when, for instance, the saturated hydraulic conductivity is largely over-or underestimated, the root zone soil moisture balance is not systematically preserved, and only the progressive and systematic correction of the key soil parameter-the saturated hydraulic conductivity-can correct and guide the model. Lu et al. [52,53] and Nie et al. [50] demonstrated the efficacy of a multiscale assimilation approach for soil moisture predictions with the dynamic calibration of soil parameters. Here, we included the assimilation of NDVI data through calibrating the main VDM model parameters, i.e., the maintenance respiration coefficients. The dynamic calibration of the m a parameters can guide the VDM, preserving the biomass balance of grass and tree species. The proposed EnKFdc approach for both LSM and VDM was the only approach to predict not only the soil moisture and the tree and grass LAI well, but was also the only approach to predict the main outputs of the coupled model, the evapotranspiration, and the carbon exchanges (Figure 8), which are strictly related to soil moisture and vegetation dynamics [57,90]. The use of the EnKF alone was not enough to obtain good predictions of the two key land surface fluxes, with errors up 70% in ET predictions, when lower m a values were prescribed; furthermore, the errors were up to 120% in F c predictions when high initial m a values were prescribed. Such high model bias in ET predictions affects soil-water balance predictions in these semiarid ecosystems because ET is the main loss term of the soil water budget, with a yearly magnitude that may be equal to the precipitation [6,91,92]. ET and F c predictions were good for the full range of model parameters when the proposed EnKFdc approach was used, showing the high performance of the approach (errors less than 4% and 15% in ET and F c predictions, respectively).
In terms of the importance of the proposed assimilation approach in relation to the seasons, the EnKFdc approach is essential for tree LAI predictions in all the seasons because the tree LAI values of the evergreen tree species are almost constant during the year, and NDVI assimilation is important for the whole year. For grass LAI predictions, the EnKFdc approach is more important in the growing season, when grass achieved maximum height, impacting evapotranspiration and carbon exchanges [14]. The EnKFdc approach becomes even more crucial for soil moisture predictions in spring and winter, when soil is wetter, and the root zone soil budget is controlled by a key soil parameter-saturated hydraulic conductivity [93]. Winter and spring are crucial seasons for water resources management in semiarid Mediterranean ecosystems, when soil moisture and ET need to be adequately predicted due to their key roles in water resources.

Conclusions
This paper proposes a multiscale assimilation approach that assimilates both radar backscatter data and grass and tree NDVI data from remote sensors for guiding the predictions of soil moisture and LAI of coupled vegetation dynamic and land surface models. It is suitable for heterogeneous ecosystems, which are typical of semiarid climates, thanks to the sufficiently high spatial (10-30 m) and temporal (~weekly) resolutions of the observations of the Sentinel 1 radar, the Landsat 8 satellite, and the Sentinel 2 satellite. The assimilation approach which is based on the ensemble Kalman filter (EnKF) was not sufficient for guiding soil moisture and LAI predictions when the model was largely biased due to the values of three key model parameters-the saturated hydraulic conductivity and the grass and tree maintenance respiration coefficients-being largely different from the calibrated values. The proposed multiscale assimilation approach is not limited to assimilating remote sensing data for model predictions using the EnKF: it uses assimilated data for dynamically updating key model parameters (the ENKFdc approach) at a longer timescale. These key model parameters are the saturated hydraulic conductivity and the grass and tree maintenance respiration coefficients, which are highly sensitive parameters of soil-water balance and biomass budget models, respectively. The use of the proposed EnKFdc approach was essential for soil moisture and grass LAI predictions, especially in the winter and spring seasons, which are key seasons for the water resources management of Mediterranean semiarid ecosystems. We demonstrated that the use of the proposed assimilation approach enabled us to predict key land surface fluxes, drastically reducing model prediction errors when parameters are wrongly estimated. From these results, we can also anticipate that the proposed assimilation approach may be even more important in ecosystems under wetter climates, where the model bias effects can be even larger on a soil water budget due to misestimates of hydraulic conductivity.