Assessing the dynamics of soil salinity with time-lapse inversion of electromagnetic data guided by hydrological modelling

Irrigated agriculture is threatened by soil salinity in numerous arid and semi-arid areas of the world, chiefly caused by the use of highly salinity irrigation water, compounded by excessive evapotranspiration. Given this threat, efficient field assessment methods are needed to monitor the dynamics of soil salinity in salt-affected irrigated lands and evaluate the performance of management strategies. In this study, we report on the results of an irrigation experiment with the main objective of evaluating time-lapse inversion of electromagnetic induction (EMI) data and hydrological modelling in field assessment of soil salinity dynamics. Four experimental plots were established and irrigated 12 times during a 2-month period, with water at four different salinity levels (1, 4, 8 and 12 dS m⁻¹) using a drip irrigation system. Time-lapse apparent electrical conductivity (σ_a) data were collected four times during the experiment period using the CMD Mini-Explorer. Prior to inversion of time-lapse σ_a data, a numerical experiment was performed by 2D simulations of the water and solute infiltration and redistribution process in synthetic transects, generated by using the statistical distribution of the hydraulic properties in the study area. These simulations gave known spatio-temporal distribution of water contents and solute concentrations and thus of bulk electrical conductivity (σ_b), which in turn were used to obtain known structures of apparent electrical conductivity, σ_a. These synthetic distributions were used for a preliminary understanding of how the physical context may influence the EMI-based σ_a readings carried out in the monitored transects as well as being used to optimize the smoothing parameter to be used in the inversion of σ_a readings. With this prior information at hand, we inverted the time-lapse field σ_a data and interpreted the results in terms of concentration distributions over time. The proposed approach, using preliminary hydrological simulations to understand the potential role of the variability of the physical system to be monitored by EMI, may actually allow for a better choice of the inversion parameters and interpretation of EMI readings, thus increasing the potentiality of using the electromagnetic induction technique for rapid and non-invasive investigation of spatio-temporal variability in soil salinity over large areas.

1 Introduction

Soil salinization may be of a primary nature, when salt accumulation arises through pedogenetic processes, or of secondary origin, due either to abiotic factors such as excessive evaporation or seawater infiltration or resulting from human intervention, chiefly the use of saline water irrigation (Geeson et al., 2002). Approximately 20 % of irrigated land (45×10⁶ ha) that produces one-third of the world's food is salt-affected (Shirvastava and Kumar, 2015), and it is estimated that soil salinity affects 1×10⁶ ha in the European Union, mainly in the Mediterranean countries (Maréchal, et al., 2008).

Effective agricultural management in many areas relies on a good understanding of the effects of irrigation on the spatial and temporal variability of soil salinity (Coppola et al., 2015). However, it is very difficult to assess soil salinity on a management scale using traditional methods. Soil salinity is traditionally assessed by measuring the electrical conductivity of a saturated soil paste (EC_e) in the laboratory; however, this technique is labour-intensive, time-consuming and costly, given the large number of soil samples that need to be collected. Alternatively, time domain reflectometry (TDR), a non-destructive electromagnetic technique, can be used in the field for simultaneous measurements of water content, θ, and bulk electrical conductivity, σ_b (Coppola et al., 2011b). While the TDR method can provide accurate information from the local measurements, the measurement support volume of sensors is limited to a few centimetres; thus extension of the information to a large area can be problematic (Coppola et al., 2016; Gonçalves et al., 2017, Dragonetti et al., 2018).

The electromagnetic induction (EMI) technique is widely used as an alternative to traditional techniques for soil salinity assessment. It allows for rapid non-invasive, reliable and repeatable measurements at a smaller cost than traditional methods. This technique measures the soil electrical conductivity, which is primarily a function of soil salinity, soil texture, moisture content and cation exchange capacity; however, in a saline soil, the salinity is generally the dominant factor responsible for the spatio-temporal variability of soil electrical conductivity (Corwin and Lesch, 2005).

In the last few decades, EMI techniques have been used increasingly to estimate soil salinity from apparent electrical conductivity (σ_a) measurements (Lesch et al., 1995; Triantafilis et al., 2000; Corwin et al., 2006; Ganjegunte et al., 2014). σ_a is the weighted average of the soil electrical conductivity distribution in the soil volume. In order to obtain the depth distribution of σ_b from σ_a, a site-specific empirical calibration between σ_a and soil salinity measured at different depths can be applied by different approaches such as multiple regression (Triantafilis et al., 2000; Amezketa, 2006; Yao and Yang, 2010; Coppola et al., 2016), modelled coefficients (Slavich and Petterson, 1990), theoretical coefficients calculated with theoretical EMI depth response functions (Cook and Walker, 1992) or empirical-mathematical coefficients (Corwin and Rhoades, 1984).

Alternatively, to assess the distribution of σ_b with depth, σ_a data collected in the field can be modelled through an inversion process. Several inversion methods have been proposed to obtain the σ_b from the measured σ_a data, including the gradient-based inversion technique (Monteiro Santos, 2004; von Hebel et al., 2014; Schamper et al., 2012; Farzamian et al., 2015a) and probabilistic inversion (Jadoon et al., 2017; Moghadas, 2019; Shanahan et al., 2015). Recently, multi-coil EMI measurements and inversion algorithms have been widely used for mapping soil salinity and sodicity distribution in quasi 2D (e.g. Goff et al., 2014; Farzamian et al., 2019; Paz et al., 2019, 2020a, b) and 3D (e.g. Huang et al., 2017). However, the potential of this method in assessing temporal variability of soil salinity has not been fully explored. Several time-lapse inversion methods for direct-current resistivity methods have been developed. These include the ratio method (Daily et al., 1992), the difference inversion (LaBrecque and Yang, 2001) and, more recently, 4D space–time algorithms (Kim et al., 2009). A number of studies have also demonstrated how the use of time-lapse inversion algorithms can reduce the inversion artefacts (e.g. Hayley et al., 2011) and improve the quantitative investigations of geophysical monitoring (e.g. Farzamian et al., 2015b). However, the usefulness of the time-lapse inversion algorithms in modelling EMI data has not been attempted yet to assess soil salinity dynamics, and only few studies have been conducted to estimate soil water content changes (Huang et al., 2016; Whalley et al., 2017). In addition, a prerequisite for such an approach concerns the reliability of the inversion of the EMI result. In fact, inverting profile-integrated EMI data to obtain the vertical distribution of σ_b is an ill-posed problem, suffering from non-uniqueness (the problem has more than one solution) and instability (incomplete data and measurement errors can lead to large changes in the parameters; e.g. Tarantola, 1987). Ill-posedness is generally treated by regularizing the inverse solution. However, different regularization schemes and parameters can have a significant impact on the results (e.g. Dragonetti et al., 2018; Zare et al., 2020); thus, inversion results of EMI data are always affected by uncertainties, which can be minimized in case of prior information from the experiment.

In this direction, preliminary numerical simulations of the same hydrological processes to be monitored by an EMI sensor, by applying real boundary conditions measured during an EMI sensor monitoring campaign, may be especially helpful to figure out the response to be expected by the sensor and its variability in the space and time and may allow for a more rational choice of the EMI inversion parameterization. In other words, hydrological simulations may help provide a priori knowledge of “where the EMI inversion has to go”. In any case, this would require the actual distribution of the hydraulic properties along the transect or, more in general, in the field to be monitored by the EMI. One can immediately realize that this is quite utopian, especially when the area to be monitored is relatively large (as previously recalled in the case of EMI measurements). By contrast, it is more common that, for the study area, one has the statistical distribution of hydraulic properties available. Thus, the statistical distribution of the hydraulic properties may be used for generating synthetic (but realistic in a probabilistic view) equiprobable realizations of the physical conditions the EMI sensor will potentially experience during the monitoring, which may be used for addressing the inversion of EMI reading.

The main objective of this paper is to propose an approach to improve the parameters' optimization and the constraints in the time-lapse EMI inversion using soil water and solute modelling. In the paper we will show how the synthetic tests may be used to guide the optimization of inversion parameters and understanding of the impact of solute concentration and water content variations on EMI σ_a readings.

The key features of this study are (i) performing a controlled irrigation experiment, allowing for the simulation of the spatial and temporal variability of soil salinity during the irrigation experiment; (ii) monitoring of σ_a using the multi-coil CMD Mini-Explorer EMI sensor which takes σ_a measurements over six different depth sensitivity ranges; (iii) running numerical simulations using a physically based hydrological model to study how well the EMI survey and time-lapse inversion can resolve the σ_b distribution in space and time in our experimental set-up and to infer the best inversion parameters; (iv) inverting time-lapse field σ_a to map the spatio-temporal variability of σ_b and to interpret them in terms of concentration distributions over time.

2 Material and methods

2.1 Experimental set-up

The experiment was carried out at the Mediterranean Agronomic Institute of Bari in south-eastern Italy. The soil is classified as Colluvic Regosol, consisting of silty loam material of an average thickness of 0.7 m lying on fractured calcarenite bedrock. Figure 1 shows the experimental set-up, consisting of four experimental plots each of 30 m length and 3.6 m width, equipped with a drip irrigation system, with nine irrigation dripper lines placed 0.40 m apart, and characterized by an inter-dripper distance of 0.20 m. Pressure self-compensating drippers were used, and the emission uniformity was over 90 %. The soil was bare during the experiment period to avoid the effect of root uptake on the interpretation of the results. The four experimental plots were irrigated with water at four different salinity levels, experimental plot 1: water at 1 dS m⁻¹ (hereafter referred to as P1); experimental plot 2: water at 4 dS m⁻¹ (hereafter P2); experimental plot 3: water at 8 dS m⁻¹ (hereafter P3); experimental plot 4: water at 12 dS m⁻¹ (hereafter P4).

Figure 1 A schematic display of the experimental set-up. Four experimental plots were designed equally and irrigated with the same amount of water at four different salinity levels.

The soil–bedrock depth was measured at 40 points by augering, and the resulting spatial distribution is shown in Fig. 2. An apparent variability in the depth to bedrock was revealed, which may have significant impacts on the spatio-temporal distribution of water and solute concentration.

Figure 2 Measured spatial distribution of the soil depth along four experimental plots, P1, P2, P3 and P4.

Irrigations were started on 30 September 2016, and all experimental plots were irrigated 12 times until 21 November 2016 with the same amount of water and delivering day. Water salinity was induced by adding calcium chloride (CaCl₂) to well water with a salinity of about 1 dS m⁻¹. Each saline application was of about 18 mm, with a Cl⁻ concentration of about 0.1 mmol cm⁻³. The volumes of supplied water were calculated according to the differences between two consecutive level measurements in a Class A evaporimeter. The details of irrigation events and precipitation information are given in Fig. 3. EMI measurements were taken four times during the experiment period in each experimental plot along three transects, 1.2 m apart and at 2 m spacing, as shown in Fig. 1. All measurements were taken 1–3 d after the irrigations, allowing relatively time-stable water contents to be assumed at each site throughout the monitoring phases.

Figure 3 The details of irrigation events and precipitation information during the experiment. The dates of EMI measurements are marked with triangles.

2.2 EMI analysis

2.2.1 Characteristics of the EMI sensor

σ_a data were collected using a CMD Mini-Explorer device (GF Instruments, Brno, Czech Republic). The characteristics of this device make it especially suitable for monitoring electrical conductivity at relatively shallow depths. The CMD Mini-Explorer was used to measure σ_a in VCP (vertical coplanar, i.e. horizontal magnetic dipole configuration) mode and then HCP (horizontal coplanar, i.e. vertical magnetic dipole configurations) mode by rotating the probe 90^∘ axially to change the orientation from VCP to HCP mode. The probe has three receiver coils with 0.32 (ρ32), 0.71 (ρ71) and 1.18 (ρ118) m distances from the transmitter coil and operates at 30 kHz frequency. With the largest coil spacing, the instrument has an effective depth investigation of 1.8 m in the HCP mode and 0.9 m in the VCP mode.

2.2.2 Electromagnetic forward model

The electromagnetic forward modelling is solved by applying the full solution of the Maxwell equations to calculate the σ_a responses of a 1D model. The response of a layered media (secondary magnetic field) excited by a small horizontal transmitter loop above the ground surface can be expressed in terms of the Hankel transform (Zhang et al., 2000):

$$\begin{array}{}\text{(1)}& H_{\mathrm{z}}={\displaystyle \frac{M}{\mathrm{4}\mathit{\pi }}}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\left[R_{\text{TE}}{e}^{-\mathit{\beta }(\mathrm{2}{z}_{\mathrm{s}}-z_{\mathrm{r}})}\right]{\mathit{\beta }}^{\mathrm{2}}{J}_{\mathrm{o}}\left(\mathit{\beta }r\right)\mathrm{d}\mathit{\beta },\end{array}$$

where M is the moment of the transmitter loop emitting at an angular frequency ω, and z_s and z_r are the heights of the transmitter and receiver loops, respectively. r is the transmitter–receiver distance, and J_o is the Bessel function of first kind and order zero.

The kernel function R_TE is defined as

$$\begin{array}{}\text{(2)}& R_{\text{TE}}={\displaystyle \frac{Z^{\mathrm{1}}-Z_{\mathrm{0}}}{Z^{\mathrm{1}}+Z_{\mathrm{0}}}},\end{array}$$

where Z₀ is the intrinsic impedance of free space, and Z¹ is the input impedance at the first layer calculated by a recursive procedure. Similar equations can be written for vertical coplanar loops. In such cases the secondary magnetic field is

$$\begin{array}{}\text{(3)}& H_{\mathrm{y}}={\displaystyle \frac{M}{\mathrm{4}\mathit{\pi }r}}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\left[R_{\text{TE}}{e}^{-\mathit{\beta }(\mathrm{2}{z}_{\mathrm{s}}-z_{\mathrm{r}})}\right]\mathit{\beta }{J}_{\mathrm{1}}\left(\mathit{\beta }r\right)\mathrm{d}\mathit{\beta },\end{array}$$

where J₁ is the first-order Bessel function.

σ_a (mS m⁻¹) is usually calculated assuming the low induction number (LIN) approximation using the formula

$$\begin{array}{}\text{(4)}& {\mathit{\sigma }}_{\mathrm{a}}={\displaystyle \frac{\mathrm{4000}}{\mathit{\omega }{\mathit{\mu }}_{\mathrm{o}}{r}^{\mathrm{2}}}}{\left({\displaystyle \frac{H_{\mathrm{s}}}{H_{\mathrm{p}}}}\right)}_Q,\end{array}$$

where μ₀ is permeability of free space, and Q denotes the out-of-phase component of the secondary to primary magnetic field coupling ratio.

2.2.3 Time-lapse inversion

With the time-lapse inversion, one seeks to calculate the temporal variation of the conductivity along a transect. The quasi-2D inversion algorithms based on Monteiro Santos (2004) with a modification of the algorithm proposed by Kim et al. (2009) were used in this study. The problem is resolved in both algorithms iteratively starting from a uniform model. Two different levels of constraints, S1 and S2, were applied. In the S1 option, the corrections to the model parameters at each iteration are calculated by solving the system of equations

$$\begin{array}{}\text{(5)}& \left({\mathbf{J}}^{\mathrm{T}}\mathbf{J}+\mathit{\lambda }{\mathbf{C}}^{\mathrm{T}}\mathbf{C}+\mathit{\alpha }{\mathbf{M}}^{\mathrm{T}}\mathbf{M}\right)\mathit{\delta }\mathit{p}={\mathbf{J}}^{\mathrm{T}}\mathit{b}-\mathit{\alpha }{\mathbf{M}}^{\mathrm{T}}\mathbf{M}\mathit{p}.\end{array}$$

In the S2 option, the corrections of the parameters at each iteration are calculated solving the equations

$$\begin{array}{l}\text{(6)}& {\displaystyle }\left({\mathbf{J}}^{\mathrm{T}}\mathbf{J}+\mathit{\lambda }{\mathbf{C}}^{\mathrm{T}}\mathbf{C}+\mathit{\alpha }{\mathbf{M}}^{\mathrm{T}}\mathbf{M}\right)\mathit{\delta }\mathit{p}={\displaystyle }{\mathbf{J}}^{\mathrm{T}}\mathit{b}+\mathit{\lambda }{\mathbf{C}}^{\mathrm{T}}\mathbf{C}\left(\mathit{p}-{\mathit{p}}_{\mathrm{o}}\right)\mathit{\alpha }{\mathbf{M}}^{\mathrm{T}}\mathbf{M}\left(\mathit{p}-{\mathit{p}}_{\mathrm{o}}\right),\end{array}$$

where δp is the vector comprising corrections of the parameters (logarithm of conductivities, p_j) of an initial model; p_o refers to a reference model; b is the vector containing the differences between the logarithm of the observed and calculated apparent conductivities. J is the Jacobian matrix with elements given by $\frac{{\mathit{\sigma }}_j}{{\mathit{\sigma }}_{ai}^c}\frac{\partial {\mathit{\sigma }}_{ai}^c}{\partial {\mathit{\sigma }}_j}$. λ is a Lagrange multiplier and determines the amplitude of the parameter corrections in the space domain, and the regularization matrix C stores the coefficients of the spatial roughness of the model parameters at time t, which is defined as

$$\begin{array}{}\text{(7)}& \mathit{\delta }{\mathrm{r}}_{\mathrm{j}}={\mathit{\delta }}_{\text{PjE}}+{\mathit{\delta }}_{\text{PjW}}-\mathrm{4}{\mathit{\delta }}_{\mathrm{P}}+{\mathit{\delta }}_{\text{PjN}}+{\mathit{\delta }}_{\text{PjS}},\end{array}$$

where the elements of matrix C are 1 or −4 according to the position of the neighbours. α is a parameter that determines the amplitude of the parameter corrections in the time domain, and M is a square matrix that accounts for the temporal continuity of the model parameters.

The elements of matrix M are defined in terms of models at time t−dt and t+dt. The model misfit is calculated using the following equation:

$$\begin{array}{}\text{(8)}& \text{misfit}=\sqrt{{\displaystyle \frac{\mathrm{1}}{N}}}{\sum }_{i=\mathrm{1}}^N\text{ln}{\left({\mathit{\sigma }}_{\mathrm{a}i}^{\mathrm{o}}-{\mathit{\sigma }}_{\mathrm{a}i}^{\mathrm{c}}\right)}^{\mathrm{2}},\end{array}$$

where N is the number of apparent conductivity value, with ${\mathit{\sigma }}_{\mathrm{a}i}^{\mathrm{o}}$ and ${\mathit{\sigma }}_{\mathrm{a}i}^{\mathrm{c}}$ representing observed and calculated σ_a, respectively.

In this algorithm, two regularizations are imposed in both space and time domains. Consequently, the spatial and temporal Lagrangian multipliers have to be optimized. The spatial Lagrangian multiplier (λ) controls the relative importance of spatial model smooth and data-response misfit and decreases gradually during the inversion process to resolve more detailed model parameters. A larger λ tends to generally produce a model with a larger misfit error but smooth variation of conductivity values. A larger λ is usually acceptable if the soil conductivity changes in a smooth manner and allows a model to be produced that is reasonably more realistic. In contrast, a smaller λ is usually required when sharper soil conductivity changes are expected in order to resolve the sharp boundaries. A suitable λ value is usually determined empirically based on the expected distribution of σ_b and by performing inversions with different values. The second regularization, the temporal Lagrangian multiplier (α), is the temporal damping factor that gives the weight for minimizing the temporal changes in the conductivity along the time axis. α is a constant value and is defined by the similarity of the two consecutive reference times. The larger the α value, the more similar the reference models are that result from the inversion; a value of zero means no temporal constraints are applied (i.e. a traditional non-time-lapse inversion). The misfit function in this algorithm is the square root of the sum of the squares of the differences divided by the number of the measurements and is expressed in millisiemens per metre (mS m⁻¹).

2.3 Synthetic experiment

2.3.1 Hydraulic properties' dataset

A large dataset of hydraulic properties was already available from the experimental site, which was obtained by previous laboratory and field hydraulic characterizations carried out during several measurement campaigns (e.g. Coppola et al., 2011a, b, 2013, 2015). A total of 70 soil samples were collected from the Ap and Bw horizon in the experimental farm. Saturated hydraulic conductivity, K₀, and water retention experimental data were measured in the laboratory by the falling head permeameter (Reynolds and Elrick, 2003) and tension table method (Dane and Hopmans, 2003), respectively. The water retention data were fitted to the van Genuchten model (van Genuchten, 1980):

$$\begin{array}{}\text{(9)}& S_{\mathrm{e}}={\displaystyle \frac{\mathit{\theta }-{\mathit{\theta }}_{\mathrm{r}}}{{\mathit{\theta }}_{\mathrm{0}}-{\mathit{\theta }}_{\mathrm{r}}}}={\left[\mathrm{1}+{\left|\mathit{\alpha }h\right|}^n\right]}^{-m},h<\mathrm{0}\end{array}$$

$$\begin{array}{}\text{(10)}& \mathit{\theta }={\mathit{\theta }}_{\mathrm{0}},h=\mathrm{0},\end{array}$$

where S_e (–) is the effective saturation, α (cm⁻¹), n (–) and m (–) are shape parameters, and θ₀ (–) and θ_r (–) are the saturated and residual water content, respectively.

The hydraulic conductivity was estimated by the van Genuchten–Mualem model (van Genuchten, 1980), with $m=\mathrm{1}-n^{-\mathrm{1}}$:

$$\begin{array}{}\text{(11)}& K_{\mathrm{r}}\left(S_{\mathrm{e}}\right)={\displaystyle \frac{K\left(S_{\mathrm{e}}\right)}{K_{\mathrm{0}}}}=S_{\mathrm{e}}^{\mathit{\tau }}{\left[\mathrm{1}-{\left(\mathrm{1}-S_{\mathrm{e}}^{m^{-\mathrm{1}}}\right)}^m\right]}^{\mathrm{2}},\end{array}$$

where K₀ is the saturated hydraulic conductivity, K_r (–) is the relative hydraulic conductivity and τ (–) is a parameter which accounts for the dependence of the tortuosity and the correlation factors on the water content.

Due to the different hydrological behaviour between the laboratory and the field (Kutilek and Nielsen, 1994), the laboratory-derived curves were scaled to the field curves by applying the procedure described in Basile et al. (2003, 2006). Statistics of the parameters are reported in Table 1. The data in the table indicate a larger variability of the hydraulic properties of the Bw horizon, compared to the Ap horizon. This is probably due to the frequent roto-tillage of the Ap horizon, inducing significantly greater homogeneity of this soil layer. Actually, even the apparent larger variability of the K₀ for the Ap horizon is much lower than the variability reported in the literature, which is generally characterized by a coefficient of variation (CVs) much larger than 100 % (Kutilek and Nielsen, 1994; Mallants et al., 1996; Coppola et al., 2011a) and may even reach values of 450 % (Carsel and Parrish, 1988). The hydraulic properties' parameters for the bedrock were those determined on the same type of bedrock by Caputo et al. (2010, 2015).

Table 1 Statistics (mean and coefficient of variation, CV) of the hydraulic property parameters of the two soil layers and bedrock.

2.3.2 Synthetic hydrological simulations

2D numerical simulations of the infiltration and redistribution process were carried out using Hydrus 2D/3D software (Šimůnek et al., 2016), by introducing the actual boundary conditions imposed during the experiment to (i) gain an insight into the spatio-temporal distribution of water and solute concentration during the experiment period, prior to interpreting the field EMI data, and to (ii) optimize the EMI inversion parameters. The Richards equation and advection–dispersion equation were used to simulate water flow and solute transport, respectively. The distributions of the hydraulic parameters (see Table 1) were used to generate several synthetic transects, with variable hydraulic properties and depth to bedrock (electrically resistive layer), in order to simulate different distributions of water contents, solute concentrations and soil depths and to explore their role in the variability of the σ_a response. Each synthetic transect consisted of an assembly of 30 interacting soil columns, each with its own hydraulic properties and depth to bedrock. Each of the 30 soil profiles included three horizons: Ap (0–15 cm), Bw and bedrock. The Bw and bedrock thickness changed in each soil profile according to the depth of the soil–bedrock interface. For hydraulic properties, the transect variability was assumed to be characterized by the variability of the parameters θ₀, α, n and K₀. The hydraulic properties and depth to bedrock were assigned randomly to each soil column by using the statistical distribution of each parameter (see Table 1). Statistical tests showed that θ₀ and n were normally distributed, K₀ and α log-normally distributed and depth to bedrock uniformly distributed. The means and the covariance matrix for θ₀, log K₀, log α, and n were computed. Due to continuous ploughing and other tillage practices (see Sect. 2.1), the Ap horizon is rather homogeneous, and therefore only the hydraulic parameters of the horizon Bw were considered to vary stochastically. The hydraulic parameters of both Ap horizon and bedrock were fixed to their average values.

A total of 20 synthetic transects with 30 random vectors of the four parameters (θ_s, α, n and K₀) were produced from the correlated multivariate distribution by generating a vector x of independent standard normal deviates and then applying a linear transformation of the form $\mathit{x}=\mathit{m}+\mathbf{L}{\text{r}}_n$, where m is the desired vector of means, and L is the lower triangular matrix derived from the symmetric covariance matrix V=LL^T decomposed by Cholesky factorization (Carsel and Parrish, 1988). θ_r was set to zero for all simulations. For solute transport, a longitudinal dispersivity of 2 cm was assumed according to a previous experiment carried out in the same field (Coppola et al., 2011b). Transverse dispersivity was assumed to be one-tenth of the longitudinal dispersivity (Mallants et al., 2011). The simulation time – according to the field experiment data – included a period of 91 d from 1 September to 30 November 2016 in which 12 irrigations were applied (a total of 210 mm), the water with an electronic conductivity (EC) of 12 dS m⁻¹ (each one of about 18 mm with a Cl⁻ concentration of about 0.1 mmol cm⁻³). The level of 12 dS m⁻¹ was considered because spatial and temporal variabilities of σ_b were expected to be larger and more apparent due to greater salinity changes. Potential evapotranspiration was estimated by a local agrometeorological station; irrigation fluxes and Cl⁻ concentration in the irrigation water were considered to be top boundary conditions. Free drainage was assumed at the lower boundary ($z=-\mathrm{150}\mathrm{cm}$). The initial soil water content value was set to 0.25 cm³ cm⁻³, based on the field measurements by a TDR probe, and the initial Cl⁻ concentration was set equal to zero.

It is worth noting that the use of synthetic transects does not aim to address the overall spatial variability of soil properties potentially observable in the investigated field but to randomly select a reasonable number of different scenarios to better understand how solute concentration and water content changes during the experiment influence σ_b distribution. They also help to identify a proper regularization strategy to invert measured σ_a data in the field. In this sense, the number (20) of transects is a trade-off between the need of accounting for the possible heterogeneity of the field and the computational challenge of carrying out too many synthetic data inversions.

For each generated transect, numerical simulations produced distributions of water contents and Cl⁻ concentrations.

2.4 Site-specific calibration $\mathit{\theta }-{\mathit{\sigma }}_{\mathrm{w}}-{\mathit{\sigma }}_{\mathrm{b}}$

The water content and Cl⁻ concentration distributions were converted to bulk electrical conductivity (σ_b) distributions by using the model proposed by Malicki and Walczak, 1999:

$$\begin{array}{}\text{(12)}& {\mathit{\sigma }}_{\mathrm{w}}={\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{b}}-\mathrm{a}}{\left({\mathit{\epsilon }}_{\mathrm{b}}-\mathrm{b}\right)\left(\mathrm{0.0057}+\mathrm{0.000071}\mathrm{S}\right)}},\end{array}$$

where ε_b (–) is the dielectric constant, which is related to the water content, and σ_w (dS m⁻¹) is the electrical conductivity of the soil solution. The latter was obtained by using a linear relationship C–σ_w for solutions at different concentrations of calcium chloride.

The parameters for the Malicki and Walczak model were calibrated through a laboratory experiment. Specifically, ε_b and σ_b for different values of σ_w were simultaneously measured on reconstructed soil samples values using TDR probes. For this purpose, four PVC cylinders (8 cm in diameter and 15 cm height) were filled with air-dried soil, reaching a dry bulk density of about 1.1 g cm⁻³, imitating the field condition. Each soil sample was wetted by adding 10 mL of CaCl₂ solution at a specified electrical conductivity: 1, 2, 4 and 6 dS m⁻¹, respectively. The cylinders were covered by 0.05 mm plastic foil before the measurement in order to avoid evaporation and to equilibrate with the air temperature of 20 ^∘C. The procedure was repeated 16 times for each soil sample to measure soil water content values ranging from air-dried to near saturation. For each wetting step, the measurements of θ and σ_b were carried out using TDR three-wire probes (10 cm long with a rod diameter of 0.3 cm and rods spaced 1.2 cm) vertically inserted in the soil columns. The following parameters of Eq. (10) were obtained: a=3.6 dS m⁻¹; b=0.11.

2.5 Schematic view of the approach used in the paper

The logical sequence of the different steps used in the proposed approach is described in Fig. 4.

Figure 4 Flow chart of the applied procedure, with the key steps in bold. The procedure is explained in detail in the text.

1. Starting from the statistical distribution of the hydraulic properties and the physical characteristics of the system under study (the latter limited to only depth to the bedrock, in this specific case), as described in Sect. 2.3.1, several synthetic transects are generated, each accounting for some of the variability in terms of hydraulic properties and physical characteristics that the EMI sensor potentially experiences during the monitoring.

2. Hydrological simulations for each of these synthetic transects are thus carried out, each producing synthetic distributions of both water contents and solute concentrations, as described in Sect. 2.3.2, which, in turn, are converted to as many synthetic distributions of σ_b by using a specifically developed $\mathit{\theta }-{\mathit{\sigma }}_{\mathrm{b}}-{\mathit{\sigma }}_{\mathrm{w}}$ calibration relationship, discussed in Sect. 2.4.

3. These σ_b distributions are used in a forward EMI modelling procedure (see Sect. 2.2.2.) to generate synthetic σ_a data for all the synthetic transects considered.

4. The σ_a, θ and σ_w (C) distributions are used to guide the inversion of EMI readings obtained during real monitoring campaigns carried out in the physical system under study in two different (but related) ways: (i) to have a priori knowledge of where a measured EMI reading may come from (e.g. at a given time, a measured σ_a distribution could come from the depth of the water or solute propagation front, from the depth to the bedrock or from water accumulation at the soil–bedrock interface); (ii) to identify the optimal regularization parameters, discussed in Sect. 2.2.3, to be used in the inversion model of the real EMI data, by looking for example to the parameter value, allowing for satisfactory results for most of the synthetic transects. For the latter purpose, we compare different inversion parameter combinations in terms of correlation (correlation coefficient, R), precision (mean square error, RMSE) and bias (mean error, ME) between the simulated and modelled σ_b.

5. Finally, with all this information at hand, one is ready to look at the real EMI datasets, by producing more realistic σ_b distributions which can now be seen with much more confidence compared to cases where EMI inversion has to be carried out without any prior information from the experiment and, thus, producing less interpretable and more uncertain σ_b distributions.

The Results and discussion section below will show the application of the approach to the system under study, with an analysis of the real EMI data carried out only in the final phase, when all the information needed are available to guide the inversion of the real EMI data and interpretation of the obtained models.

3 Results and discussion

3.1 Synthetic spatio-temporal distributions of solute water content and solute concentration

This section shows how the synthetic simulations described in Sect. 2.3.2 can be used to analyse the sensitivity of the EMI response to both the hydrological behaviour and physical characteristics of the system under study.

As an example, Fig. 5 depicts the depth to the bedrock for one of the synthetic transects described in Sect. 2.3.2. For the same transect, Fig. 6 shows, for 4 selected days, the spatial distribution of water content obtained from the Hydrus 2D/3D simulations with irrigation water at 12 dS m⁻¹. The soil shows very high values of water content on the selected days. This is because the water was supplied to the soil only 1–3 d before the selected days. On the other hand, lower values were obtained at larger depths within the bedrock. The water content distributions show a significant lateral heterogeneity, while the temporal variations are very small. The spatial variations of water content distribution can be partly explained by the variability of the hydraulic properties. However, soil depth proved to be the dominant factor in determining water content lateral variability. In fact, high inverse correspondence of soil depth with water content is revealed in the vertical profiles where bedrock is superficial (e.g. profiles 3, 4 and 5) or deep (i.e., profiles 19, 20 and 21). The correlation between soil depth and water content averaged along each profile is very high (r=0.88) for the four dates, thus confirming that the depth of the bedrock is the main factor governing the water distribution.

Figure 5 Simulated spatial distribution of the soil depth for the selected scenario.

Figure 6 Soil water content distribution simulations for the selected scenario: (a) 17 October, (b) 26 October, (c) 14 November and (d) 23 November 2016.

In terms of the temporal variations, the water content distribution did not change significantly, and the average soil water content of the whole soil profile was found to be similar on the selected dates, with an average in the range 0.20–0.22 cm³ cm⁻³ and a standard deviation of 0.02 cm³ cm⁻³. This was expected as the experimental field was irrigated regularly, and the four selected dates refer to approximately the same time after an irrigation application (1–3 d). However, small differences of water content may be observed near the soil surface, which may well be explained by the evaporation process taking place during the 1–3 d after irrigation and mostly involving the shallower soil layer. For example, the near surface shows higher water content (i.e. 0.25–0.30 cm³ cm⁻³) on 26 October because the irrigation took place 1 d before on 25 October. On the other hand, lower water content is evident in the near surface on 17 October due to a 3 d gap between the water irrigation on 14 October and the simulation date. Figure 7 shows the spatio-temporal distribution of Cl⁻ concentration obtained for the same synthetic transect of Fig. 6. Compared to the water content, Cl⁻ concentrations show a significantly greater evolution over time, with a slow and steady Cl⁻ concentration increase along the soil profile due to the 12 injections of saline water. On average, the front of chloride deepens slowly at a fairly constant rate but with a lateral variability, which is related to the lateral variability of water contents and hydraulic properties. The lateral variability is quite low on the first days of solute applications and becomes evident only when the solute migrates deeper into the soil profile. The lateral variability of Cl⁻ concentration is mainly related to the depth to bedrock, as can be immediately observed by comparing the solute distributions and the transect morphology shown in Fig. 5. The depth of the soil–bedrock interface as well as soil hydraulic properties condition the spatial distributions of water contents, which, in turn, influence the Cl⁻ concentration distribution.

Figure 7 Soil Cl⁻ concentration distribution simulations for the selected scenario: (a) 17 October, (b) 26 October, (c) 14 November and (d) 23 November 2016.

3.2 Simulated spatio-temporal distribution of σ_b

Figure 8 shows the spatio-temporal distribution of σ_b for the same synthetic transect, obtained by converting water contents and solute concentrations by applying Eq. (10). Figure 8 shows a resistive zone beneath a conductive zone. The conductivity of the resistive zone varies slightly spatially and temporally; however the conductivity of this zone is generally within 25 mS m⁻¹. The resistive zones in the maps correspond to the bedrock in the study area (see Fig. 5). The conductivity of the upper layers changes significantly both spatially and temporally from an average conductivity of 50 mS m⁻¹ on 17 October to more than 100 mS m⁻¹ on 23 November. The time evolution of the conductive zone is evident and, as the water content does not change significantly over time, is mostly related to the chloride propagation during the simulation. The lateral variation at any time, by contrast, is largely ascribable to bedrock topography.

Figure 8 σ_b distribution simulations for the selected scenario: (a) 17 October, (b) 26 October, (c) 14 November and (d) 23 November 2016.

3.3 Time-lapse synthetic σ_a data

The spatio-temporal distribution of σ_b shown in Fig. 8 was used to generate the time-lapse synthetic σ_a data. The generated σ_a data for the model obtained on 23 November (see plot d in Fig. 8) are shown in Fig. 9. Profile ρ32 shows the greatest σ_a values in each orientation, while the minimum σ_a was recorded on profile ρ118, indicating a conductive zone over a resistive zone, which is expected from the model shown in Fig. 8d. In addition, significant lateral σ_a change is evident along the transect, with strongest fluctuations in ρ32 in both orientations. This is because of the strong lateral σ_b variations at the near surface due to the saline water infiltration and heterogeneity of the subsurface. The general behaviour of the synthetic σ_a data suggests that the field data might have strong lateral σ_a changes along the transect, and a careful processing of data is required (e.g. filtering of σ_a data should be avoided).

Figure 9 σ_a data generated from the forward calculation of the σ_b distributions on 23 November, shown in Fig. 8d. (a) HCP-mode and (b) VCP-mode configurations are displayed.

3.4 Optimizations of the inversion parameters

In this section we show how well the developed inversion method can resolve the σ_b distribution from generated synthetic σ_a distribution. Firstly, we investigated the influence of the spatial smoothing parameter (λ) and the inversion algorithm S1 (Eq. 5) and S2 (Eq. 6) in the Supplement in resolving the σ_b distribution. In this regard, the generated σ_a data were inverted using both S1 and S2 and different values of λ in the 0.01 to 10 range. Figure 10 shows an example of the σ_b distribution models after inverting the synthetic σ_a data presented in Fig. 9 This example was selected due to the larger lateral and vertical contrast.

Figure 10 σ_b distributions obtained from the inversion of the synthetic σ_a data (23 November), shown in Fig. 9, using four different combinations of two algorithms, S1 and S2, and λ parameters.

Ideally, the obtained σ_b distribution should be very similar to the one shown in Fig. 8d. However, looking closely at the obtained σ_b distribution, we observe that all of them are different, to some extent, to the one shown in Fig. 8d. First of all, we note that the model obtained using λ values greater than 1 (not shown here) and regardless of the choice of inversion algorithm (i.e. S1 and S2) shows a high misfit error and also significantly over-smoothed σ_b distribution. This is not surprising as the sharp vertical spatial variability of σ_b due to the saline water irrigation and soil heterogeneity as well as the shallow bedrock cannot be well resolved using high λ (over-smoothed parameters). Consequently, a high λ is not a wise choice when sharp vertical and lateral conductivity contrasts are expected in the field. The obtained model using S2 algorithm and moderate λ values (i.e. 0.5) also yields a very smooth model. In contrast, the obtained models using S1–λ=0.05, S1–λ=0.5 and S2–λ=0.05 do a better job in resolving near-surface anomalies. However, the conductivity distribution at depth is not well recovered using the S1 algorithm.

In Fig. 11, we plotted the synthetic σ_b distribution data against the obtained modelled σ_b distributions using different inversion parameters shown in Fig. 10, and we calculated the statistical scores R, RMSE and ME to further investigate the impact of the inversion algorithm and parameters in resolving the σ_b distributions. Firstly, we observed that both S1 and S2 algorithms underestimated the σ_b to some extent, judging from the ME values. The S2–λ=0.5 and S1–λ=0.05 show weaker statistical scores, with higher RMSE and ME and relatively lower R. The over-smoothed impact of S2–λ=0.5 and the resulting highest ME indicates that the S2 algorithm with moderate to high λ values is not a rational choice when large spatial and vertical variations of σ_b are expected. On the other hand, the lowest R and highest RMSE, obtained for S1–λ=0.05, suggest that the S1 algorithm with very small values of λ is not able to predict the spatial variation of σ_b well. The use of the S1 algorithm with very small values of λ usually applies an insignificant spatial constraint that may result in more inconsistency between synthetic σ_b and modelled σ_b distributions. The S2–λ=0.05 and S1–λ=0.5 show better statistical results, with higher R and lower RMSE and ME, which makes them a better choice for the inversion of field data. While both inversion parameters set present almost the same RMSE, the S2–λ=0.05 yields a higher R, suggesting that the S2–λ=0.05 can better resolve the spatial σ_b distributions. This is expected from a comparison of S2–λ=0.05 and S1–λ=0.5 results in Figs. 10 and 11, where a relatively lower correlation is evident between synthetic σ_b and modelled σ_b at lower ranges of σ_b (Fig. 11) located at depths of more than 70 cm (Figs. 10 and 11) when we used the S1–$\mathit{\lambda }=\mathrm{0.5}.$ The difficulty of resolving a resistive zone at depth and beneath a conductive zone is indeed expected. In fact, the sensitivity of the EMI signals is very limited over the resistive zone, and therefore the resistive zone cannot be well resolved. The condition will be worse in our study as a resistive zone is located beneath a conductive zone: the EMI response of the subsurface will be dominated by the influence of the near-surface conductive zone. In addition, five of the six depths of investigation of the CMD Mini-Explorer are limited to the first 1 m, and, as a result, a lower resolution is expected at greater depths. The S2 algorithm did a better job in resolving the resistive zone at depth. This is because the S2 algorithm constrains the value of each cell to the reference model during the inversion process, which limits the large variations of each cell during the inversion process.

Figure 11 Synthetic versus modelled σ_b distribution for four different combinations of S and λ inversion parameters. Different colours refer to different depths. Main statistics are also reported.

In terms of recovering the absolute values of soil electrical conductivity, it appears that all obtained models underestimate the conductivity of the anomalies near the surface. This can be explained by the over-parameterized inverse problem and the effects of smoothing from regularization applied in the inversion algorithm, as well as the impact of the resistive zone beneath the conductive zone. In addition, the six measurements per site with a site spacing of 1 m are not sufficient for recovering sharp σ_b variability along the transect. Measuring σ_a at different heights as well as smaller site spacing enables more σ_a data to better resolve changes which occur over short depth and length increments.

In terms of the influence of the temporal smoothing parameter (α), we explored different values (results not shown here). We noticed that the values larger than 0.1 over-smoothed the expected temporal variation, and thus the detailed variations cannot be resolved. We conclude that, for our study, a value of 0.05 is the best choice for α in resolving the temporal variation of σ_b.

We repeated the same analyses using 20 different synthetics transects, consisting of a random assembly of 30 interacting soil columns, each with its own hydraulic properties and depth to bedrock and different saline treatment (results not shown here). The results of our analysis show that the algorithm S2 is the best choice in the presence of a resistive zone at depth and small values of λ and α work best for dealing with the sharp lateral and temporal expected changes that occurred during the experiment. Based on the synthetic tests, the spatio-temporal algorithm S2 described in Eq. (6) with λ and α values of 0.05 and 0.05 respectively were selected to invert time lapse actual σ_a datasets measuring during the experiment period.

3.5 Inversion of the real time-lapse σ_a field data

Using the optimized inversion parameters obtained in Sect. 3.4., we inverted time-lapse σ_a data collected over the four experimental plots, P1, P2, P3 and P4. Figure 12 shows the σ_a data of the middle transect in each experimental plot (Fig. 1) for the last date of monitoring, i.e. 23 November. The σ_a data show a relatively similar pattern in both VCP and HCP modes, with greater σ_a values on ρ32 and ρ71 and the minimum σ_a values recorded on ρ118, indicating a conductive zone over a resistive zone. In addition, greater lateral σ_a changes in the VCP mode are evident along the four transects, with noticeable fluctuations in P4, suggesting greater lateral σ_b variations at the near surface. The σ_a data obtained from plot P1 and P2 show a lower range of conductivity, varying in the 20–40 mS m⁻¹ and 20–50 mS m⁻¹ ranges, respectively. P3 and P4 represent a higher range of conductivity, with the former in the 20–70 mS m⁻¹ range and the latter in the 20–80 mS m⁻¹ range.

Figure 12 σ_a distribution of the middle transect in four experimental plots, P1, P2, P3 and P4, for 23 November. Panels (a), (b), (e) and (f) show the VCP configuration, and panels (c), (d), (g) and (h) show the HCP configuration.

Figure 13 shows the time-lapse inversion results for the same transects for four dates: 17 and 26 October and 14 and 23 November 2016. The same colour scale was used for all models, allowing for a comparison of spatio-temporal variation of σ_b along all profiles at different times. The corresponding model responses for the last measurements on 23 November are shown in Fig. 12. The misfit errors for the P1, P2, P3 and P4 are 1.43, 1.21, 2.15 and 2.95 mS m⁻¹, respectively, indicating a good fit between data and model responses. Slightly higher misfit errors for P3 and P4 are probably due to greater range of σ_a, as well as the larger lateral variations of σ_a. Looking at EMI models at different times and along all four experimental plots, we identify two distinct zones: a resistive zone at a depth of more than 50 cm beneath a conductive zone. The conductivity of the resistive zone varies laterally and temporally along all four transects; however, the conductivity of this zone rarely exceeds 25 mS m⁻¹. The resistive zones in the maps shown correspond to the bedrock in the study area (Fig. 2). The conductivity of this zone is slightly lower in P3 and P4, which is probably due to the shallower bedrock in these plots. In contrast, the conductive zone at the near surface shows significant spatio-temporal conductivity changes depending on the conductivity of the injected water.

Figure 13 Spatio-temporal distribution of σ_b along the middle transect in four experimental plots, P1, P2, P3 and P4, and for four selected dates: 17 October, 26 October, 14 November and 23 November 2016.

The time-lapse models obtained from plot P1 show the minimum spatio-temporal conductivity changes in this zone, with conductivity varying in the 30–60 mS m⁻¹ range. The P2 plot shows larger spatial and temporal conductivity changes compared to P1, as expected, with conductivity varying between 30 and 80 mS m⁻¹. The conductivity of this zone decreases slightly on 14 November. This is probably due to the impact of heavy rain (31 mm) on 7 November, inducing salt leaching and reducing the soil conductivity. The P3 and P4 plots show stronger conductivity changes both spatially and temporally, with conductivity varying between 50 and 100 mS m⁻¹. The first models, obtained for 17 October, shows the minimum conductivity among the four selected dates. This is consistent with the simulated distributions shown in Fig. 8, and it is probably due to the saline front being still undispersed in the initial propagation phase (Fig. 7).

As the saline water is continuously added to the soil surface, the Cl⁻ concentration increases in the soil and also deepens slowly in the soil profile, thus increasing soil electrical conductivity. Consequently, the higher soil conductivities seen on 23 October (Fig. 13, P3 and P4) are due to higher Cl⁻ concentration and its propagation to deeper layers. The conductivity models obtained for 14 November show a decrease in soil conductivity in both P3 and P4, although there were six irrigations after 23 October. On the other hand, the conductive zone extends to deeper soil. This is not surprising as we already expected it from the simulations (see Fig. 8). This behaviour is probably due to the rainfall event on 7 November (as noted on P2) which diluted the Cl⁻ concentration in the soil. This suggests that the EMI surveys and data modelling detected the expected change in Cl⁻ concentration well.

Finally, the models obtained for the data collected on 23 November, depicted in Fig. 13, show the maximum soil conductivity and extension of the conductive zone among the selected dates in P3 and P4. A comparison of these models with our simulations results discussed in Sects. 4.1 and 4.2 shows that the increase of Cl⁻ concentration in soil after 12 irrigation events, as well as the redistribution of Cl⁻ during the experiment period, is the main reason for increase in soil conductivity. Judging from the water content distribution maps, obtained from numerical simulations and shown in Fig. 6, we expect very small temporal variations of water content in all experimental plots. Consequently, the temporal variations of Cl⁻ concentration and distribution are expected to be the key factor in temporal variations of soil electrical conductivity.

From a management perspective, the discussion above suggests that using the described approach regularly after irrigation applications may allow for monitoring of salt propagation and redistribution. As water content values and patterns are expected to be quite similar at similar times after each irrigation event, changes in bulk electrical conductivities obtained by an EMI sensor will be closely related to changes in salt concentrations.

4 Conclusions

In this study, we carried out a time-lapse EMI survey over four experimental plots irrigated with water at four different salinity levels during 3 months. We examined how well the time-lapse EMI measurements and a time-lapse inversion algorithm can be used to monitor soil salinity variability in space and time through performing simulation experiments and inversion processes. Based on our detailed simulations and synthetic tests as well as the interpretation of the time-lapse models, the following main conclusions can be drawn:

1. The numerical simulations performed in this study allow for predictions of the spatial and temporal variability of soil salinity and water content during the irrigation experiment prior to the modelling of field σ_a data. This improves our understanding of how soil salinity and water content changes during the experiment influence σ_b distribution in time, which may be crucial for interpreting EMI models in terms of soil salinity and water content distribution. They also provided a synthetic time-lapse EMI dataset very similar to the field conditions, allowing for the optimization of data processing and allowing us to find the best inversion approach for the field experiment. From the synthetic analysis, we also established that the EMI measurements do not return enough subsurface information to resolve the expected sharp conductivity changes in this specific experiment.

2. A comprehensive investigation of our results and a joint interpretation of numerical simulations and time-lapse models reveal that the soil Cl⁻ concentration change is the key factor responsible for σ_b changes when the EMI surveys are repeated after the irrigation at the same time. In fact, repetition of the EMI surveys after the irrigation has three main advantages: (i) the soil is wet and conductive, and, in such a condition, the signal-to-noise ratio is usually better and EMI data will be more reliable; (ii) the water content distribution will not change significantly, allowing the impact of soil salinity changes in time to be studied better; (iii) the sensitivity of σ_b to the Cl⁻ concentration in wet soils is higher, and thus EMI inversion results may be used to interpret salt propagation with more confidence.

3. A previously developed least-squares 4D space–time domain inversion algorithm was implemented in this study to invert entire time-lapse EMI datasets simultaneously for monitoring soil salinity for the first time. The regularizations in this algorithm are introduced in both space and time domains to improve the stability of the inversion process and to reduce the inversion artefacts. Using a synthetic test, the applicability of the algorithm has been examined, and we showed that this approach can provide information about the conductivity variability in space and time. More synthetic tests are required to further investigate the efficiency of the inversion algorithm under different scenarios; however, we anticipate that the algorithm can be used for other EMI monitoring surveys.

4. Performing EMI surveys over four experimental plots irrigated with water at four different salinity levels, we evaluated the potential of EMI surveys for monitoring the dynamics of soil salinity due to irrigation. Comparing the EMI models obtained from four experimental sites, we showed that the σ_b variations are consistent with the expectations related to the amount of salt and water irrigated at each plot during the experiment. The water content did not change significantly during the EMI measurement campaigns; hence, the temporal variations of σ_b as well as their difference at each plot are mainly related to the soil salinity distributions. In other words, the results indicate that the applied experimental methodology is strongly capable of giving information on the spatial and temporal variability of soil salinity. Further investigations have to be conducted to use EMI sensors for monitoring salinity under significant water content changes over time. In this case, discriminating the role of water content changes and salinity changes in the σ_a response may be quite complicated if not impossible.

5. The EMI method provides enormous advantages over traditional methods of soil sampling because it allows for in-depth and non-invasive analysis, covering large areas in less time and at a lower cost. However, a proper interpretation of the EMI inversion models in terms of soil process is usually difficult, owing to the fact that the soil electrical conductivity is a complex function of soil properties, which may vary significantly over space and time. Thus, retrieving soil properties from EMI data requires appropriate understanding of site-specific soil processes. Our study shows that the independent interpretation of time-lapse EMI data without hydrologic insight and understanding of soil processes may be misleading. This fact highlights the necessity of collaboration of geophysicists, soil scientists and hydrologists to construct a hydrologic conceptual model which can explain the salinity and water process.