Estimating Pore Water Electrical Conductivity of Sandy Soil from Time Domain Reflectometry Records Using a Time-Varying Dynamic Linear Model

Despite the importance of computing soil pore water electrical conductivity (σp) from soil bulk electrical conductivity (σb) in ecological and hydrological applications, a good method of doing so remains elusive. The Hilhorst concept offers a theoretical model describing a linear relationship between σb, and relative dielectric permittivity (εb) in moist soil. The reciprocal of pore water electrical conductivity (1/σp) appears as a slope of the Hilhorst model and the ordinary least squares (OLS) of this linear relationship yields a single estimate (1/σp^) of the regression parameter vector (σp) for the entire data. This study was carried out on a sandy soil under laboratory conditions. We used a time-varying dynamic linear model (DLM) and the Kalman filter (Kf) to estimate the evolution of σp over time. A time series of the relative dielectric permittivity (εb) and σb of the soil were measured using time domain reflectometry (TDR) at different depths in a soil column to transform the deterministic Hilhorst model into a stochastic model and evaluate the linear relationship between εb and σb in order to capture deterministic changes to (1/σp). Applying the Hilhorst model, strong positive autocorrelations between the residuals could be found. By using and modifying them to DLM, the observed and modeled data of εb obtain a much better match and the estimated evolution of σp converged to its true value. Moreover, the offset of this linear relation varies for each soil depth.


Introduction
Salinization reduces crop productivity, decreases profitability, and causes land scarcity [1]. Thus, it decreases the world's agricultural productivity and causes a global income loss of US$ 12 billion per year [2]. Extracting soil solution by suction or using saturated paste conductivity measurements are the common methods to determine the electrical conductivity of soil pore water (σ p ) as an indicator of the soil salinity; however, they are labour-and cost intensive. There is no evidence that all ions are collected in the sample extract [3]. For soil salinity assessment, it is important to look for practical methods that evaluate the soil salinity state temporally and spatially. These methods help to correctly evaluate soil salinity evolution and reasonably predict its values [4][5][6][7][8][9]. In recent times, soil electromagnetic sensors have been used to estimate bulk electrical conductivity (σ b ). Then, methods are required to transform σ b to σ p [3,6,10].
where σ p is the pore water electrical conductivity (dS/m); ε p is the relative dielectric permittivity of the soil pore water (dimensionless), ε b is the relative dielectric permittivity of the bulk soil (dimensionless, relative dielectric permittivity is dimensionless since it is a ratio of permittivity of medium to the permittivity of free space), σ b is the bulk electrical conductivity (dS/m), ε σb=0 is the relative dielectric permittivity of the soil when the bulk electrical conductivity is 0 (dimensionless). However, ε σb=0 appears as an offset of the linear relationship between ε b and σ b . The Hilhorst model [3] concluded that his method could be validated for water contents between 0.10 and saturation and for a conductivity of the pore water up to 0.3 S m −1 . He found that ε σb=0 depends on soil type and varies between 1.9 and 7.6. He recommended using 4.1 as a generic offset. Many studies applied the deterministic Hilhorst model [3] in their experiments to convert σ b into σ p but they did not use the same offset to achieve their study objective. For example, some studies concluded their work by using different offsets (within the range of 3.67 to 6.38) according to the soil type [10]. The producer of capacitance soil moisture sensors 5TE [16] recommends the use of an offset ε σb=0 of 6 while another study found that an offset ε σb=0 = 6 does not present a good linear relationship between ε b and σ b [17]. The WET sensor (Delta-T Device Ltd., Cambridge, UK) is a frequency domain dielectric sensor. It has been designed to estimate the σ p based on the Hilhorst model [3] and incorporate the standard offset ε σb=0 = 4.1 of the model in the software of the device. By applying the Hilhorst model [3] in a saline gypsum-influenced soil, the accuracy of the WET sensor in predicting σ p was very poor when using the offset model =4.1 [18]. Another study used a WET sensor for experimental measurements in the laboratory using four different soils (sand, sandy loam, loam, and clay) [9] and found that the offset depends on both soil type and σ p , where it becomes larger for larger σ p . Moreover, oscillator frequency and sensor circuitry could affect the estimation of ε b and water content (θ) [19].
Sensors 2018, 18, x FOR PEER REVIEW 2 of 12 recent times, soil electromagnetic sensors have been used to estimate bulk electrical conductivity (σb). Then, methods are required to transform σb to σp [3,6,10]. According to Wyllie and Southwick [11], three conductance pathways, see Figure 1, contribute to the σb of a soil: (i) solid phase pathway through soil particles that are continuous contact with one another, (ii) liquid phase pathway through dissolved ions in the soil water inhabiting the large pores, and (iii) a liquid-solid interphase pathway through exchangeable cations like surfaces of clay minerals. Electrical conductivity (EC) in the liquid phase (σp) is used to estimate the soil salinity, a high EC refers to a high concentration of soluble salts, and vice versa. The σp could be estimated if the relationship between σp, σb, and water content (θ) is fixed [12][13][14]. The discovered linear correlation between the soil relative dielectric permittivity (εb) and σb values [15] enabled Hilhorst [3] to convert σb to σp by using a theoretical model. According to Hilhorst, σp can be determined from the equation: where σp is the pore water electrical conductivity (dS/m); εp is the relative dielectric permittivity of the soil pore water (dimensionless), εb is the relative dielectric permittivity of the bulk soil (dimensionless, relative dielectric permittivity is dimensionless since it is a ratio of permittivity of medium to the permittivity of free space), σb is the bulk electrical conductivity (dS/m), εσb=0 is the relative dielectric permittivity of the soil when the bulk electrical conductivity is 0 (dimensionless).
However, εσb=0 appears as an offset of the linear relationship between εb and σb. The Hilhorst model [3] concluded that his method could be validated for water contents between 0.10 and saturation and for a conductivity of the pore water up to 0.3 S m −1 . He found that εσb=0 depends on soil type and varies between 1.9 and 7.6. He recommended using 4.1 as a generic offset. Many studies applied the deterministic Hilhorst model [3] in their experiments to convert σb into σp but they did not use the same offset to achieve their study objective. For example, some studies concluded their work by using different offsets (within the range of 3.67 to 6.38) according to the soil type [10]. The producer of capacitance soil moisture sensors 5TE [16] recommends the use of an offset εσb=0 of 6 while another study found that an offset εσb=0 = 6 does not present a good linear relationship between εb and σb [17]. The WET sensor (Delta-T Device Ltd., Cambridge, UK) is a frequency domain dielectric sensor. It has been designed to estimate the σp based on the Hilhorst model [3] and incorporate the standard offset εσb=0 = 4.1 of the model in the software of the device. By applying the Hilhorst model [3] in a saline gypsum-influenced soil, the accuracy of the WET sensor in predicting σp was very poor when using the offset model =4.1 [18]. Another study used a WET sensor for experimental measurements in the laboratory using four different soils (sand, sandy loam, loam, and clay) [9] and found that the offset depends on both soil type and σp, where it becomes larger for larger σp. Moreover, oscillator frequency and sensor circuitry could affect the estimation of εb and water content ( ) [19]. There are three elementary causes responsible for why the deterministic system and control theories do not produce a totally sufficient means of performing this analysis and design: (i) many effects are left unknown since the objective of the model is to represent the main modes of system response, (ii) deterministic models are driven not by only our own control inputs but also disturbances which we can neither control nor model deterministically, and (iii) sensors do not offer exact readings of chosen quantities but present their own system dynamics and distortions as well and these devices are noise corrupted [20]. Despite the importance of computing σ p from σ b , a good method for doing so remains elusive (Campbell [16], personal communication).
Solute transport and water flow in the unsaturated zone are normally derived from the classical Richards equation and the convection-dispersion equations. Deterministic explanations of these equations are important aspects of research; due to soil heterogeneity at a variety of spatial scales, these equations for predicting actual field-scale processes are being increasingly questioned [21]. Therefore, some researchers working on soil heterogeneity concluded that for the evolution of soil water and solutes, it is more desirable to use stochastic models rather than constant values, where the parameters of stochastic transport models are treated as random variables with discrete values assigned according to a given probability distribution [21][22][23][24][25][26]. Among stochastic models, many studies used Kalman filtering in hydrological applications. A Kalman filter is an optimal recursive data processing algorithm that recursively couples the most recent measurements into the linear model to update the model state output [27]. Under the assumption that the linear system is a stochastic process with Gaussian noises, it produces the best estimation with minimum mean square error and it has been widely used in hydrological models to optimally merge information from the model simulations and the independent observations with appropriate modeling [28][29][30][31][32].
In previous work, we installed frequency domain reflectometry (FDR) sensors (5 TE), which are commercially available from METER Group, Inc. USA, in field conditions at different depths where the soil is heterogeneous to estimate σ p [22]. We used ε b and σ b observations to modify the Hilhorst deterministic model [3] to a stochastic model using a time-varying dynamic linear model and Klaman filter before studying the linear relationship between them.
In this study, we used Time Domain Reflectometry (TDR) sensors (FP/mts), which are commercially available from Easy Test, Poland, to measure ε b and σ b in laboratory conditions where the soil is homogeneous. Then, we tried to use the Hilhorst model [3] to convert σ b to σ p . Later, we could show the weakness of applying the deterministic Hilhorst model [3] even in homogeneous soils. Thus, we are aiming to adapt this approach to a stochastic model under laboratory conditions. Thus, we used one homogeneous soil type to accurately estimate the changes in σ p over time and to conclude whether the model offset is constant or if it changes in one soil profile.

The Column Experiment
To achieve the objective of this study, we used two soil columns with a height of 55 cm provided by a sprinkler, see Figure 2. The lower boundary was controlled using a vacuum pump at a constant pressure head of −30 hPa. The columns were packed with a density of 1.4 g/cm 3 . The substrate was sand, 80% of which was fine sand. The water content during packing was approximately 4 m 3 /m 3 . The TDR and soil temperatures sensors were installed in four depths: 7, 21, 35, and 48 cm. Since the soil is sand, the soil relative dielectric permittivity (ε b ), bulk electrical conductivity (σ b ), and temperature were measured every 5 min to obtain enough observations for modeling. Additionally, porous suction cups for taking soil solution samples were installed at each depth to validate the results of our model. The lower boundary of the column uses a membrane to let the water drain. Drainage water was collected in a bottle under −30 hPa vacuum, which is supplied in the range from −20 to −30 hPa. The sprinkler is 5 cm above the soil surface and allows water to drop through small nozzles. Five irrigation events using KCl solution with different electrical conductivities were applied. The first three events were irrigated with 20 dS/m of KCl, then the fourth, and fifth events with 30 dS/m of KCl. The flux was approximately 1 l/h. The columns were free of salt at the beginning and before the irrigations events started. The TDR probes are FP/mts commercially available from Easy Test, Poland, and have been calibrated in air and deionized water. The temperature probes are Thermistors of the type 2k252 (type Fenwal UUA 32J3) with a range of −20 up to 60 °C. Soil temperature data (Tsoil) were used to estimate the relative dielectric permittivity of the soil pore water directly (εp): To apply the dynamic linear model and the Kalman filter, a time series of the variable of interest is needed [27]. In our study, time series of εb, σb, and εp are required to estimate σp. Therefore, we used five irrigation events with two levels of KCL solution to obtain the variation of these variables over time for each depth. In total, 289 observations were made of σb, εp, and εb for each soil depth and these were used to estimate both the offset εσb=0 of the modified Hilhorst model [3] and the evolution of σp at its corresponding depth, of which 144 observations were used to validate their forecasts.

Time-Varying Dynamic Linear Model
In general, the state space model is identified by two assumptions, (i) there is a hidden or latent process xt called the state process. The state process is assumed to be a Markov process, where past and future values of xt are independent, conditional on the present xt, ({xs, S > t}, and {xs, S< t} are independent on the xt), (ii) the observations, yt are independent given the states xt. This means that the dependence among the observations is generated by states. The dynamic linear model (DLM) or Additionally, porous suction cups for taking soil solution samples were installed at each depth to validate the results of our model. The lower boundary of the column uses a membrane to let the water drain. Drainage water was collected in a bottle under −30 hPa vacuum, which is supplied in the range from −20 to −30 hPa. The sprinkler is 5 cm above the soil surface and allows water to drop through small nozzles. Five irrigation events using KCl solution with different electrical conductivities were applied. The first three events were irrigated with 20 dS/m of KCl, then the fourth, and fifth events with 30 dS/m of KCl. The flux was approximately 1 l/h. The columns were free of salt at the beginning and before the irrigations events started. The TDR probes are FP/mts commercially available from Easy Test, Poland, and have been calibrated in air and deionized water. The temperature probes are Thermistors of the type 2k252 (type Fenwal UUA 32J3) with a range of −20 up to 60 • C. Soil temperature data (Tsoil) were used to estimate the relative dielectric permittivity of the soil pore water directly (ε p ): To apply the dynamic linear model and the Kalman filter, a time series of the variable of interest is needed [27]. In our study, time series of ε b , σ b , and ε p are required to estimate σ p . Therefore, we used five irrigation events with two levels of KCL solution to obtain the variation of these variables over time for each depth. In total, 289 observations were made of σ b , ε p , and ε b for each soil depth and these were used to estimate both the offset ε σb=0 of the modified Hilhorst model [3] and the evolution of σ p at its corresponding depth, of which 144 observations were used to validate their forecasts.

Time-Varying Dynamic Linear Model
In general, the state space model is identified by two assumptions, (i) there is a hidden or latent process x t called the state process. The state process is assumed to be a Markov process, where past and future values of x t are independent, conditional on the present x t , ({x s , S > t}, and {x s , S< t} are independent on the x t ), (ii) the observations, y t are independent given the states x t . This means that the dependence among the observations is generated by states. The dynamic linear model (DLM) or linear Gaussian state space model, in its simple form, employs a first-order, p-dimensional vector autoregression as the state equation: Sensors 2018, 18, 4403

of 12
We do not observe the state vector x t directly, but only a linear transformed version of them with noise added, say: y t is an m-dimensional vector, representing the observation at time t, A t is a q × p measurement or observation matrix. Equation (4) is called the observation equation, in which v t , w t are the Gaussian white-noise errors. The evolution variances are V t , W t and can be estimated from available data using maximum likelihood or Bayesian techniques. In this study, we modified the deterministic Hihlorst model (1) to a stochastic one. The model (1) has the variables σ b , ε p , σ p , ε b , and ε σb = 0 . The σ p and ε σb = 0 are unobserved and they need to be estimated by the state Equation (3) as x t , while σ b , ε p , and ε b are observed by the sensors (ε p is calculated from Equation (2) using soil temperature sensor data) and represented by observation Equation (4) as y t .
The R [33] package, dlm [34], provides an integrated environment for Bayesian inference using DLM, and the package includes functions for Kalman filtering and smoothing, as well as for maximum likelihood estimation.

Deterministic Model
The offset of the Hilhorst model [3] can be calculated from Equation (1): We derived the offset (ε σb = 0 ) from this linear model after using measurements of ε b and σ b . For example, applying the ordinary least squares (OLS) on measurements of ε b and σ b obtained from soil column 2 data during the third irrigation at a depth of 21 cm, Table 1 shows that the offset of the linear relationship between ε b -σ b is 9.41. Further, the single estimate of the slope ( 1/σ p ) of the regression parameter vector (1/σ p ) for the entire data set is very small. Thus, the estimated soil pore water electrical conductivity (σ p ) is too high compared with the EC meter value, see Table 2. Afterward, we applied the Durbin-Watson test in order to test if there was any autocorrelation between the residuals of the regression. Table 3 shows that there is an extremely strong and positive autocorrelation, meaning that the result of that regression is not valid.   For each irrigation event, we obtained one solution sample at each depth by using porous suction cups. Unfortunately, some samples did not have enough solution to measure their electrical conductivity using the EC meter device. Table 2 shows the values of EC measured by the EC meter device. The table shows eight EC values from the EC meter according to the depth and irrigation event number for each soil column. Due to the variability in the water flow in unsaturated soil, we observed in our experiment a variation in the time needed to collect the solution sample. More time was required to collect enough solution for the EC meter device when a greater number of ions gathered in the sample resulting in a high EC value of the sample. Therefore, there is a difference in the EC values between the soil columns at the same depth, see Table 2. We applied a modified Hilhorst model on the eight time-series data corresponding to Table 2 (depth, irrigation event, and soil column) to compare our finding of σ p obtained from our modified model to the values of σ p obtained by the EC meter device, see Table 2.
The reason for choosing 1 l/h for the irrigation rate and 5 min for the irrigation interval is visualized in Figure 3. At each depth, we could see how the bulk electrical conductivity responds to the irrigation event.  For each irrigation event, we obtained one solution sample at each depth by using porous suction cups. Unfortunately, some samples did not have enough solution to measure their electrical conductivity using the EC meter device. Table 2 shows the values of EC measured by the EC meter device. The table shows eight EC values from the EC meter according to the depth and irrigation event number for each soil column. Due to the variability in the water flow in unsaturated soil, we observed in our experiment a variation in the time needed to collect the solution sample. More time was required to collect enough solution for the EC meter device when a greater number of ions gathered in the sample resulting in a high EC value of the sample. Therefore, there is a difference in the EC values between the soil columns at the same depth, see Table 2. We applied a modified Hilhorst model on the eight time-series data corresponding to Table 2 (depth, irrigation event, and soil column) to compare our finding of σp obtained from our modified model to the values of σp obtained by the EC meter device, see Table 2.
The reason for choosing 1 l/h for the irrigation rate and 5 min for the irrigation interval is visualized in Figure 3. At each depth, we could see how the bulk electrical conductivity responds to the irrigation event.

Time-Varying Linear Dynamic Model (LDM)
The deterministic Equation (5) can be modified into the time-varying DLM for observation and unobservable (state) models. In this case, the observation data are the soil relative dielectric permittivity (εb), bulk electrical conductivity (σb), and the relative dielectric permittivity (εp), while the unobservable data are the offset (εσb=0) and pore water electrical conductivity (σp). Equation (4) can be modified to the time-varying DLM as follows: • The observation equation can be obtained by modifying the Hilhorst model [3] (written in Equation (5)) into a stochastic equation, in accordance with Equation (4) as follows:

Time-Varying Linear Dynamic Model (LDM)
The deterministic Equation (5) can be modified into the time-varying DLM for observation and unobservable (state) models. In this case, the observation data are the soil relative dielectric permittivity (ε b ), bulk electrical conductivity (σ b ), and the relative dielectric permittivity (ε p ), while the unobservable data are the offset (ε σb=0 ) and pore water electrical conductivity (σ p ). Equation (4) can be modified to the time-varying DLM as follows: • The observation equation can be obtained by modifying the Hilhorst model [3] (written in Equation (5)) into a stochastic equation, in accordance with Equation (4) as follows: Here, we consider ε σb = 0 as a constant. The actual value is related only to its past value. The slope 1/σ p changes over time and its actual value is related to its past value plus the Gaussian white-noise errors (w t ). We applied the equation in reverse order to estimate the state variables (ε σb = 0 and σ p ) at all time points from a complete series of the soil relative dielectric permittivity (ε b ). This process is known as smoothing.
An example of the evolution of ε b , ε p , and σ b data needed for the Hilhorst model [3] is shown in Figure 4. By applying the Equations (6) and (7) using DLM and the Kalman filter on the eight time-series data, we see in Figure 5 the observed and predicted time series of the soil relative dielectric permittivity (ε b ). The predicted and observed values of ε b agree reasonably well. The mean absolute prediction error (MAPE) for the time series never exceeded 0.02.
The state equation (unobservable data) in Equation (3) Here, we consider εσb = 0 as a constant. The actual value is related only to its past value. The slope 1/σp changes over time and its actual value is related to its past value plus the Gaussian white-noise errors ( ). We applied the equation in reverse order to estimate the state variables (εσb = 0 and σp) at all time points from a complete series of the soil relative dielectric permittivity (εb). This process is known as smoothing.
An example of the evolution of εb, εp, and σb data needed for the Hilhorst model [3] is shown in Figure 4. By applying the Equations (6) and (7) using DLM and the Kalman filter on the eight timeseries data, we see in Figure 5 the observed and predicted time series of the soil relative dielectric permittivity (εb). The predicted and observed values of εb agree reasonably well. The mean absolute prediction error (MAPE) for the time series never exceeded 0.02.  Since the prediction of the soil relative dielectric permittivity (ε b ) is valid, the estimation of the electrical conductivity of pore water (σ p ) and the offset ε σb = 0 , see Equation (7), are also valid because they are used in the prediction of the soil relative dielectric permittivity (ε b ) and have converged to their true values. The evolution of σ p over time obtained by DLM is presented in Figure 6; it shows the importance of using DLM because it obtained all the changes of σ p over time and not a single value of σ p for the entire data set. Another interesting aspect is that Figure 6 shows the changes in the model offset for each irrigation event at each depth. This finding is very important since it shows that the offset does not depend on the soil type [10,[14][15][16] nor on the soil type and salinity [17] when two columns with the same type soil are used, as in this study. Moreover, in Figure 6 we put the corresponding value of σ p measured by the EC meter device for each depth according the irrigation event and soil column number.
Comparing the mean evolution of σ p values obtained from our modified Hilhorst model, see Figure 6, with the single corresponding EC value obtained from porous suction cups and measured by the EC meter device, see Table 2, we found that they agree very well (R 2 = 72%). From these results, three advantages are evident when using DLM and a Kalman filter to estimate σ p in two homogenous soil columns; first, we observed that the offset value of the Hilhorst model does not depend on the soil type and σ p and it changes in the same soil profile. Secondly, we obtained the changes in the estimated σ p over time and not just a single value as a coefficient for the entire data set. Third, the estimated changes in σ p occur instantly and save time and labor costs.

Conclusions
In this study, we applied the ε b -σ b linear relationship to homogeneous soil column data obtained from TDR sensors. We found an extremely strong positive autocorrelation between the residuals of the regression analysis. When residuals are correlated, the least squares method is not the most efficient model coefficient estimator. By modifying the regression by a time-varying dynamic linear model (DLM), the match between the observed and modeled data of ε b is significantly improved and the estimated evolution of σ p converges to its true value. Moreover, in this study, we used two homogeneous soil columns with the same condition to show that the offset of the Hilhorst model [3] is not constant, as suggested for all moist soil or, as others suggested, that it is soil-type-dependent [10,[14][15][16] or soil-type-and salinity-dependent [17]. We repeated the experiment to show that the offset changes even in the same soil type and the same conditions. A dynamic linear model enables the capture of the offset changes and it shows the importance of calculating it simultaneously when estimating σ p using the Hilhorst model. The next promising step would be programming and inserting these models into the TDR software in order to estimate the soil pore water electrical conductivity (σ p ) from senor records directly.