Improving soil moisture profile prediction from ground-penetrating radar data : a maximum likelihood ensemble filter approach

Introduction Conclusions References


Introduction
Understanding the dynamics of soil moisture at the root zone is essential for hydrological, meteorological and agricultural researches.The water content at this zone influences on most important processes of the hydrological cycle such as infiltration, runoff and evaporation as well as partitioning of energy at the land surface into sensible and latent exchange with the atmosphere (Vereecken et al., 2008;Lambot et al., 2009).The availability of the root zone water is the main factor that controls the separation of the rainfall into runoff and infiltration.Disregarding the spatial pattern of antecedent soil moisture may cause significant errors on runoff prediction (e.g.Merz and Bardossy, 1998;Minet et al., 2010).In agriculture, information on the root zone soil moisture is crucial for an optimal management and irrigation practices for its substantial impact on the crop production (Vereecken et al., 2008).As a result, development and integration of measurement techniques for quantitative characterization of the root zone soil moisture is an urgent need.
During last two decades, GPR has become a popular noninvasive technique that is widely applied in many engineering applications (e.g.bridge and road evaluation, buried pipe location, crack inspection and landmine detection) (Slater and Comas, 2009).In soil and water sciences, GPR has been used to provide reliable, highly spatially-resolved soil moisture at field scale (Huisman et al., 2003;Robinson et al., 2008).The theoretical foundation underpinning the applications of GPR for soil moisture estimation is the overwhelmingness of the water permittivity to the other soil components, which makes soil water content mainly governs the electromagnetic wave propagation in the soil.Several methods were developed to interpret the GPR signals into soil water content.Detailed reviews on these methods can be found in Huisman et al. (2003); Slater and Comas (2009).Generally, these methods can be classified into three approaches.Firstly, the water content is derived from the propagation velocity of the ground wave, which travels from transmitting to receiving antennas through the soil surface (Grote et al., 2003;Galagedara et al., 2005;Mangel et al., 2012;Pan et al., Figures Back Close Full 2012; Steelman and Endres, 2012).The second approach estimates the soil moisture from the surface reflection coefficient, which is calculated as the ratio between the reflection amplitude from the soil surface and that from the calibrating perfect electric conductor (PEC) (Serbin andOr, 2004, 2005).However, due to the significant assumptions related to electromagnetic wave propagation and only using a part of information contained in the GPR data, the accuracy of these approaches are very limited.Recently, Lambot et al. (2004b) developed a more advanced method known as full-wave inversion of GPR data, which permits to use all GPR information for the soil moisture estimation.The method is based on an analytical antenna model, which represents an exact solution of Maxwell's equations in far-field conditions.The model allows to reproduce accurately the GPR signals reflected from a planar multi-layered medium.Very good results were obtained for the case of a 1-layered soil model (e.g.Minet et al., 2011;Tran et al., 2012).However, for 2-layered or soil moisture profile models, the soil moisture estimation at the lower layers was not adequately accurate due to the ill-posedness of the inverse problem caused by too many unknown parameters and low sensitivity of GPR data with the electrical properties at the deeper layers (Lambot et al., 2004a;Minet et al., 2011).Efforts have been made to increase the accuracy of the soil moisture profile estimation by using the time-lapse GPR data to constrain a soil hydrodynamic model (Kowalsky et al., 2005;Lambot et al., 2006;Jadoon et al., 2008;Looms et al., 2008;Lambot et al., 2009).A joint inversion procedure based on the integrated geophysical and hydrodynamic models was demonstrated to optimally estimate the soil unsaturated hydraulic parameters.However, in case the soil moisture profile slightly varies with the soil depth, some parameters were suffering from relatively large uncertainties due to their low sensitivity to the radar data (Lambot et al., 2006).
With the development of the data assimilation algorithms, there have been increasingly intensive researches on assimilation of remote-sensed data to improve soil moisture profile prediction (Loew, 2008).For instance, Hoeben and Troch (2000); Walker et al. (2001) applied the Extended Kalman filter (EKF) to retrieve the soil moisture profile from the near-surface soil moisture measurements.Reichle et al. (2002); Das and Introduction

Conclusions References
Tables Figures

Back Close
Full  2012) assimilated the soil moisture data from remote sensing measurements with the ensemble Kalman filter (EnKF).Application of the variational assimilation methods for updating the soil moisture profile was performed by Reichle et al. (2001); Sabater et al. (2007).Montzka et al. (2011) employed the particle filter technique to simultaneously update the soil hydraulic properties and soil moisture profile with the measurements from geophysical and satellite measurements.Most of these studies used only the surface soil moisture, which is obtained by processing the remotely-sensed data, to update the whole soil moisture profile.This reduces the improvement of the assimilation due to lack of information for the assimilation (Pauwels et al., 2007;Das et al., 2008).In addition, errors caused by the conversion from the remote-sensed data to surface soil moisture also contribute to the smaller effectiveness of the assimilation.There are only few researches that directly use the geophysical data for the assimilation (e.g.Reichle et al., 2002;Huang et al., 2008;Rings et al., 2010).Reviews on the assimilation of geophysical techniques to soil moisture profile are given in Loew (2008); Reichle (2008); Vereecken et al. (2008).
In this paper, we propose a new sequential assimilation procedure to raise the accuracy of the soil moisture profile prediction using time-lapse GPR data.The hydrodynamic model, Hydrus-1D ( Šimunek et al., 2009), was employed to simulate the vertical dynamics of the soil moisture in the root zone.The GPR model works as an observation operator to relate the soil moisture profile with the GPR observation.The assimilation was performed within the maximum likelihood ensemble filter (MLEF) framework developed by Zupanski (2005), for which the problem of nonlinear observation operator is solved much more effectively than the EnKF techniques.The method estimates the optimal state as the maximum of the probability density function (PDF) instead of the minimum variance like in most of the other ensemble data assimilation methods (Zupanski et al., 2008).Direct assimilation of GPR data is a prominent advantage of our approach.It avoids solving the time-consuming inverse problem as well as the estimation errors of the soil moisture caused by inversion.In addition, instead of using Introduction

Conclusions References
Tables Figures

Back Close
Full only surface soil moisture, the approach allows to use the information of the whole soil moisture profile, which is reflected via the ultra-wideband (UWB) GPR data, for the assimilation.The use of the UWB antenna in this study is also an advantage as it provides more information about soil moisture profile with a better depth resolution compared to other classical remote sensing techniques (Lambot et al., 2004a).Consequently, this approach is expected to provide better estimation than the use of the surface soil moisture only.The assimilation procedure was numerically tested by synthetic simulations to solve the "wrong" initial conditions.The relationship between the assimilation and update interval was also investigated.To the best of our knowledge, this is the first attempt to directly assimilate GPR measurements for updating the soil moisture profile.

Materials and methods
Figure 1 shows the assimilation procedure that we used to update the state of the soil moisture profile using GPR data.The procedure consists of 3 main components: the hydrodynamic model, the observation operator (electromagnetic model) and the MLEF data assimilation algorithm.

The hydrodynamic model
In this study, the one-dimensional vertical flow in an homogeneous soil column was simulated by the water flow module of Hydrus-1D model ( Šimunek et al., 2009).The model numerically solves Richards's equation, which governs the movement of water in the unsaturated zone: where h is the water pressure head, θ is the volumetric water content, t is time, z is the spatial coordinate taken positive upward, and K (θ) is the unsaturated hydraulic Introduction

Conclusions References
Tables Figures

Back Close
Full conductivity function.C(h) = ∂θ(h)/∂h is the differential water capacity with θ(h) being the water retention curve.The unsaturated hydraulic conductivity and water retention curves are described by the Mualem-van Genuchten model (Mualem, 1976;van Genuchten, 1980): where θ r and θ s are, respectively, the residual and saturated water contents, K s is the saturated hydraulic conductivity, α is the inverse of the air-entry value, n is a pore-size distribution index and l is the pore connectivity parameter.In Hydrus model, the numerical solution of Eq. ( 1) is obtained by using the standard Galerkin-type linear finite element schemes.The atmospheric boundary with surface layer and free drainage were, respectively, selected for the upper and lower boundary conditions.For the potential evaporation, we used the Penman-Monteith equation (Monteith, 1981).

The observation operator
In this study, we updated the prediction of the soil moisture profile by directly assimilating the GPR data.The radar electromagnetic model (Lambot et al., 2004b) worked as an observation operator to relate the soil moisture profile with the GPR data.We Introduction

Conclusions References
Tables Figures

Back Close
Full adopted the GPR measurement system as described in Lambot et al. (2004b).We simulated radar data for a Vivaldi antenna in the frequency range 1-3 GHz with a frequency step of 6 MHz.The distance between the antenna and medium was 37 cm.GPR data are represented by the Green's function G ↑ xx (f , b), which is defined as the backscattered, x-directed electric field at the antenna phase center for a unit-strength, x-directed electric source situated at the same position above the multilayered medium.It is an analytical solution of Maxwell's equations for wave propagation in 3-D multilayered media, which is derived by a recursive scheme to compute the transverse electric and magnetic global reflection coefficients of the multilayered medium in the spectral domain given the parameter vector b (Lambot et al., 2004b(Lambot et al., , 2007)): where subscript 0 refers to the upper half-space (free-space), z 0 is the distance between the antenna phase center and the first medium interface, R TM and R TE are, respectively, the transverse magnetic (TM) and transverse electric (TE) global reflection coefficients accounting for all reflections in the multilayered medium, k ρ is the radial component of the polar coordinate system in the spectral domain.Γ is the vertical wavenumber defined as with the magnetic permeability µ, dielectric permittivity ε and electrical conductivity σ.For the free-space layer 0 (upper half-space), we have k 2 0 = ( ω c ) 2 with c being the free-space wave velocity.We refer to Lambot et al. (2004b) for more detailed information about the Green's function evaluation.
The vector b = [ε n , σ n , h n ], n = 1, . . ., N contains the layer thickness (h n ), and constitutive parameters governing wave propagation in the medium, namely, dielectric permittivity (ε n ) and electrical conductivity (σ n ).The relationship between the soil volumetric water content (θ n ) and the permittivity was formulated by the model of Ledieu et al. (1986) as: Introduction

Conclusions References
Tables Figures

Back Close
Full where we fixed a = 0.1181 and b = −0.1841,which are suitable parameters for a wide variety of soils (Ferr é et al., 1996).The relationship between the soil moisture and electrical conductivity was formulated by the model of Rhoades et al. (1976): where σ w = 0.075 S m −1 is the electrical conductivity of the soil water (Jadoon et al., 2008;Minet et al., 2011), and σ s is the electrical conductivity of the dry soil.Table 1 presents the σ s and two empirical parameters c and d for the loamy sand, silt and clay soils.Due to the high electrical conductivity of the dry soil (σ s = 4.39 × 10 −2 S m −1 ), the electrical conductivity of the clay soil is relatively large, which increases the losses in GPR wave propagation.Consequently, for a given soil moisture profile, we can obtain the dielectric permittivity and electrical conductivity profiles by using Eqs.( 6) and ( 7), and afterwards, the GPR data reflected from that profile as well from Eq. ( 5).
For the fact that the assimilation procedure works with real number, we used the absolute values of the complex GPR data in the frequency domain.It is clear that the observation operator relating the Green's function and soil moisture is a complex nonlinear operator.It is also worth noting that in realistic measurements, the antenna effects are filtered out from the radar data to obtain the Green's function using Eq. ( 1) in Lambot et al. (2004b).

The maximum likelihood ensemble filter
MLEF is an alternative deterministic ensemble filter technique based on control theory (Zupanski et al., 2008).It is based on a combination of the maximum likelihood and ensemble data assimilation.The MLEF is a posterior maximum likelihood approach, Introduction

Conclusions References
Tables Figures

Back Close
Full in which the optimal analysis state is obtained as the maximum of the PDF, and determined by minimizing the cost function derived from a multivariate posteriori PDF.
The approach share the idea to use the minimization of the cost function to derive the analysis model state with the variational data assimilation.However, as an ensemblebased approach, the minimization is performed in an ensemble-spanned subspace, instead of the full model space like in the variational approach.Since the MLEF uses the maximum likelihood estimation to obtain the analysis model state, it allows to effectively solve the nonlinearity of both model and observation operator.The following paragraphs present the basic steps of the MLEF algorithm.Detailed explanation can be found in Zupanski (2005); Zupanski et al. (2008).
Generally, similar to the other assimilation algorithms, the MLEF includes two steps, the forecast and the analysis: Forecast: the forecast step uses the model M (Hydrus-1D) to propagate the system state in time: where x is the N S × 1 state variable vector with N S being the number of state variables, t and t + 1 represent the current and next time steps, respectively.Superscript f and a stand for the forecast and analysis.ω t is the model error vector, which was neglected in our analysis.The i th column of the square-root N S × N E forecast covariance matrix P 1/2 f at time t + 1 is calculated from the ensemble forecast as below: where p a i ,t is the i th column of the square-root, N S ×N E analysis error covariance matrix P Analysis: as soon as the observation is available, the analysis step is performed to update the forecast state variable by the maximum likelihood approach.Accordingly, the MLEF seeks the optimal state variable that maximizes the posterior probability distribution, or in the other words, minimizing the cost function given the Gaussian PDFs for the observation and forecast errors: where the increment vector x − x f can be expressed as a linear combination of the forecast ensemble perturbations P 1/2 f : where w = (w 1 , w 2 , • • •, w N E ) T is the weighting coefficient vector, y denotes the measurement vector (GPR data) and H(x) the observation operator (radar electromagnetic model).R represents the N O × N O measurement covariance matrix with N O being the length of measurement dataset.Superscript T represents the transpose operator.
Our objective is to find the optimal weighting coefficient vector w that minimizes the cost function (Eq.10).In this study, we followed the non-differentiable minimization approach in Zupanski et al. (2008), which estimates the unknown variables by the generalized nonlinear conjugate-gradient optimization algorithm employing the generalized first derivative, calculated as: in which the i th column of the N O × N E observation perturbation matrix Z, z i is calculated as:

Conclusions References
Tables Figures

Back Close
Full Since the optimal preconditioning is defined as an inverse square-root Hessian matrix, a changing variable is introduced using the inverse square-root Hessian matrix as a multiple factor: with the Hessian matrix being calculated ∇ 2 G J(x f ) as: where I is an N E × N E identity matrix and Z(x f ) is calculated by substituting x = x f into Eq.( 13).Now, the optimization of w is implemented via the changing variable ξ.If the observation operator is linear and the changing variable (Eq.14) is employed, the solution of the optimization problem (Eq.10) is obtained in a single step of minimization iteration (Zupanski, 2005).
Once we obtain the analysis state variable at time t + 1, the square-root analysis covariance matrix (P 1/2 a,t+1 ) is updated by: where x opt a represents the optimal analysis state variable.The columns of P 1/2 a are then used in Eqs.(8, 9) for the next analysis cycle.

Numerical simulation
The proposed assimilation procedure was validated by performing numerical simulations, for which the true system was exactly known.We considered a synthetic homogeneous soil column with a depth of 80 cm.This soil column was discretized into 32 equidistant elements.We analyzed the effects of the soil type on the data assimilation Introduction

Conclusions References
Tables Figures

Back Close
Full by considering 3 soil types, namely, loamy sand, silt and clay.Table 1 presents the 6 parameters of the Mualem-Van Genuchten's equation to construct the water retention curve and calculate the unsaturated hydraulic conductivity.These parameters were obtained from Schaap et al. (2001).The simulation period was set at 720 h (30 days).For realistic simulations, the hourly rainfall, temperature, humidity and wind speed data were taken from a meteorological station in Louvain-la-Neuve (Belgium) from 1 April to 30 April 2011.Figure 2 presents the rainfall and potential evaporation data (ET 0 ) in this period.The figure shows that there are 3 main rainfall events during the simulation period in which 2 large events occurred at the beginning and end of the period with the maximum rainfall being 0.46 cm.
The potential evaporation fluctuates proportionally with the variation of the temperature.It is relatively small at time from 0 to 400 h but becomes larger from 400 to 720 h due to the increase of the temperature.The minimum and maximum temperature are 2.5 and 26.5 • C, respectively.The average temperature over the period is around 13.8 • C.
In this study, we used the assimilation procedure to solve the problem of the unknown initial conditions.We assumed that the "true" initial profile was constant with θ t=0 (z) = 0.2 cm 3 cm −3 , while that of the "forecast" was 0.3 cm 3 cm −3 .From the "true" initial state, the 720 hourly "true" soil moisture profiles were generated using the Hydrus-1D model.Figure 3a, c, e presents the variation of the "true" soil moisture at several depths versus time over the simulation period for the loamy sand, silt and clay soils.In general, there is a good correlation between the forcing data (rainfall and evaporation) and moisture of all types of soil (with time lag).The surface soil moisture (z = 0 cm) is the most sensitive to the variation of the rainfall.The peak of the surface soil moisture rapidly appears right after the peak of the rainfall.During dry period, it sharply reduces due to the evaporation and infiltration processes.When the soil depth increases, the time lag between the rainfall and soil moisture peaks gradually increases.This is explained by the fact that water requires more time to move from the surface to lower soil layers.The soil moisture at the deeper layers also varies more smoothly with time due to the Introduction

Conclusions References
Tables Figures

Back Close
Full water retention effect of the soil.The effects of the evaporation on the soil moisture is only observed at the near surface depths (z = 0, 2.5, 5 cm) during the dry period.The temporal variation of the soil moisture of the three soil types are different.Compared to the silt and clay soils, the loamy sand dries out more quickly, and the differences of the water content among the soil depths are also smaller.This is because this soil type has a low water-holding capacity.By contrast, the clay soil shows its ability to well retain water.As the rainfall even occurs, the water content at the near surface layers (from 0 to 10 cm) increases rapidly but it slowly moves down to the lower layers.As a result, the water content at depths 40 and 80 cm is constant at the initial state during the simulation period.
For each "true" soil moisture profile, the frequency-domain synthetic GPR data were generated using the forward radar model.In this model, the soil moisture profile was modeled by 32-planar layers with an equal thickness of 2.5 cm under the air layer and above a lower half space.
Figure 3b, d, f presents the "observed" frequency-domain absolute values of GPR data (Green's function) corresponding to the three soil types (loamy sand, silt and clay).For each soil type, there are 720 GPR datasets generated from 720 hourly soil moisture profiles.Each dataset has 334 "observed" values of absolute GPR data in a frequency range of 1-3 GHz and a step of 6 MHz.It is worth noting that in realistic practice, in addition to the reflections from the medium (Green's function), the measured GPR data also include the GPR effects.Therefore, the "observed" GPR data in this study correspond to the measured GPR data after filtering the antenna effects.Generally, the figure shows that the absolute value of the Green's function is higher at the higher frequencies, implying that the reflections from the surface layer are stronger than those from the deeper layers.Comparing Fig. 3b, d, e and 3c, d, f, we clearly see that the temporal variation of the GPR data agrees well with that of the soil moisture.
The absolute Green's function increases when the soil is wet and vice versa.The 3 peaks of the GPR data occur at approximately the same time as those of the surface soil moisture.This is attributed to the fact that the amplitude of the GPR reflections Introduction

Conclusions References
Tables Figures

Back Close
Full is stronger as there is an higher contrast in the dielectric permittivity (and the corresponding soil moisture) between the soil and the air layers.The good correspondence between the GPR and soil moisture data indicates a great potential of GPR to correct the prediction of the soil moisture profile.The figure also indicates that at the high frequency part (2-3 GHz), the absolute Green's function of the clay is larger than that of the silt and and loamy sand.This can be explained by the higher water content at the upper part of the clay soil moisture profile.
As for the measurement error covariance, we assumed that only the elements in the main diagonal (i.e.variances) of the observation error covariance matrix R are different from zero.These elements were calculated corresponding to a constant variance of each element in the the soil moisture profile, σ 2 = 0.013 2 .
We began to assimilate the GPR data to update the state of the soil moisture profile at a time step of 80 h. Figure 4a, b, respectively, presents the synthetic and forecast soil moisture profiles at initial time (0 h) and at time that the GPR data assimilation begins (80 h) for the loamy sand, silt and clay soils.Both figures show that the soil moisture profiles corresponding to the three soil types are different due to their different hydraulic parameters, though the initial conditions and forcing data are identical.Under impact of the same amount of rainfall, the discrepancy between the upper and lower parts of the clay soil moisture profile at time step 80 h is much larger than those of the loamy sand and silt soils.
The figure also shows that, compared to the silt and clay soils, after 80 h simulation, the gaps between the "true" and forecast soil moisture profiles of the loamy sand is the smallest.This indicates that the predicted soil moisture profile of the loamy sand soil can converge to the "true" profile without data assimilation if the simulation time is long enough.
In this study, we used an ensemble including 10 members to estimate the analysis soil moisture.The initial ensemble was generated from the initial state with the perturbations being assumed as white-noise Gaussian random variables.The standard deviation of the initial perturbations was set to 20 % of the initial state.Introduction

Conclusions References
Tables Figures

Back Close
Full To illustrate the effect of the update interval on the assimilation, we performed the data assimilation every 5, 30, and 50 h.

Assimilation results
Figure 5 compares the "true" and forecast soil moisture profiles of the loamy sand, silt and clay at several time steps with an update frequency of 5 h.For comparison purpose, in addition to the GPR data assimilation, the results obtained by the open-loop prediction were also presented.The open-loop is the prediction in which the state system is propagated using the initial conditions and forcing data without data assimilation (Walker et al., 2001).
As for the effectiveness of the assimilation, the figure shows that the performance of the GPR assimilation is much better than the open-loop prediction for all soil types.However, this performance is different for loamy sand, silt and clay soils.At time step 80 h (1st assimilation cycle), the assimilation significantly improves the forecast soil moisture profile.The updated soil moisture profiles of the 3 soil types approach the "true" state much more closely than the open-loop profiles.At time step 180 h (21th assimilation cycle), while the discrepancies between the updated and "true" soil moisture profile are clearly observed for the clay soil, they are invisible for the loamy sand and relatively small for the silt soil.With respect to the open-loop prediction, there are still large gaps between the open-loop forecast and "true" soil moisture profile of the silt and clay soils.Contrastingly, the open-loop moisture profile of the loamy sand rapidly approaches the "true" state, indicating the wrong initial conditions can be effaced for this soil if the simulation time is long enough.Until a time step of 280 h (41th assimilation cycle), for the loamy sand and silt soils, the updated and "true" soil moisture profiles approximately identical.For the clay soil, a good agreement between the updated and "true" soil moisture profile is found from the surface to a depth of 20 cm.At the deeper Introduction

Conclusions References
Tables Figures

Back Close
Full layers, the errors are around 0.008 cm 3 cm −3 .At this time step, the differences between the open-loop prediction and "true" state are negligible for the loamy sand, while these differences are still large for the silt and clay soils.At the end of the simulation period (time step = 720 h), it is impossible to separate the "true" soil moisture profile of the loamy sand from the updated and open-loop ones.A perfect agreement between the updated and "true" soil moisture is also observed for the silt soil.However, the errors of the open-loop prediction for this soil type are relatively large (around 0.04 cm 3 cm −3 ), though they are smaller than at time step 280 h.As for the clay soil, no much improvement is observed for both of the assimilation and open-loop prediction, compared to those observed at time step 280 h.
The figure also presents that the upper part of the soil moisture profile is better corrected than its lower counterpart, which is clearly seen for the clay soil.This can be explained by the fact that the GPR data for the assimilation mostly reflect the variations of the soil moisture at near soil surface due to the electrical and dielectric losses.The correction of the deeper soil moisture profile is mainly based on the hydrodynamic interactions between the upper and lower layers.Consequently, for the soils with a good downward drainage (loamy sand and silt in this study), the variations of the soil moisture at the upper layers quickly influence on the lower layers.This leads to an easy correction of the wrong initial problem.By contrast, the variation of the soil moisture at the upper layers takes much longer time to propagate downward for the clay soil.Additionally, compared to the loamy sand and silt soils, the clay soil has the higher electrical conductivity, which causes the larger electrical losses and, therefore, reduces the GPR penetrating depth.As a result, it is more difficult to correct the whole soil moisture profile.
The comparison between the "true" and the forecast with and without GPR assimilation was quantitatively evaluated by the root mean square error (RMSE), which is formulated as: Introduction

Conclusions References
Tables Figures

Back Close
Full Figure 6 presents the RMSE of the soil moisture profile obtained by the open-loop forecast and GPR data assimilation as a function of the simulation time for the 3 soil types.The update interval is 5 h.Comparing Fig. 6a and 6b, we observed that, for all soil types, the RMSE of the open-loop prediction is much higher than that of the data assimilation.The smallest difference between the open-loop and GPR assimilation was found for the loamy sand due to its ability to self-converge to the "true" state.The RMSE of the open-loop prediction for this soil type ranges from 3 to 8 times larger than that of the GPR data assimilation.This quantity is 9-13 times for the clay, and 14-93 times for the silt soil.
As for the relationship between the RMSE and soil type, the figure shows that the temporal variation of the RMSE both for the open-loop and GPR assimilation relies on the hydraulic parameters of the soil.The rapid change of the water pressure head with the soil moisture helps the RMSE to quickly reduces with the simulation time, i.e. the wrong initial conditions is easy to be corrected.Indeed, for the loamy sand, the RMSE both for the open-loop and GPR assimilation sharply reduces with time and converges to the "true" soil moisture profile (with RMSE < 0.002) at time step 575 h (for the openloop prediction) and 165 h (for the GPR assimilation).For the silt soil, while the RMSE for the open-loop prediction steadily decreases from 0.087 to 0.039, that for the GPR data assimilation quickly reduces and approximately equal to the "true" state at time step 250 h.The RMSE for the clay soil remains stable over the simulation period for the case of the open-loop prediction, indicating that the wrong initial conditions cannot be self-corrected in this simulation period.After the first assimilation cycle at time step 80 h, the RMSE reduces suddenly from around 0.1 to 0.01 and continue to decrease to 0.0076 at time step 150 h.From time step 150 h to the end of the simulation period, the Introduction

Conclusions References
Tables Figures

Back Close
Full of the soil moisture profile has a low effect on its lower part due to the slow downward movement of the wetting front.

Effect of the update interval
The effects of the update interval on the GPR data assimilation is shown in Fig. 7, which plots the temporal variation of the RMSE for the GPR assimilation corresponding to update frequencies of 5, 30 and 50 h.It is worth noting that, for each soil type, the RMSE of all update frequencies at the first assimilation cycle are identical because they started assimilating at the same time step (80 h).The figure shows that the RMSE of the updated soil moisture profile for all update frequencies is very small, indicating that the wrong initial condition can be effectively corrected by assimilating GPR data even with the update interval up to more than 2 days.
As for the relationship between the assimilation and update interval, for the three soil types, the figure shows a trend that the RMSE increases as the update interval reduces.However, this dependence is different for different soil types.For the loamy sand, the obtained results indicate that for this "wrong" initial condition problem, increasing the update interval does not much improve the effectiveness of the GPR data assimilation.The difference of the RMSE of the soil moisture profile obtained from the update frequencies of 5 and 50 h is lower than 5 × 10 compared to that observed for the loamy sand and silt soil.As a result, increasing update interval enables us to more closely follow the soil moisture dynamics.In addition, as mentioned, for its high drainage ability, the discrepancy between "true" and forecast states of the loamy sand considerably reduces with the simulation time even without data assimilation, whilst this discrepancy slowly reduces for the silt and levels off for the clay soil.Consequently, while increasing update interval insignificantly impacts on the effectiveness of the GPR data assimilation for the loamy sand, it corrects better the predicted soil moisture profile for the silt and clay soils.

Conclusion
In this paper, we explored the potentials of assimilating UWB GPR data based on the MLEF technique.Resorting to an accurate radar model for wave propagation in planar layered media, the approach directly assimilates the GPR data to update the soil moisture profile.This is different from the common approaches which assimilate the surface soil moisture obtained from remote sensing techniques to update the whole soil moisture profile.Our approach allows us to avoid solving the difficult, time-consuming inverse problem to estimate the soil moisture from GPR data.We validated our approach using synthetic experiments for the 3 typical soil types (loamy sand, silt and clay) in which the initial conditions were assumed to be wrong, as they are usually not known in practice.The obtained results demonstrated that the RMSE of the prediction without the GPR data assimilation is 3-8 times larger than that with the GPR data assimilation for the loamy sand, 9-13 times for the clay, and 14-93 times for the silt soil.
The results also reveal that the effectiveness of the GPR data assimilation depends on the hydraulic properties of the soil type.Due to its high drainage capacity, the updated soil moisture profile of the loamy sand quickly converges to "true" state.In contrast, the updated moisture profile of the clay soil takes much longer time to approach the "true" profile, especially at the deep layers.Introduction

Conclusions References
Tables Figures

Back Close
Full With respect to the relationship between the assimilation and update interval, the obtained results show that for the wrong initial condition problem, increasing update interval slightly improves the updated soil moisture profile for the loamy sand with the difference of the RMSE between updated frequencies of 5 and 50 h lower than 5×10 −4 , while it significantly increases the effectiveness of the GPR data assimilation for the clay soil with the maximum difference of the RMSE among the update frequencies of 0.0012.The success of the proposed approach appears to be promising for using GPR data to improve the real-time soil moisture profile prediction at the field scale, which is very important in water and agricultural management.However, this study only solved the problem of wrong initial condition.The soil hydraulic parameters were assumed to be exactly known, which does not usually occur in the practice.In addition, the model errors were ignored.Our next research will concentrate on the assimilation of GPR data to simultaneously update both soil hydraulic parameters and state variables for different soil types under the model errors.Full  Full Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Mohanty (2006); De Lannoy et al. (2007); Huang et al. (2008); Crow et al. (2008); Draper et al. ( Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | at time t.N E denotes the number of members of the ensemble.It is worth noting that at the initial time step, the analysis model state (x a t=0 ) and its associated uncertainty (p a i ,t=0 ) are needed to a priori determined for the fact that these values are not available.Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | −4 .As time goes by, these differences gradually reduce and at the time step of around 440 h, the RMSE of all update frequencies is approximately equal.When it comes to the silt soil, the impact of the update interval becomes clearer.The maximum difference of the RMSE between the update frequencies of 5 and 50 h increases to 0.0012 (at time step 130 h) and until the end of the simulation period, the difference is still visible.The strongest effects of the update interval on the GPR data assimilation is found for the clay soil.The difference of the RMSE between update frequencies of 5 and 50 h reduces slowly from 0.002 at time step 130 h down to 0.0018 at time step 680 h.These facts can be explained by observing Fig. 3.The figure shows that the temporal variation of the surface layers of the clay soil, which is the most sensitive layers with the GPR signal, is very large 1599 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

Fig. 1 .
Fig. 1.Assimilation procedure using GPR data to update the state of the soil moisture profile.

Fig. 1 .Fig. 2 .
Fig. 1.Assimilation procedure using GPR data to update the state of the soil moisture profile.

Fig. 4 .Fig. 6 .
Fig. 4. Synthetic (left) and forecast (right) soil moisture profiles at initial time (0 h) and at time that the GPR data assimilation begins (80 h) corresponding with the loamy sand, silt and clay soil types.

Fig. 6 .Fig. 7 .
Fig. 6.The RMSE of the open-loop forecast (left) and updated (right) soil moisture profiles for clay (blue), loamy sand (red) and silt (black) soils as a function of the simulation time.

Fig. 7 .
Fig. 7.The RMSE of the soil moisture profile obtained by the assimilation of GPR data with different update frequencies, namely 5, 30 and 50 h.The three soil types were compared, namely, (a) loamy sand, (b) silt and (c) clay.

Table 1 .
VanGenuchten's hydraulic parameters and Rhoades's petrophysical parameters for the loamy sand, silt and clay soils.