A Polarimetric First-Order Model of Soil Moisture Effects on the DInSAR Coherence

Changes in soil moisture between two radar acquisitions can impact the observed coherence in differential interferometry: both coherence magnitude |γ| and phase φ are affected. The influence on the latter potentially biases the estimation of deformations. These effects have been found to be variable in magnitude and sign, as well as dependent on polarization, as opposed to predictions by existing models. Such diversity can be explained when the soil is modelled as a half-space with spatially varying dielectric properties and a rough interface. The first-order perturbative solution achieves–upon calibration with airborne L band data–median correlations ρ at HH polarization of 0.77 for the phase φ, of 0.50 for |γ|, and for the phase triplets Ξ of 0.56. The predictions are sensitive to the choice of dielectric mixing model, in particular the absorptive properties; the differences between the mixing models are found to be partially compensatable by varying the relative importance of surface and volume scattering. However, for half of the agricultural fields the Hallikainen mixing model cannot reproduce the observed sensitivities of the phase to soil moisture. In addition, the first-order expansion does not predict any impact on the HV coherence, which is however empirically found to display similar sensitivities to soil moisture as the co-pol channels HH and VV. These results indicate that the first-order solution, while not able to reproduce all observed phenomena, can capture some of the more salient patterns of the effect of Remote Sens. 2015, 7 7572 soil moisture changes on the HH and VV DInSAR signals. Hence it may prove useful in separating the deformations from the moisture signals, thus yielding improved displacement estimates or new ways for inferring soil moisture.

soil moisture changes on the HH and VV DInSAR signals.Hence it may prove useful in separating the deformations from the moisture signals, thus yielding improved displacement estimates or new ways for inferring soil moisture.

Introduction
The techniques repeat pass radar interferometry and differential interferometric SAR (DInSAR) are routinely applied to derive deformations of the Earth's surface, digital elevation models or vegetation heights [1].A possible influence of soil moisture changes on the observables in DInSAR has been surmised for decades [2].Such an impact can impair the estimation of the aforementioned parameters but might also open the possibility for soil moisture retrieval by this technique.Since the first conjecture of the presence of such effects on the phase φ and the coherence |γ| [2], several explanations about the origins have been advanced, and both observational studies (e.g., [3,4]) and dedicated experiments (e.g., by Rudant et al. [5], Morrison et al. [6]) have been conducted.The empirical evidence tended to be inconclusive and some results were inconsistent.This lack of agreement might be due to different study designs-e.g., radar wavelength and soil properties-and shows that the prevalence and size of these effects on the observables remains poorly understood.
This knowledge gap also encompasses the predictive power of the proposed explanations.Their applicability has recently been analysed empirically at L band by Zwieback et al. [7]: for the two airborne campaigns studied, the inferred dependence of φ on soil moisture changes ∆m v was not consistent with the penetration depth hypothesis [8] or soil swelling, but only with a dielectric volume scattering mechanism.More than 20 years ago the possibility of volume scattering within the soil was considered by Gabriel et al. [2], but they preferred to attribute the observed phase patterns to deformations as the soils studied were prone to swelling; the observed sign of φ was mostly consistent with such swelling.The sign of the phase change was also studied by Rudant et al. [5] in two laboratory experiments: they interpreted the phase signal not due to deformation as one of dielectric origin.A numerical model of a heterogeneous soil (which gives rise to scattering) based directly on Maxwell's equations predicts this sign of the phase dependence [9], as does the analytical approach by De Zan et al. [10].They considered volume scatterers embedded in a soil bounded by a smooth surface: in the first-order scattering solution the scattered field could be expressed as the superposition of the return from all the particles, the phase and magnitude of each depending on its position and the bulk dielectric constant of the soil.The complex coherence γ followed by ensemble averaging assuming no correlation between the returns from different particles; it did not depend on the polarization.
At higher frequencies the opposite sign of the phase dependence has been observed in laboratory experiments.The homogeneous samples of soil (S-band measurements by Yin et al. [11]) and pure sand (C-band measurements by Morrison et al. [6]) were put into containers and the surface was smoothed.The observed phase dependence on soil moisture was found to by consistent with the Small Perturbation Model by Yin et al. [11].That is, the signals could be explained by scattering from a slightly rough surface.
As both surface scattering and sub-surface volume scattering seem the most promising hypotheses for explaining the effects on γ for different non-swelling soils and radar frequencies, we want to determine whether these can be combined in an analytical model based directly on Maxwell's equations [12].The relative importance of the two scattering contributions is intended to account for the variations in sign and magnitude of the observed soil moisture effects.To this end, we propose to extend De Zan et al. [10] and model the soil in a first order Born approximation, where the zero-order scenario corresponds to a flat interface and a homogeneous soil.The first order expansion incorporates the roughness of the surface and spatial variations in permittivity.The resulting parametric solutions allow us to make predictions about the DInSAR observables-in particular the phase φ-for all polarizations, and compare them to observations.The evaluation of the proposed model focusses on the consistency of the predictions for different observables and polarizations, as well as the patterns in the data that the model cannot explain.

Proposed Model
The proposed model combines scattering from two origins, i.e., spatial variations of the soil permittivity and the rough surface.The first-order impact of the permittivity fluctuations on the scattered E field will be studied first, followed by an analysis of the surface roughness.Subsequently we will combine these two and derive explicit expressions for the observables based on parameterizations of the permittivity fluctuations.We will represent E as a phasor throughout, and employ the exp iωt sign convention.

Permittivity Fluctuations
The soil is represented by the half space b in Figure 1a.Its dielectric properties are characterized (i) a mean scalar permittivity ε and (ii) a fluctuating part ε f prq with mean xε f prqy " 0. The mean permittivity ε is mainly governed by the amount of liquid water present in the soil but also e.g., the ice content and soil texture [13,14].This dependence can be described by dielectric mixing models such as those by Hallikainen et al. [15], Peplinski et al. [16].
Let a plane wave of angular frequency ω impinge upon this half space.Following Tsang et al. [17], let E and E b be the electric field above (region a, V a ) and inside (region b, V b ) the medium, respectively.These electric field vectors are each split into two parts: a perturbation term and an unperturbed component E 0 or E 0 b .These latter correspond to the E-field in the absence of permittivity fluctuations.The following set of integral equations is readily obtained from Maxwell's equations (see [17] Chapter 6.2.3): where Qprq " ω 2 µε f and the fields with superscript 0 denote the solutions in the absence of fluctuations, whereas G 0 ij represents the dyadic undisturbed Green's function for a source in region j and an evaluation point in region i.Equation (1a) relates the measurable field above the interface Eprq to the unperturbed field E 0 prq and to the field E b pr b q at all points inside the half-space.The latter has to fulfill the consistency relation of Equation ( 1b).
An approximate solution, which does not obey energy conservation, can be obtained using a perturbative approach.In the first-order Born approximation the unknown field E b on the right-hand side is replaced by the unperturbed one E 0 b .This field corresponds to a plane wave, whose wave vectors both above and below the interface are shown in Figure 1a, and whose magnitude is determined by the Fresnel transmission coefficient for refraction into the lower half-space.For a smooth interface in the far-field, using the plane-wave approximation for the phase, i.e., where W rm ´1s can be considered a constant [17].The rank one dyad P h projects any vector onto the horizontal polarization vector in medium b and maps it to the horizontal polarization vector in medium a; P v does likewise for the vertical polarization.The factors T h and T v are the two-way amplitude Fresnel coefficients (see Supplement, section 1), and e is the unit polarization vector of the incident field in region b.Let q m " rq HH m q ?2HV m q V V m s T be the scattering vector in the lexicographic basis [18] for acquisition m.It is a column vector as indicated by the transposition operation T and encodes the components of the scattered E-field Eprq for different transmit polarizations e [18].In interferometry two such acquisitions, e.g., m and n, are combined: the position of the antenna can change r m ‰ r n , as can the dielectric properties of the half-space b, ε m ‰ ε n .The interferogram is formed by cross-multiplying the scattering vectors of the two acquisitions to yield the interferometric covariance matrix rC mn s v " xq m q : n y [1].Its predicted value according to the model under the Born approximation follows by rearranging (1) using (2) and grouping the results according to both incident and scattered polarization components.Denoting the acquisition with an additional index, we get with where the dielectric fluctuations are encoded in xQ m pr b,m qQ n pr b,n q ‹ y.These variations represent small-scale variations in the dielectric constant, due to e.g., pores, clumps or stones.If the fluctuations are assumed to be uncorrelated in space, this reduces to Q mn pr b,m qδpr b,m ´rb,n q, which eliminates one of the integrals in Equation ( 3).The polarimetric behaviour is only governed by the transmissivity matrix rT s in Equation (3); a more general representation could be obtained using an effective covariance matrix rT s 1 encoding both transmission and scattering effects.

Surface Term
When the interface between regions a and b is not smooth, backscattering occurs even in the absence of inhomogeneities.The perturbative analysis can be extended to the surface roughness (characterized by a non-zero characteristic height η): the first order scattered field consists of (i) the volume contribution excited by the zero-order fields Equation (3); (ii) the first order contribution from the surface (Small Perturbation Model, SPM) [19].The latter approximation is valid for small roughness heights xk a η cos θ i y ! 1 and small slopes Bη By ! 1.The scattering vector q s,m at acquisition m is [18] q s,m " Z » -- where Z is a factor that does not depend on the dielectric constant, but on the roughness.Assuming that Z is constant and that the volume and surface contributions are uncorrelated, the total scattering covariance matrix rC mn s " rC mn s v `rC mn s s " rC mn s v `xq s,m q : s,n y.This assumption is commonly done in both polarimetric [20,21] and interferometric modelling [18,22] of surface and "volume" terms.Note that the surface contribution rC mn s s can in principle also be modelled differently, e.g., by the Kirchhoff Geometric Optics model, a constant or by semi-empirical formulae [23].

Propagation Phase
The integration kernel in the volume contribution of Equation (1) modulates both the magnitude and the phase of the return from any point in the halfspace.The propagation phase ϕ m in Equation (2), i.e., its phase term, depends on the two-way propagation of the plane wave in half space a via the real part of its wavevector k a,m and similarly for half space b via k b,m (see Figure 1a) where the second line uses the coordinate definitions of Figure 1b and the scenario is kept two dimensional for simplicity.y i is an internal variable determined by Snell's law, and w m " " sin 2 θ i ` 2 ´aε m ´sin 2 θ i ¯ı0.5 -it corresponds to the index of refraction in the absence of absorption-such that As described in more detail in the supplement (section 2, Fig. S1), the interferometric phase of a point scatterer is given by the difference in the propagation phases and evaluates to where the second line is obtained by linearization with respect to the position of the antenna (F B " k 0 B 2 K R ´1 ! 1, where B K is the projection of the baseline vector perpendicular to the line of sight).φmn is the part due to different background dielectric constants, whereas φmn is a correction due to an antenna offset in height ∆H and in the horizontal direction ∆y a .
For scatterers close to the interface it is expedient to simplify the expressions even further by linearization with respect to the position.Particularly relevant are the changes with respect to z at constant range R m -the vertical interferometric wavenumber `Bφ Bz ˘Rm -and with respect to R m at constant height-the range wavenumber ´Bφ BRm ¯z.The total differentials evaluated just below the surface-the antenna positions being held constant-of the two terms in Equation ( 9) are derived in the supplement: The interferometric wavenumbers follow from Equations ( 10) and (11), cf. also Dall [24] for the vertical one: Note that the range wavenumber does not depend on the permittivity and is identical to the one above the interface [25], thus opening the possibility to do range spectral filtering [26].The linear approximation for a difference in φ of two scatterers separated by ∆z is appropriate if F z " k 0 p∆zq 2 R ´1 ! 1.
In astronomic interferometry the bispectrum is analysed when no adequate phase reference can be established [28].It is defined as a triple product of complex visibilities [29], which correspond to the interferometric covariances rC mn s.We propose to adapt this bispectrum to the bicoherence Γ mno in InSAR by forming the analogous triple product of the normalized covariance, i.e., the coherence Its argument is the phase triplet [10] (or closure phase) Ξ mno " φ mn `φno ´φmo .For N acquisitions there are only `n´1 2 ˘independent triplets, as opposed to `n 2 ˘phases [28].The lost information consists of the phase offset of each slave acquisition [30], which implies that these phase triplets are insensitive to deformations, atmospheric influences and general phase offsets.

Parametric Solutions
In order to arrive at closed-form expressions, we propose to represent the dielectric fluctuations as simple functions of depth.The magnitude of these fluctuations consequently does not depend on the horizontal position but only on z.Using these representations, along with the linearized phase terms of Equations ( 12) and ( 13) in Equation ( 1) will yield such analytic expressions.Let δ 0 be a representative penetration depth and z 1 " zδ ´1 0 (k 1 " kδ 0 ) the dimensionless depth (a dimensionless wavenumber).The ensemble average in Equation (3) will be assumed real and parameterized as where apz 1 q is a polynomial of degree A, which we represent in the Laguerre basis.Note that this assumes (i) invariance of the dielectric fluctuations between acquisitions, i.e., only the mean ε changes; and (ii) a lack of spatial correlation.These assumptions mean that these fluctuations are like white noise added to the background dielectric constant.For simplicity only the cases A " 0, i.e., apz 1 q " a 0 , and A " 1, i.e., apz 1 q " a 0 `a1 p1 ´z1 q, will be considered [31].These two cases are illustrated in Figure 2. The representation as a zero-order polynomial (A " 0) implies constant heterogeneities with depth, which is shown in the left panel.The linear dependence with depth (A " 1) is illustrated for a 1 ą 0 in the right panel: the soil becomes increasingly homogeneous as Figure 2. Illustration of the depth dependence of the dielectric fluctuations apz 1 q of Equation ( 16) for the parametric functions with A " 0 (constant, top panel) and A " 1 (linear dependence, bottom panel).In the left-hand diagrams, the graininess of the soil corresponds to the magnitude of the dielectric fluctuations, whose spatial variation is shown on the right.(a) A " 0; (b) A " 1.
The assumption of a real and acquisition-invariant profile of the spatially uncorrelated dielectric fluctuations renders Equation (3) similar to the standard first-order scattering model employed in interferometry, with apz 1 q corresponding to the structure function, as employed-also as a series expansion-in coherence tomography [27].
In the following it will be assumed that the B K term in Equation ( 12) is negligible, i.e., that the phase contribution due to different elevations is much smaller than π, k 0 B K R ´1δ 0 ! 1.We also assume that range spectral filtering [26] has been applied so that the range wavenumber vanishes.These assumptions are made because the data used in this study were acquired from the same track.In this case the exponential term in Equation (3) reduces to 2i ´k1 b,m ´k1 ‹ b,n ¯¨ẑ z 1 where ẑ is a unit vector in the vertical direction; this is the same functional form as in De Zan et al. [10].The integration in the horizontal direction (which should include the antenna pattern) can be absorbed into the normalization constant W 1 , as it cancels when the coherence Equation ( 14) is formed.Owing to the lack of spatial correlation of Q as parameterized by apz 1 q in Equation ( 16), one of the integrals in Equation (3) drops out For A " 0 the integral in Equation ( 17) yields whereas A " 1 gives The terms a 0 and a 1 capture not only the spatial variability of the dielectric fluctuations but also their absolute size.This size is usually not known in practice, and also presumably difficult to determine in-situ, which is discussed in more detail in Section 5.3.Similarly, the absolute value of the surface contribution depends on the roughness properties, which are generally not known a-priori.For practical reasons, we thus introduce a weighting parameter, the volume-to-surface ratio f .It expresses the relative importance of the two contributions (the absolute values are not necessary if only the coherence is looked at) at a reference moisture value: where f 0,s and f 0,v are normalization magnitudes, corresponding to the element of rC mm s s or rC mm s v at a reference moisture state in a specified polarization.As the scaling of rC mn s is immaterial for the evaluation of γ, a 0 will be set to one.The case of A " 0 and f Ñ 8 (pure volume) corresponds to the scenario studied by De Zan et al. [10]: its coherence is independent of frequency (and thus the free space wavelength) for frequency-independent ε as the half space has no intrinsic length scale.
The simulated impact of a varying volume-to-surface contribution f on the coherences is illustrated in Figure 3a.For a fixed change in soil moisture of ∆m v = 0.1 m 3 ¨m´3 , an increase in f makes the absolute values |γ| decrease (i.e., the soil decorrelates more) and the magnitude of the phases grow.This increase in sensitivity is due to the fact that the sensitivity of the phase to soil moisture changes is much smaller (and of opposite sign) for the surface term than for the bulk volume.A similar difference in sensitivity is present for |γ|, which is predicted to be equal to one for the surface, which also implies vanishing phase triplets.The volume term, by contrast, leads to decorrelation: |γ| is predicted in Figure 3b to decrease monotonically as the absolute change in soil moisture |∆m v | increases, whereas the magnitude of the interferometric phase is expected to increase.This relation is non-linear: the phase φ saturates as ∆m v increases, and it is not a function of ∆m v alone but depends on both values of m v .This is illustrated for two different dielectric mixing models in Figure 3c,d, whereas Figure 3e,f shows the impact of varying a 1 (A " 1).v and m o v , using A " 0 (solid line) and A " 1 (a 1 " ´0.7, dashed line).

Calibration and Model Assessment
In order to assess the accuracy and validity of the model, we compare its predictions of the radar observables to airborne radar observations.These model predictions are based on the soil moisture that was measured in-situ at the time of the radar acquisitions.As the model contains unknown parameters, these will be estimated in a calibration step.The calibrated model can subsequently be directly compared to the radar measurements.

Study Site
The available data were acquired within the Canadian Experiment for Soil Moisture in 2010 [32] and cover the time period from 2-14 June 2010.The test site near Kenaston, Saskatchewan, Canada (51 ˝30 1 N, 106 ˝18 1 W) covers an area of around 270 km 2 .Rainfed agricultural fields, pastures and grassland predominate, but at least 1.5% of the area is covered with open water surfaces.This percentage is likely to have been even larger due to the wet conditions before and during the campaign [33].The majority of the agricultural fields had been tilled prior to the measurements and several were covered with crop residues to varying degrees [32].

Soil Moisture Data
Environment Canada operates a network of permanent stations where soil moisture is recorded by Stevens Hydra Probe II sensors at different depths [33], from which we use the data measured by those installed vertically at 0-5 cm depth.Due to battery failure, the presence of standing water in some fields and limited image extent, only the data of ten fields are available, of which we use the same seven as in Zwieback et al. [7]; the remainder are put aside for future inversion studies.The sensor outputs are converted to values of volumetric soil moisture m v (m 3 ¨m´3 ) using a calibration procedure based on soil texture.
The m v time series of these sensors show similar temporal patterns, but both the mean values and dynamic ranges differ, the latter by about a factor of 10.An illustration of this diversity is given in Figure 4 for two different fields.

Radar Data and Processing
The fully polarimetric (quadpol) UAVSAR data (L-band: λ " 0.24 m) have a resolution of 1.7 m (0.8 m) in range (azimuth) [34].They were acquired in irregular intervals between one and three days with a nominal baseline of 0 m and from a height of 13 km.The incidence angle ranges between 22 and 66 ˝from the near to the far range of the image [33].Using all possible combinations of the six SLC images, 15 interferograms were formed.As in Zwieback et al. [7], a rectangle (extent: 50 m in range, 100 m in azimuth) was manually delineated in the immediate surroundings of the location of each probe; it was chosen so as to maximize the homogeneity within this region of interest (ROI) and avoid pools of water.At each of the four corners of these rectangles, the interferograms were averaged (number of looks L = 172, area: 20 m 2 ) to estimate the coherence.These subROIs are denoted by the field number (three digits) and a letter a-d for each of the corners.For each of the fields, we reference the phase with respect to the closest stable scatterer as in Zwieback et al. [7], where it was observed for the same data that the results are not sensitive (less than 10% change on average) to the choice of reference target.
As the flight tracks during the acquisition deviated from the nominal ones-the UAVSAR team aims at a track accuracy of 10 m and achieved one of 5 m in this campaign-the interferometric phase depends on the height, thus introducing an additional source of uncertainty.In the worst case scenario for an incidence angle of θ = 30 ˝, the height of ambiguity H a « 190 m (the height corresponding to a phase difference of one cycle).For a height difference between the field and the phase reference of 5 m, we thus expect a maximum phase error of 10 ˝, which is about one order of magnitude smaller than the maximum phase differences expected from the model or observed in the data.An example of these observed phases is shown in Figure 4: as the soil becomes more moist between 6 and 9 June, the measured phases with respect to a drier acquisition on 5 June take on values of around 70 ˝or 3π 8 rad.They subsequently decline, and so does the soil moisture, whereas the absolute coherences increase again.The radar observations in field 206 (Figure 4b) exhibit similar temporal behaviour, but the magnitude is smaller by a factor of 2. The soil moisture observations also display a smaller dynamic range (by about 75%).

Rationale
In order to assess the validity of the proposed model of Equation ( 20), we compare its predictions with these L-band measurements.The predictions are governed by two parameters that are not known a priori: f and a 1 .Assuming the model can capture all the influences on the coherence-thus neglecting phenomena such as vegetation dynamics-we estimate these via minimization of misfit functions: these measure the difference between the predicted observables (based on prescribed soil moisture) and the measured ones.This estimation is done for each subROI separately.The overall aim of the calibration is to first estimate the parameters, and then assess the quality of fit.

Minimization of Misfit
We focus on the two normalized observables of Section 2.4, i.e., γ and Γ, and separate the magnitude and argument of these complex quantities.For the coherence we propose to minimize a generalized root mean square (RMS) deviation µ γ pζq, where ζ weights the relative importance of magnitude and argument: where the weights w φ mn are the reciprocals of the estimated variances (based on the Cramer-Rao bound of the Gaussian speckle model, see Bamler and Hartl [1]) and are normalized such that ř w mn " 1.This normalization condition is also enforced for the uniform weights w |γ| mn .The predicted values are denoted by hats and depend on the parameters f and a 1 , as well as on the prescribed soil moisture values.Note that the cos dependence of the phase misfit is not susceptible to phase wrapping and that the metaparameter ζ weights the influence of φ and γ, with ζ " 0.5 corresponding to a roughly equal importance.
The bicoherence misfit function is structured along the same lines; the weights w Ξ m 0 ,n,o are obtained assuming that the deviations of the phases φ mn are uncorrelated for different acquisitions, and the w Γ m 0 ,n,o are again uniform.The first acquisition m 0 in the bicoherence is held fixed as otherwise redundancy would be introduced [28].

Evaluation
In order to assess the impact of different parameterizations, such as the dependence on the dielectric mixing models, we employ bootstrap resampling: Bias Corrected and accelerated (BCa) confidence intervals are derived using case resampling [35].These quantities are estimated by taking random subsamples from the data, with which one can then approximate the distribution of the quantity of interest.These are the difference of the estimates of a parameter (e.g., the volume-to-ground ratio f ) between two parameterizations.
For any parameterization of the model, the minimized value of the misfit function µ is an indicator of the discrepancy between model and observations.Apart from this value, the quality of fit will also be assessed by the correlation ρ and the sensitivity ratio s. ρ is the Pearson's correlation coefficient between predictions and measurements, whereas the sensitivity ratio spζq is given by the ratio of the standard deviation of the predictions to that of the data.A value s >1 thus corresponds a larger spread of the predicted values than the observed ones, whereas s <1 indicates a comparatively larger scatter in the data.This can be due to a larger m v sensitivity of the observations than predicted by the model, but also due to influences not represented by the model.In particular, it is expected that a poor model fit will lead to small values of s upon calibration.

Examples
Such differences between the measured DInSAR observables and their predictions are found in the data, e.g., in the examples shown in Figure 5. Panel (a) shows the predicted φ V V of both a volume-only (f Ñ 8) and the calibrated version (f estimated to yield f using ζ " 0).For both versions, there is an approximately linear relation between predictions and observations of the correct sign but the observed values are greater in magnitude than the predictions (s " 0.56).These simulated values of the two parameterizations cannot be distinguished as f " 1.Such a difference does occur for φ HH in Figure 5b, where the sensitivity ratio s " 2.93 is reduced, upon calibration, to 0.76.Similar variations in the sensitivity are also found for the phase triplets Ξ, e.g., in Figure 5c,d).

Calibration with A " 0
The calibration results of all the fields for the simple case of A " 0-i.e., the fluctuation term Q is independent of depth-are summarized for the HH polarization in Taylor plots in Figure 6.The location of each point encodes the correlation ρ via its angle and the relevant misfit µ via its distance from the origin; the colour depends on the field according to the topmost legend and the sensitivity ratio s is represented by the symbol as indicated by the legend underneath.The three panels (a) to (c) show the evaluation results for the phase φ, i.e., the correlation and sensitivity ratio between predicted and measured as well as the phase misfit µ γ p0q.The evaluation measures of the same models but applied to the coherence magnitude |γ| are displayed underneath.In panel (a) the Peplinski-based model was calibrated on phase data alone, i.e., ζ " 0. The subROIs results tend to cluster with respect to the field, achieving a median correlation ρ of 0.71.Typical misfits µp0q vary between 0.1 and 0.25, with only field 201 (light blue) showing worse fits.For this field, the model cannot replicate the interferometric observables based on the measured soil moisture data.This is even more evident for the coherence |γ| in panel (d), as negative correlations are obtained.This field shows considerably larger variations in the interferometric observables than the model can predict: even in the pure volume case, the sensitivity ratio s ă 1.Such low values of s for the volume-only case are not observed in the other fields.The exceptional properties of this field 201 are also reflected in the results obtained using ζ " 0.5 (roughly equal weight on φ and |γ|).As seen in Figure 6b, the overall fits and correlations are similar to the phase-only case using ζ " 0. The phase-only calibration ζ " 0 using the Hallikainen mixing model contrasts with the Peplinski model in two important aspects.Firstly, according to Figure 7d the model cannot reproduce the magnitude of the observed phases, i.e., s ă 1.Secondly, the median correlation of the phase increases to 0.77, and the one of the coherence magnitude from 0.5 to 0.58.There is no clear trend in the misfits µ γ p0q.They are significantly different between the two mixing models (at a significance level α " 0.05) in 45 %, as shown in Figure 8a.Non-significant differences occur e.g., for fields 109 and 136, corresponding to confidence intervals that overlap with µ γ p0q Pep ´µγ p0q Hal " 0, i.e., equal misfit of the two mixing models.The Taylor plots for the VV phase-only calibration (i.e., ζ " 0) in Figure 9 reveal similar patterns for the Pep and Hal mixing model: more than 90 % of the correlations are between 0.3 and 0.95 (for φ: median 0.78 and 0.70 with Hal and Pep, respectively).The two versions of the model cannot always reproduce the observed magnitudes of the interferometric observables.As seen in Figure 7b,e, this occurs for the same fields as in the HH polarizations.The misfits µ γ p0q ă 0.4, except for field 201, with significant differences between Pep and Hal occurring in 36% of the cases (Figure 8b).The volume-to-ground ratios f estimated independently for HH and VV phases do not show a consistent offset in Figure 8c.Taylor plots for the phase φ and the coherence γ at VV polarization; cf. Figure 6 for a description of these plots.(a The model of Equation (20) predicts zero contribution in the cross-pol channel from both the surface and the volume term.As the observed signals are similar to the ones found in HH and VV [7], we study them using a volume-only model.In this approach the surface contribution is not considered (f Ñ 8), and the volume model of Equation ( 17) is adopted (using the appropriate Fresnel terms).Its correlations (median: 0.73) and misfits for the HV phase are plotted against those of the calibrated HH model applied to φ HH in Figure 7c,f.The correlations tend to be larger for the HH channel-however more than 85% exceed 0.5 in the HV channel-and the misfits µp0q smaller.

Calibration with A " 1
The inclusion of an additional term in the expansion of the dielectric fluctuations, i.e., letting A " 1, always reduces the misfit of the best-fitting model in both HH (Figure 10a) and VV (Figure 10b) channels.This reduction is less than 5% for all samples except two.Note that the misfit cannot increase as the A " 1 model includes the A " 0 model as a special case.The correlations ρ increase for more than 90% of the samples (see Figure 10c,d).

Bicoherence
The calibration procedure based on the bicoherence is achieved by minimizing the misfit µ Γ of Equation (22); the relevant Taylor plots for HH are shown in Figure 11.The subfigures (a) and (b) reveal that for A " 0, the correlations ρ are typically between 0.3 and 0.75, with median values of 0.56 and 0.58 for the Hal and Pep mixing model, respectively.The correlations and also the misfits display as much intra-field as inter-field variability.Panel (c) reveals that when a 1 is estimated as well (A " 1), the correlations tend to increase (median: 0.56) along with the sensitivity ratios.The great majority of these sensitivity ratios s are less than one, even if the volume-to-surface ratio f Ñ 8, in which case only 22% exceed 1 for the Hallikainen model.
The correlations obtained at VV (see Figure 12) are comparable, with a median value of 0.59 for both mixing models.The misfits, by contrast, are smaller by around 30%.

Discussion
The comparison of the measured observables with the measured soil moisture changes and the model predictions has revealed certain agreements and inconsistencies.These will subsequently be discussed and assessed with respect to the uncertainties of the radar observations and the soil moisture measurements.The differences between the polarizations will be analysed first, followed by the impact of different parameterizations of the model and possible extensions or generalizations.

Polarimetry
In the proposed model, both the term due to the dielectric heterogeneities Equation (3) and the one due to the surface scattering Equation ( 5) predict zero backscattering in the cross-pol channel HV.This prediction, which is a general feature of such first-order expansions of dielectrically isotropic scattering scenarios [23], is known to fail in practice, and has thus led to the introduction of both more sophisticated analytical and simpler semi-empirical models [36,37].The discrepancy between predictions and observations is found to also apply to the interferometric coherence: soil moisture changes impact the DInSAR signals at HV [7].The measured phases and their observed dependence on soil moisture changes cannot be explained by noise or the uncertainty due to the phase referencing.Rather, they can also be described by the volume scattering of Equation (3) provided that the HV channel is treated like the co-polar ones.The correlations ρ between predicted and observed phases obtained in this way (see Figure 7) are of comparable size to the ones in the co-polar channels.This similarity suggests that it is the changing propagation within the soil that is crucial [10] rather than the more local properties that give rise to scattering.For many practical purposes it might not be necessary to introduce more complex models (such as higher-order expansions); rather, this semi-empirical extension of treating the HV like the HH channel might suffice.
These values of the phase correlations ρ at HH and VV indicate that the model can capture the most salient features of the observed phases.These are thus considered to be larger than decorrelation noise and the uncertainties due to phase referencing (see Section 3.1) and their sign and general dependence on moisture changes can be described by the model.However, these ρ values do not reveal the mismatch between the sensitivities to soil moisture changes of the model and the data.For four fields the observed ones tend to exceed the maximum sensitivities obtainable using the Hal mixing model: these are achieved with the volume-only version as f Ñ 8, see Figure 3a.Such a case is displayed in Figure 5a: as the sensitivity cannot be further increased and f is the only parameter that can be adapted, the optimized and the volume-only model coincide.The tendency of a reduction of sensitivity with decreasing f becomes evident in one of the counterexamples, shown in Figure 5b.
Such a discrepancy in the magnitude between the modelled and the observed impact of soil moisture can be caused by inconsistent definitions and calibrations of the input data (e.g., the soil moisture), or by deficiencies in the model itself, cf.Section 5.2.In particular in this data set, the dynamic ranges of the measured soil moisture time series differ by about one order of magnitude despite similar soil textures and the absence of vegetation [32].Manual soil moisture sampling on some of the days indicated a smaller inter-field variability of the soil moisture dynamics than the observations by the permanent probes used in this study.The interferometric observables show a smaller variability in magnitude than the in-situ measurements, cf. Figure 4 and Zwieback et al. [7].More generally, such discrepancies are not uncommon in backscatter-based estimation of soil moisture: this is why the correlation coefficient remains one of the most popular accuracy metrics [38], and the differences in both mean and dynamic range have to be considered when comparing and merging these data, as well as for model integration [39,40].
For many practical purposes, the link between soil moisture and the phase will be more important than the one with the coherence |γ|, not least because there are many processes that lead to decorrelation over time [41,42].The data set used in this study is characterized by a repeat period of 1-3 days, which is smaller than that of most satellite missions.The impact of the soil moisture on the coherence is evident in the timeseries of Figure 4, where the coherences increase again after a rain event has temporarily led to wet conditions and decorrelation.Except for field 201, for which negative correlations ρ between predicted and modelled |γ| are obtained, the ρ are similar to those observed for the phase φ, see Figures 6 and 9.The sensitivity of these predictions to soil moisture is, like for φ, smaller than expected (s ă 1) in certain cases.
Such a small sensitivity is also observed for the phase triplets.As for the comparison with φ, the correlations with the observations vary typically between 0.3 and 0.9.These larger sensitivities observed in the data correspond to negligible surface contributions (as in Figure 5d), as the surface term leads to zero closure phases for all values of m v , resulting in identical predictions for the volume-only and the calibrated models.

Parameterization and Uncertainties of Soil Moisture Dependence
The dependence on the parameterization is particularly pronounced for the link between the dielectric constant and m v , i.e., the mixing model and the calibration of the soil moisture data.The difference between the Hal and the Pep mixing model is also evidenced by the internal sensitivity analysis of Figure 3: compared to the Hal model, the Pep model predicts (i) a higher sensitivity of φ to ∆m v ; (ii) increased decorrelation for fixed |∆m v |; and (iii) in general larger magnitudes of the phase triplets Ξ.The main difference between these two models is the loss tangent, and thus the absorption, both of which are smaller for the Pep model [16].For a fixed value of m v , the wave penetrates deeper into a soil whose dielectric properties are governed by the Pep model.There is consequently a proportionally larger amount of scattering from heterogeneities deeper in the soil; when m v changes, these heterogeneities add contributions with a larger change in propagation phase than the ones closer to the surface.This leads to (i) larger changes in the interferometric phase φ (as the deeper parts have a larger contribution); (ii) lower values of |γ| (due to the larger effective diversity in propagation phases); and (iii) larger phase triplets (as they vanish for surface-scattering and increase in magnitude with decreasing loss tangent).
The overall similarity in the results despite the difference in the sensitivities can be partially explained by the compensating effect of a varying f .A relative increase in volume scattering acts similarly to a decrease in absorption: it also leads to larger phase sensitivities, more decorrelation and larger phase triplets.Empirically, the calibrated f are consistently larger for the Hal model than for Pep.However, the compensation of the difference in the loss tangent is not perfect, as more than 40% of the samples show significant-at α " 0.05-differences in the misfit (Figure 8a,b), i.e., the differences are larger than the uncertainties.There is no clear tendency towards any of the two, but the majority of fields can be modelled more accurately by the Hal model, which on average also achieves higher correlations.This points towards more complex disparities between the two mixing models and the predictions obtained with them.The deviations in the latter, such as the nonlinear components, are seen in Figure 3c,d to increase with differences in soil moisture.The latter figure indicates the differences in the predictions for Ξ, the impact of which on the misfit is seen to be similar to the one of φ in Figures 11 and 12.
The small impact of extending the expansion of the fluctuations to A " 1 on the quality of fit and the correlations is also consistent with such a balancing effect of f .A negative a 1 corresponds to an increase in the fluctuations Q mn pz 1 q with depth.Deeper parts of the soil thus exert a relatively larger influence on the scattering-similar to the Pep model compared to the Hal mixing model, or to an increase in f .The corresponding predictions of Figure 3e,f for φ and Ξ, respectively, reflect this similarity: the sensitivity to changes in soil moisture tends to increase when a 1 <1.Upon extension of the expansion of the fluctuations, the empirically determined misfits µ γ p0q always decrease.This decrease is always below 5%, and the size of the misfits remains stratified according the field.

Extensions and Underdetermination
The parameterization of the scattering contribution as a function of depth-in the present model it is encoded in Equation ( 16)-thus holds potential for improving the model fit.An analogous depth variation appears conceivable for the mean dielectric constant [10], but without direct observation (e.g., using tomographic techniques), there is little physical basis for either.Rather, it would presumably lead to an overparameterization, given the typically available data and the compensating effects mentioned above.
Arguably more critical issues from the point of view of structure of the model-as they concern barely justified assumptions implicit in Equation ( 16)-are (i) the postulated invariance of the permittivity fluctuations with the acquisition; and (ii) their lack of spatial correlation.Both were dictated mostly by practical reasoning: in particular (i) could-if relaxed-easily explain any interferometric phase φ as Q mn pzq would be represented by an arbitrary complex number for m ‰ n; similar observations apply to |γ| if m " n is also taken into account separately.The inclusion of spatial correlations (ii) would, in the isotropic case, mainly impact the overall backscatter from these fluctuations rather than their polarimetric properties (see [17], chapter 6.2.3).Anisotropic correlations, on the other hand, could possibly affect the latter, especially the HV component.As these fluctuations are somewhat germane to the model (they are not amenable to direct measurements except perhaps at the pore scale) and the relaxation of assumptions greatly increases the possible set of measurements explainable by the model, such an extension would suffer from underdetermination that could not easily by resolved.
Similar reasoning applies to the use of more general surface and volume scattering models.For instance, the use of discrete particles embedded in the soil as opposed to spatially varying dielectric properties, could yield similar results.Such discrete particles form the core of the model by De Zan et al. [10], but their properties (and thus e.g., the polarimetric characteristics of the scattered fields) were not specified.More general models would also permit a relaxation of certain limitations such as the lack of energy conservation, the small RMS roughness height or the lack of correlation between surface and volume scattering.The latter might be particularly critical when the expansions are extended to include more terms.Such higher-order models could also introduce multiple scattering, e.g., within the subsurface or of the surface [43].These higher-order contributions could also give rise to soil moisture effects in the interferometric returns from the soil.For instance, the double bounce between e.g., a vegetation stalk and a subsurface heterogeneity involves the same propagation phase within the soil as the direct return from the heterogeneity.It will thus introduce a similar dependence on soil moisture.

Conclusions
We introduced a fully polarimetric model that describes the impact of soil moisture changes on the measurements obtained using radar interferometry.It incorporates both volume scattering effects due to heterogeneities and a surface contribution in a first-order Born approximation of Maxwell's equations.This combination is intended to make the model capable of accounting for varying contributions of these two sources, which have previously only been studied separately.The volume contributions are found to be dominant in an airborne L-band data set, when the model is calibrated using the Hallikainen mixing model to estimate their relative importance.The estimation of the relative size of the volume contributions was hampered by uncertainties due to the calibration of the soil moisture probes.Regarding the temporal dynamics, the uncertainties in the interferometric phase were considered to be an order of magnitude smaller than the observed soil moisture signals.The predominantly positive correlations between these observations and predictions of the calibrated model indicate that the model can capture some of the most salient features of the impact of soil moisture changes on the coherence |γ| (median correlation at HH: 0.50), the phase φ (0.77), and the phase triplets Ξ (0.56).However, the presence of these effects in the HV channel cannot be reconciled with the model.The description of these effects at HV is thus an open question for which the inclusion of anistropic heterogeneities or higher order expansions seem promising.
At all polarizations the model predictions depend on the parameterization of the depth-dependence of the volume scattering and the soil moisture-dependence of the dielectric constant of the soil, in particular the absorptive properties.The variation of these parameters impacts the predictions in a way similar to a change of the relative importance of volume and surface scattering.This parameter can thus partially compensate these differences, and similar qualities of fit can be achieved.These findings indicate that even though the model predictions are sensitive to these parameterizations, data-driven calibration of the model might render it robust enough to be used in future studies.Possible future accuracy assessments might elucidate the impact of the radar frequency and the soil properties on the validity of the model.
The qualitative agreement, as expressed by the correlations, points towards the possible use of such a model to estimate soil moisture.It might also be instrumental in accounting for changes in soil moisture in the estimation of deformations-a task for which the phase triplets are most promising, as they are insensitive to movements.They can thus potentially be used to remove the soil moisture effects from the estimated deformations due to, for instance, volcanic activity or groundwater-related subsidence.

Figure 1 .
Figure 1.Description of the geometric setup: the circle at in the soil (halfspace V b ) is illuminated by the antenna at height H in V a .The coordinate vectors with respect to the origin (black circle at the surface) are shown in the left-hand panel, along with the wave vectors k in red.The panel on the right illustrates the coordinate definitions used in Section 2.3.(a) Geometry of the half-space; (b) Geometry for evaluating the phase.

Figure 3 .
Figure 3. Model predictions for different observables (HH polarization) using the Peplinski mixing model Pep [16] or the Hallikainen [15] model Hal.In panel a) the volume-to-ground ratio f varies, in the remaining ones the soil moisture varies.(a) Coherence values (Pep, A " 0) for m m v = 0.2 (0.3) m 3 ¨m´3 for the master (slave) and varying f ; (b) Coherence values (Pep, A " 0) for f = 100, a master m m v of 0.2 m 3 ¨m´3 and varying slave m n v ; (c) Interferometric phases φ mn (A " 0) for f = 25 and different m m v and slave m n v " m m v ∆m v ; (d) Interferometric phase triplets (A " 0) Ξ mno for f = 25 and m m v = 0.3 m 3 ¨m´3 and different m n v and m o v ; (e) Interferometric phases φ mn (Hal model) for volume-only and different m m v and slave m n v , using A " 0 (solid line) and A " 1 (a 1 " ´0.7, dashed line); (f) Interferometric phase triplets (Hal) Ξ mno for f = 25 and m m v = 0.3 m 3 ¨m´3 and different m n v and m o v , using A " 0 (solid line) and A " 1 (a 1 " ´0.7, dashed line).

Figure 4 .
Figure 4. Time series of measured phases φ and absolute coherences |γ| in two fields, where the master acquisition m is 5 June.The soil moisture values m v measured in the respective field are shown in the second panel; their dynamic ranges differ by about one order of magnitude.(a) ROI 307a; (b) ROI 206d.

Figure 7 .
Figure 7.The first two columns contain scatter plots of the phase correlation and misfit for and VV between volume-only and calibrated versions of the model.The third column shows a comparison between HV and HH fit using the Hallikainen mixing model.The colours of the symbols represent the agricultural fields as in Figure 6.(a) Phase slope HH: Peplinski model; (b) Phase slope VV: Peplinski model; (c) Correlation φ: HH and HV (volume-only); (d) Phase slope HH: Hallikainen model; (e) Phase slope VV: Hallikainen model; (f) Misfit φ: HH and HV (volume-only).

Figure 8 .
Figure 8. Bootstrap confidence intervals of differences in µ γ or f between two parameterizations.For each quantity and each ROI, the estimated difference is shown by the symbol (a diamond if it is significantly different from zero at α " 0.05, a circle otherwise), along with the 95% confidence intervals.(a) HH, A " 0, ζ " 0; (b) VV, A " 0, ζ " 0; (c) Hal, A " 0, ζ " 0.