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

The vertical profile of shallow unsaturated zone soil moisture plays a key role in many hydro-meteorological and agricultural applications. We propose a closed-loop data assimilation procedure based on the maximum likelihood ensemble filter algorithm to update the vertical soil moisture profile from time-lapse ground-penetrating radar (GPR) data. A hydrodynamic model is used to propagate the system state in time and a radar electromagnetic model and petrophysical relationships to link the state variable with the observation data, which enables us to directly assimilate the GPR data. Instead of using the surface soil moisture only, the approach allows to use the information of the whole soil moisture profile for the assimilation. We validated our approach through a synthetic study. We constructed a synthetic soil column with a depth of 80 cm and analyzed the effects of the soil type on the data assimilation by considering 3 soil types, namely, loamy sand, silt and clay. The assimilation of GPR data was performed to solve the problem of unknown initial conditions. The numerical soil moisture profiles generated by the Hydrus-1D model were used by the GPR model to produce the “observed” GPR data. The results show that the soil moisture profile obtained by assimilating the GPR data is much better than that of an open-loop forecast. Compared to the loamy sand and silt, the updated soil moisture profile of the clay soil converges to the true state much more slowly. Decreasing the update interval from 60 down to 10 h only slightly improves the effectiveness of the GPR data assimilation for the loamy sand but significantly for the clay soil. The proposed approach appears to be promising to improve real-time prediction of the soil moisture profiles as well as to provide effective estimates of the unsaturated hydraulic properties at the field scale from time-lapse GPR measurements.


Introduction
Understanding the dynamics of soil moisture at the shallow unsaturated zone is essential for hydrological, meteorological and agricultural research.The water content at this zone influences the most important processes of the hydrological cycle as well as partitioning of energy at the land surface into a sensible and latent exchange with the atmosphere (Vereecken et al., 2008;Lambot et al., 2009).The availability of the unsaturated 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., 2011a).In agriculture, information on the unsaturated zone soil moisture is crucial for optimal management and irrigation practices toward a tangible impact on the crop production (Vereecken et al., 2008).As a result, development and integration of measurement techniques for quantitative characterization of the shallow unsaturated 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 land mine 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 the main characteristic that 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) and 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., 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.Lambot et al. (2004b) developed a more advanced method known as full-wave inversion of GPR data, which permits one to use all GPR information for 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 one 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., 2011b;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 illposed inverse problem caused by too many unknown parameters and the low sensitivity of GPR data with the electrical properties at deeper layers (Lambot et al., 2004a;Minet et al., 2011b).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., 2006Lambot et al., , 2009;;Jadoon et al., 2008Jadoon et al., , 2012;;Looms et al., 2008;Dagenbach et al., 2013).A joint inversion procedure based on integrated geophysical and hydrodynamic models was demonstrated to optimally estimate the soil unsaturated hydraulic parameters.However, in cases where the soil moisture profile slightly varies with the soil depth, some parameters suffered 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) and Walker et al. (2001) applied the extended Kalman filter (EKF) to retrieve the soil moisture profile from the nearsurface soil moisture measurements.Reichle et al. (2002), Das andMohanty (2006), De Lannoy et al. (2007), Huang et al. (2008), Crow et al. (2008) and Draper et al. (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) and 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) and 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 and petrophysical relationships worked 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 the 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 only surface soil moisture, the approach allows one 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 a soil moisture profile with better depth resolution in comparison 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.

Data assimilation procedure
Figure 1 shows the assimilation procedure that we used to update the state of the soil moisture profile using GPR data, which is outlined below: 1. Specify the initial soil moisture state and initial state ensemble based on the a priori knowledge on the moisture conditions.These values will work as the analysis state and analysis state ensemble at time t = 0.This proposed data assimilation scheme consists of 3 components: the hydrodynamic model, the observation operator and the MLEF data assimilation algorithm, which are presented in the following paragraphs.

The hydrodynamic model
In this study, the one-dimensional vertical flow in an homogeneous soil column was simulated by the water flow module of the 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 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) and petrophysical relationships worked as an observation operator to relate the soil moisture profile with the GPR data.We 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 antenna phase center is located at 7.5 cm above the antenna aperture.The frequency dependence of the phase center is accounted for by the transfer functions in the far-field antenna model (Lambot et al., 2004b;Jadoon et al., 2011).
The distance between the antenna aperture and medium was 26.5 cm, for which the far-field assumptions can be satisfied.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 freespace 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 ).For relating the soil moisture profile with the GPR data, petrophysical relationships are necessary to convert the soil moisture to these constitutive properties.In this study, the relationship between the soil volumetric water content (θ n ) and the dielectric permittivity was formulated by the model of Ledieu et al. (1986) as 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., 2011b), 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 as from Eq. ( 5).
For the fact that the assimilation procedure works with real numbers and the phase of the GPR data does not vary, given the fixed antenna height for all measurements, we used the absolute values of the complex GPR data in the frequency domain as the observation data.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 the far-field antenna model (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, 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 ensemble-based 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 one 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) and 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 is Gaussian random variables with zero mean and covariance matrix, Q, assumed to be the same form as the forecast error covariance one, P f , but with a smaller value: The ith column of the square-root N S × N E forecast error covariance matrix P 1/2 f at time t + 1 is calculated from the ensemble forecast as seen below: where p a i,t is the ith column of the square-root, N S ×N E analysis error covariance matrix P 1/2 a 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 determine the fact that these values are not available.

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 other words, minimizing the cost function given by 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. 11).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 ith column of the N O × N E observation perturbation matrix Z, z i is calculated as 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.( 14).Now, the optimization of w is implemented via the changing variable ξ .If the observation operator is linear and the changing variable (Eq.15) is employed, the solution of the optimization problem (Eq.11) 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) and ( 10) 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 80 equidistant elements.We analyzed the effects of the soil type on the data assimilation by considering 3 soil types: 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 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 the 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 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 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 retain water well.Even as the rainfall 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 frequencydomain synthetic GPR data were produced using the forward radar model.In the electromagnetic model, the soil moisture profile was modeled by 80-planar layers with an equal thickness of 1 cm under the air layer and above a lower half space.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 soil moisture profile, σ 2 = 0.013 2 ."Observed" GPR data were generated by  adding Gaussian random perturbations to the synthetic GPR data.The mean of the perturbations is set to zero and its standard deviation is equal to the variance components of the observation error covariance matrix.figure.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 with c, 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 is stronger as there is a 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 high frequency (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.We began to assimilate the GPR data to update the state of the soil moisture profile at a time of 80 h. Figure 4a and  b, respectively, present the synthetic and forecast soil moisture profiles at initial time (0 h) and at the 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 6 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.
To illustrate the effect of the update interval on the assimilation, we performed the data assimilation every 10, 30, and 60 h.

Assimilation results
Figure 5 compares the "true" and forecast soil moisture profiles of loamy sand, silt and clay at several time steps with an update interval of 10 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 80 h (1st assimilation cycle), the assimilation significantly improves the forecast soil moisture profile.The soil moisture profiles of the 3 soil types updated by the GPR data assimilation approach the "true" state much more closely than the open-loop profiles.At time 100 h (3th assimilation cycle), while the discrepancies between the assimilated and "true" soil moisture profile are clearly observed for the clay soil, they are very small for the loamy sand and 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 of 170 h (10th assimilation cycle), for the loamy sand and silt soils, the assimilated and "true" soil moisture profiles are 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 16 cm.For the deeper layers, the errors are higher (around 0.5-2 × 10 −2 cm 3 cm −3 ).At this time, 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 = 720 h), it is impossible to separate the "true" soil moisture profile of the loamy sand from the assimilated 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 4 × 10 −2 cm 3 cm −3 ), though they are smaller than at time 170 h.With respect to the clay soil, while no improvement is observed for the open-loop prediction, the agreement between the estimated soil moisture profile with assimilation and the synthetic one increases with absolute errors lower than 8 × 10 −4 cm 3 cm −3 .The figure also presents that the upper part of the soil moisture profile is better corrected than its lower counterpart, which is more 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 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 10 h.Comparing Fig. 6a and b, 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 openloop prediction for this soil type ranges from 1.6 × 10 −3 to 1.3×10 −2 larger than that of the GPR data assimilation.This quantity is 3.8-8 × 10 −2 for the silt, and 8-9.7 × 10 −2 for the clay 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 < 2 × 10 −3 ) at time step 610 h (for the open-loop prediction) and 180 h (for the GPR assimilation).For the silt soil, the RMSE for the open-loop prediction steadily decreases from 8.7 × 10 −2 to 3.8 × 10 −2 , while the GPR data assimilation quickly reduces and is approximately equal to the "true" state at time 160 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 1×10 −1 to 2.2×10 −2 and continue to decrease to 5 × 10 −3 at time 210 h.From time 210 h to the end of the simulation period, the RMSE reduce slightly from 5 to 3×10 −3 for the fact that the corrected soil moisture at the upper part 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 intervals of 10, 30 and 60 h.It is worth noting that, for each soil type, the RMSE of all update intervals 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 soil moisture profile obtained by GPR data assimilation for all update intervals is relatively small compared to the open-loop prediction, indicating that the wrong initial condition can be effectively corrected by assimilating GPR data even with the update interval up to 2.5 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 with decreasing update interval.However, this dependence is different for different soil types.For the loamy sand and silt soils, the obtained results indicate that for this "wrong" initial condition problem, decreasing update interval does not much improve the effectiveness of the GPR data assimilation.The difference of the RMSE obtained from the update intervals of 10 and 60 h varies from 5 × 10 −4 to 3 × 10 −3 for loamy sand, and from 7 × 10 −4 to 2 × 10 −3 for silt soils.The stronger effects of the update interval on the GPR data assimilation is found for the clay soil.The difference of the RMSE between update intervals of 10 and 60 h varies in the range 5-6 × 10 −3 .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 are the most sensitive layers with the GPR signal, is very large compared to that observed for the loamy sand and silt soil.As a result, decreasing 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 levels off for the clay soil.Consequently, while decreasing update interval insignificantly impacts on the effectiveness of the GPR data assimilation for the loamy sand, it better corrects 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 for all soil types, the RMSE of the prediction without the GPR data assimilation is much larger than that with the GPR data assimilation.
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.
With respect to the relationship between the assimilation and update interval, the obtained results show that for the wrong initial condition problem, decreasing update interval slightly improves the updated soil moisture profile for the loamy sand and silt soils with the difference of the RMSE between update intervals of 10 and 60 h lower than 3 × 10 −3 , 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 intervals of 6 × 10 −3 .
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.Our next research will concentrate on the assimilation of GPR data to simultaneously update both soil hydraulic parameters and state variables by combing the MLEF algorithm with the state augmentation technique.Moreover, the study ignored the frequency dependence of the electrical conductivity and dielectric permittivity, which cause errors for the observation operator.In order to reduce these errors, the frequency dependence of the electrical properties should be taken into account in realistic applications.This problem can be solved by using the conceptual dielectric mixing models (e.g., complex refractive index model (CRIM)) to relate the effective complex permittivity with the soil moisture, instead of the empirical formulas (Eqs.6 and 7) like in this synthetic study.The soil in the dielectric mixing models is considered to be a mixture of three components, namely, air, water and soil matrix.The frequency dependence of the soil electrical properties is accounted for via the free water component, which varies with frequency by the Debye's equation (Debye, 1929).

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

FigureFig. 4 .Fig. 4 .
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).The measurement errors are represented by the noises in the

Fig. 5 .Fig. 5 .
Fig. 5. Comparison of the open loop and data assimilation with the synthetic soil moisture profile at times 80, 100, 170, and 720-hour.The assimilation was performed every 10 hours.The three soil types were accounted for: (a, d, g, j) loamy sand, (b, e, h, k) silt and (c, f, i, l) clay.

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

Fig. 7 .
Fig. 7.The RMSE of the soil moisture profile obtained by the assimilation of GPR data with intervals, namely 10, 30 and 60 hours.The three soil types were compared, namely, (a) loamy s (c) clay.

Fig. 7 .
Fig. 7.The RMSE of the soil moisture profile obtained by the assimilation of GPR data with different update intervals, namely 10, 30 and 60 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.