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 semiarid 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 15 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 two-month period, with water at four different salinity levels (1, 4, 8 and 12 dSm-1) using a drip irrigation system. Time-lapse apparent electrical conductivity (σa) data were collected 6 times during the experiment period using a 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 20 and redistribution process. The obtained potential spatio-temporal distributions of water content, solute concentration and bulk electrical conductivity (σb) assisted in understanding of how solute concentration and water content changes during the experiment influence σb distribution as well as optimizing the time-lapse inversion parameters for resolving σb changes in this experiment. Finally, we inverted the time-lapse field σa data and interpreted the results in terms of concentration distributions over time. Our investigation shows that EMI measurements and suitable modelling techniques allow for rapid and non-invasive 25 investigation of spatio-temporal variability in soil salinity over large areas. Their effectiveness and relatively low cost make them appealing for management of water irrigation in salinity-threatened regions of the world.

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 million hectares in the European Union, mainly in the Mediterranean countries (Toth 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. However, it is very difficult to assess soil salinity on a management scale using traditional 35 methods. Soil salinity is traditionally assessed by measuring the electrical conductivity of a saturated soil paste (ECe) 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 40 volume of sensors is limited to tens of centimetres, thus extension of the information to a large area can be problematic (Coppola et al., 2015;Farzamian et al., 2015a).
The electromagnetic induction (EMI) technique is widely used as an alternative to traditional techniques for soil salinity assessment. It allows 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, 45 and cation exchange capacity; however, in a saline soil, the salinity is generally the dominant factor responsible for the spatiotemporal 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-50 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 55 process. Several inversion methods have been proposed to obtain the σb from the measured σa data, including the gradientbased inversion technique (Monteiro Santos, 2004;von Hebel et al., 2014;Schamper et al., 2012;Farzamian et al., 2019) 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., 2020a) and 3D (e.g. Huang et al., 2017). However, the potential of this method in 60 assessing temporal variability of soil salinity has not been fully explored. Several time-lapse inversion methods for directcurrent resistivity methods have been developed. These include the ratio method (Daily et al., 1992), the difference inversion (LaBrecque and Yang, 2001) and more recently, 4-D 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 dynamic and only few studies have been conducted to estimate soil water content changes (Huang et al., 2016;Whalley et al., 2017). Besides, a prerequisite for such an approach concerns the reliability of the inversion of the EMI result. In fact, inverting EMI data to obtain the 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-70 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. Farzamian et al., 2019), 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, numerical simulations of the same hydrological processes to be monitored by an EMI sensor 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 75 EMI inversion parameterization. In addition, such numerical simulations assist in better understanding of the complexity of the spatial and temporal variability of σb due to simultaneous variations of solute concentration and water content in irrigated lands.
In this study, we performed an irrigation experiment with the objective of evaluating the potential of the EMI method and a time-lapse inversion algorithm for monitoring spatial and temporal variability of soil salinity. The key features of this study 80 are: i) performing a very controlled irrigation experiment in an experimental site, allowing the simulation of the spatial and temporal variability of soil salinity during the irrigation experiment; ii) monitoring of σa using a multi-coil CMD Mini-Explorer EMI sensor which takes σa measurements over 6 different depth sensitivity ranges; iii) running numerical simulations by 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 85 the spatio-temporal variability of σb and to interpret them in terms of concentration distributions over time.

Background information
The experiment was carried out at the Mediterranean Agronomic Institute of Bari (MAIB) 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 90 bedrock.
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;. 70 soil samples were collected from the Ap and Bw horizon in the experimental farm of the MAIB. Saturated hydraulic conductivity, K0, and water retention experimental data were measured in the laboratory by the falling head 95 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): https://doi.org/10.5194/hess-2020-514 Preprint. Discussion started: 16 October 2020 c Author(s) 2020. CC BY 4.0 License.
The hydraulic conductivity was estimated by the van Genuchten-Mualem model (van Genuchten, 1980), with m=1-1/n: where K0 is the saturated hydraulic conductivity, Kr (-) 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. 105 Due to the different hydrological behaviour between laboratory and field (Kutilek and Nielsen, 1994), the laboratory-derived unsaturated 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 Ks for the Ap horizon 110 is quite lower than the variability reported in the literature, which is generally characterized by coefficients of variation (CV) much larger than 100% (Kutilek and Nielsen, 1994;Mallants et al., 1996;Coppola et al., 2011) 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;. Figure 1a 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 (EU) 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 dSm -1 (hereafter referred to as P1); 120 experimental plot 2: water at 4 dSm -1 (hereafter P2); experimental plot 3: water at 8 dSm -1 (hereafter P3); experimental plot 4: water at 12 dSm -1 (hereafter P4).

Experimental set-up 115
The soil-bedrock depth was measured in 40 points by augering and the resulting spatial distribution is shown in Fig. 1b. The figure shows an apparent variability in the depth to bedrock, which may have significant impacts on spatio-temporal distribution of water and solute concentration. 125 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 (CaCl2) to well water having a salinity of about 1 dSm -1 . The volumes of supplied water were calculated according to the differences between two https://doi.org/10.5194/hess-2020-514 Preprint. Discussion started: 16 October 2020 c Author(s) 2020. CC BY 4.0 License. consecutive measurements acquired by a Class A evaporimeter. The details of irrigation events and precipitation information are given in Fig. 2. 130 σ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 135 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. EMI measurements were taken 6 times during the experiment period (see Fig. 2) in each experimental plot along three transects, 1.2 m apart and at 2 m spacing, as shown in Fig. 1a. All measurements were taken 1-3 days after the irrigations, allowing relatively time-stable water contents to be assumed at each site throughout the monitoring phases. 140

Synthetic hydrological simulations
Preliminary 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 to the spatio-temporal distribution of water and solute concentration during the experiment period, prior to interpreting 145 the field EMI data and to ii) optimize the EMI inversion parameters. The Richards equation (RE) and Advection-Dispersion Equation (ADE) were used to simulate water flow and solute transport, respectively. A large data set of hydraulic properties is known for the farm where the experimental field is embodied as discussed in section 2.1. 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 150 soil depths and to explore their role on the variability of the σa response. Each synthetic transects 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 θ0, α, n and K0. The hydraulic properties and depth to bedrock were assigned randomly to each 155 soil column by using the statistical distribution of each parameter (see Table 1). Statistical tests showed that θ0 and n were normally distributed, K0 and α log-normally distributed, and depth to bedrock uniformly distributed. The means and the covariance matrix for θ0, log K0, log α, and n were computed. Due to the continuous ploughing and other tillage practices (see section 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. 160 Thus, for each saline treatment, 20 synthetic transects with 30 random vectors of the four parameters (θs, α, n and K0) 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 x=m+Lrn, 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 165 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 included a period of 91 days from 01 September 2016 to 30 November 2016 in which 12 irrigations were applied (a total of 210 mm). Four different salinity levels (i.e. 1, 4, 8 and 12 dSm -1 ) were applied in the simulations imitating the field experiment condition. Potential evapotranspiration was estimated by a local agrometeorological station; irrigation fluxes and Clconcentration in the irrigation water were 170 considered as top boundary conditions. Free drainage was assumed at the lower boundary (z=-150 cm). Initial soil water content value was set to 0.25 cm 3 cm -3 , based on the field measurements by TDR probe, and the initial Clconcentration 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 175 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 only 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 Clconcentrations. These, in 180 turn, were converted to σb distributions by using a specifically developed θ-σw-σb calibration relationship obtained directly in situ according to the procedure described in section 3.2. Finally, the simulated (synthetic) σb distributions were converted to σa distributions by using an electromagnetic forward model (see section 3.3). These σa values, in turn, were inverted again by using two time-lapse inversion algorithms, S1 and S2, described in section 3.4 to obtain distributions of σb to be compared to those obtained from hydrological simulations to identify the best regularization parameters for the inversion of the actual σa 185 data measured during the experiment.

Site-specific calibration θ-σw-σb
The water content and Clconcentration distributions were converted to bulk electrical conductivity (σb) distributions by using the model proposed by Malicki and Walczak, 1999: where εb (-) is the dielectric constant, which is related to the water content, and σw (dSm -1 ) 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 195 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 -3 , imitating the field condition. Each soil sample was wetted by adding 10 ml of CaCl2 solution at a specified electrical conductivity: 1, 2, 4 and 6 dSm -1 , 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-dry to near saturation. For each wetting 200 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. 3 were obtained: a = 3.6 dSm −1 ; b = 0.11.

Electromagnetic Forward Model
The electromagnetic forward modelling is solved by applying the full solution of the Maxwell equations to calculate the σa 205 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 Hankel transform (Zhang et al., 2000): where M is the moment of the transmitter loop emitting at an angular frequency , zs and zr are the heights of the transmitter and receiver loops, respectively. is the transmitter-receiver distance and a is the Bessel function of first kind and order zero. 210 The kernel function YZ is defined as: where 2 is the intrinsic impedance of free space, O is the input impedance at the first layer calculated by a recursive procedure. Similar equations can be written for vertical coplanar loops. In such case the secondary magnetic field is: where O is the first-order Bessel function.
σa (mSm -1 ) is usually calculated assuming the low induction number (LIN) approximation using the formula: where µ0 is permeability of free space and Q denotes the out-phase component of the secondary to primary magnetic field coupling ratio. 220

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 225 calculated by solving the system of equations: In the S2 option, the corrections of the parameters at each iteration are calculated solving the equations: where is the vector comprising corrections of the parameters (logarithm of conductivities, pj) of an initial model; po refers 230 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 . is a Lagrange multiplier and determines the amplitude of the parameter corrections in the space domain and the regularisation matrix C stores the coefficients of the spatial roughness of the model parameters at time t which is defined as: where N is the number of apparent conductivity value, with H a and H OE a 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 (λ and α respectively) have to be optimized. The spatial Lagrangian multiplier in this algorithm is not constant but is a spatially dependent variable and decreases gradually during the inversion process to resolve more detailed 245 model parameters. Generally larger values will produce smoother inversion results. A suitable damping factor value is usually determined empirically based on the expected distribution of σb and by performing inversions with different values. The second regularization is the minimization of the model roughness along the time axis. In contrast, the temporal Lagrangian multiplier is a constant value and is defined by the similarity of the two consecutive reference times. The larger the α value, the more similar are the reference models that result from the inversion; a value of zero means no temporal constraints are applied (i.e. 250 a traditional non-time-lapse inversion).

Results and Discussion
In terms of simulations of water flow and solute transport, one example of the scenarios simulated by water at 12 dSm -1 is presented here, where spatial and temporal variabilities of σb are expected to be larger and more apparent due to greater salinity changes. In addition, results of the field data for the middle transects of each experimental plot (i.e. P1, P2, P3, and P4) will 255 be discussed here. All results refer to four selected dates: 17 and 26 Oct, 14 and 23 Nov of 2016.

Simulated spatio-temporal distributions of solute water content and concentration
Figure 3 depicts depth to the bedrock for one of the selected scenarios, described in section 3.1, and simulated with water irrigation at 12 dSm -1 . Figure 4 shows the spatio-temporal distribution of water content obtained from the Hydrus 2D/3D simulations for the same transect. The soil shows very high values of water content on the selected days. This is because the 260 water was supplied to the soil only 1-3 days 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 265 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 balance.
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 in the shown dates, with an average in the range 0.20-0.22 cm 3 cm -3 270 and a standard deviation of 0.02 cm 3 cm -3 . This is expected as our experimental field was irrigated regularly, furthermore, note that the four selected dates refer to approximately the same time after irrigation. Nevertheless, the differences in simulation time and irrigation time differences (1-3 days) may be responsible for the small differences of water content observed near the soil surface where the evaporation process takes place and even one day of difference can have an effect. For example, the near-surface shows higher water content (i.e. 0.25-0.30 cm 3 cm -3 ) on 26 Oct, because the irrigation took place one day before 275 on 25 Oct. On the other hand, lower water content is evident in the near-surface on 17 Oct due to a three days gap between the water irrigation on 14 Oct and simulation date. Figure 5 shows the spatio-temporal distribution of Clconcentration obtained for the same synthetic transect shown above.
Compared to the water content, Clconcentrations show a significantly greater evolution over time, with a slow and steady Clconcentration increasing along the soil profile due to the eight injections of saline water from 17 Oct to 23 Nov (each one of 280 about 18 mm with a Clconcentration of about 0.1 mmol cm -3 ). 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.
We recall that simulations were all carried out by considering a dispersivity of 2 cm in all the 30 profiles making up the https://doi.org/10.5194/hess-2020-514 Preprint. Discussion started: 16 October 2020 c Author(s) 2020. CC BY 4.0 License.
transect. The lateral variability is quite low in the first days of solute applications and becomes evident only when the solute migrates deeper into the soil profile. The lateral variability of Clconcentration is mainly related to the depth to bedrock, as 285 can be immediately observed by comparing the solute distributions and the transect morphology shown in Fig. 3. The depth of the soil-bedrock interface as well as soil hydraulic properties conditions the spatial distributions of water contents which influence the Clconcentration distribution. Figure 6 shows the spatio-temporal distribution of σb for the same synthetic transect, obtained by converting water contents 290 and solute concentrations (Eq. 3). The figure 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 mSm -1 . The resistive zones in the maps correspond to the bedrock in the study area (see Fig. 3). The conductivity of the upper layers changes significantly both spatially and temporally from an average conductivity of 50 mSm -1 on 17 Oct to more than 100 mSm -1 on 23 Nov. The time evolution of the conductive zone is evident and, as the water content does not change significantly over time, 295 is mostly related to the chloride propagation during the simulation. The lateral variation at any time, by contrast, is largely ascribable to bedrock topography.

Time-lapse synthetic σa data
The spatio-temporal distribution of σb, shown in Fig. 6, was used to generate the time-lapse synthetic σa data. The corresponding generated σa data for the model obtained on 23 Nov (see Fig. 6d) is shown in Fig. 7. Profile ρ32 shows the 300 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. 6d. 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 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 305 (e.g. filtering of σa data should be avoided).

Optimizations of the inversion parameters
We evaluated how well the developed inversion method can resolve the σb distribution from generated σa data. Firstly, we investigated the influence of the spatial smoothing parameter (λ) and the inversion algorithm S1 (Eq. 8) and S2 (Eq. 9) in resolving the σb distribution. In this regard, the generated σa data were inverted using both S1 and S2 and different values of λ 310 in the 0.01 to 10 range. Figure 8 shows an example of the σb distribution model after inverting the synthetic σa data presented in Fig. 7. This example was selected due to the larger lateral and vertical contrast.
Ideally, the obtained σb distribution should be very similar to the one shown in Fig. 6d. 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. 6d. 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 315 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, 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 (Fig. 7a). In contrast, the obtained models using S1 -λ=0.05, S1 -λ=0.5 and S2 -λ=0.05 do 320 a better job in resolving near-surface anomalies. However, the conductivity distribution at depth is not well recovered using the S1 algorithm. 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 325 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.
In terms of recovering the absolute values of soil electrical conductivity, it appears that all obtained models underestimate the 330 conductivity of the anomalies near the surface. This can be explained by 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 6 measurements per site with site spacing of 1 m is 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. 335 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). 340 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. 9 with λ and α values of 0.05 and 0.05 respectively were selected to invert time lapse actual σa datasets measuring during the experiment period.

Inversion of field time-lapse EMI data 345
Using the optimized inversion parameters obtained in section 4.4, we inverted time-lapse σa data collected over the four experimental plots, P1, P2, P3, and P4. Figures 9 show the time-lapse inversion results for the middle transect in each experimental plot (Fig. 1) for four selected dates: 17 and 26 Oct; and 14 and 23 Nov 2016. The same colour scale was used for all models, allowing comparison of spatio-temporal variation of σb along all profiles at different times.
Looking at EMI models at different times and along all four experimental plots, we identify two distinct zones: a resistive zone 350 at depth 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 mSm -1 . The resistive zones in the maps shown correspond to the bedrock in the study area (Fig. 1b). 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 near-surface shows significant spatiotemporal conductivity changes depending on the conductivity of the injected water. 355 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 mSm -1 range. The P2 plot shows larger spatial and temporal conductivity changes compared to P1, as expected, with conductivity varying between 30 and 80 mSm -1 . The conductivity of this zone decreases slightly on 14 Nov. This is probably due to the impact of heavy rain (31 mm) on 7 Nov, inducing salt leaching and reducing the soil conductivity. The P3 and P4 plots show stronger conductivity changes both spatially and temporally with conductivity varying 360 between 50 and 100 mSm -1 . The first models, obtained for 17 Oct, shows the minimum conductivity among the four selected dates. This is consistent with the simulated distributions shown in Fig. 6 and it is probably due to the saline front being still undispersed in the initial propagation phase (Fig. 5).
As the saline water is continuously added to the soil surface, the Clconcentration 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 Oct (Fig.  365 9, P3 and P4) are due to higher Clconcentration and its propagation to deeper layers. The conductivity models obtained for 14 Nov show a decrease in soil conductivity in both P3 and P4, although there were 6 irrigations after 23 Oct. 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. 6). This behaviour is probably due to the rainfall event on 7 Nov (as noted on P2) which diluted the Clconcentration in the soil. This suggests that the EMI surveys and data modelling detected well the expected change in Clconcentration. 370 Finally, the models obtained for the data collected on 23 Nov, depicted in Fig. 9, 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 4.1 and 4.2 show that the increase of Clconcentration in soil after 12 irrigation events, as well as the redistribution of Clduring the experiment period, are the main reason for increase in soil conductivity. Judging from the water content distribution maps, obtained from numerical simulations and shown in Fig. 4, we expect very small temporal variations 375 of water content in all experimental plots. Consequently, the temporal variations of Clconcentration and distribution is 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 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 380 sensor will be closely related to changes in salt concentrations.

Conclusion
In this study, we carried out a time-lapse EMI survey over four experimental plots irrigated with water at four different salinity levels during three 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 385 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 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 390 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 condition, allowing the optimization of data processing and 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. 395 2-A comprehensive investigation of our results and joint-interpretation of numerical simulations and time-lapse models reveals that the soil Clconcentration 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/noise ratio is usually better and EMI data will be more reliable; ii) the water content distribution will not change significantly, allowing to better study the impact 400 of soil salinity changes in time; iii) The sensitivity of σb to the Clconcentration 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 4-D space-time domain inversion algorithm was implemented in this study to invert entire time-lapse EMI data sets 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 405 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 410 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 415 methodology is strongly capable of giving information on 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 on 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 in-depth and 420 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. 425 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.