A New Method to Determine the Spatial Sensitivity of Time Domain Reﬂectometry Probes Based on Three-Dimensional Weighting Theory

: The time domain reﬂectometry (TDR) method has been widely used to measure soil water content for agriculture and engineering applications. Quick design and optimization of the probe is crucial to achieving practical utilization. Generally, the two-dimensional weighting theory, calculation of the spatial sensitivity of TDR probes in the plane transverse to the direction of electromagnetic wave propagation, and relevant numerical simulation techniques can be used to solve any issues. However, it is di ﬃ cult to tackle speciﬁc problems such as complex probe shape, end e ﬀ ect, and so forth. In order to solve these issues, a method including a three-dimensional weighting theory and the relevant numerical simulation technique was proposed and veriﬁed to conﬁrm the feasibility of this method by means of comparing the existing experimental results and the computational values. First, a shaft probe was used to determine the impact of the shaft on the e ﬀ ective dielectric constant of the probe. Then, three-rod probes were calibrated by a sample with a special shape and water-level variations around the probe using the proposed method to determine the values of the apparent dielectric constant. Besides, model boundary size and end e ﬀ ect were also considered in the computation of dielectric constants. Results showed that compared with the experimental and computational data, the newly proposed method calculated the measurement sensitivity of the shaft probes well. In addition, it was observed that experiment dielectric constant values were slightly di ﬀ erent from computational ones, not only using a vertical probe but also horizonal probe. Moreover, it was also found that there was a slight inﬂuence of sample shape and end e ﬀ ect on the apparent dielectric constant, but model boundary size has a certain impact on the values. Overall, the new method can provide beneﬁts in the design and optimization of the probe. the RMSE values between the experimental data of the probe B and data computed by 3D or 2D models were 1.609 and 2.115, respectively, which indicates that data obtained by the 3D models were more representative of the experimental data than those derived by the 2D models. This may be because the 3D method can encompass the spatial e ﬀ ect, which does not occur in the 2D method. In this regard, the 3D method has an advantages over the 2D method.


Introduction
The time-domain reflectometry (TDR) technique has been widely applied to measure soil water content in agriculture and engineering practice. A number of TDR probes have been designed to satisfy special circumstances, so probe patterns are diverse and include the rod probe, plane probe, shaft probe, and noninvasive probe, etc. [1][2][3][4]. Meanwhile, there are many methods for the design and optimization of the probes. These methods can be roughly separated into two classes: experimental calibration and numerical computation. Numerical computation can quickly compute and analyze the feasibility in the process of designing the probes. Whether it works or not depends on the theory used to support the computation. However, as per the author's understanding, there are still some issues that cannot be solved by existing theory, such as complex probe shape that cannot be simplified to a plane, nonuniformity or layering of the soil medium, and special sample shape; thus, new theories are required to solve these issues.
Initially, a balanced transmission line with two conductors carrying opposite polarities was used to examine the voltage distribution around the rods [5]. Then, the theory of the sampling volume dependent on the weighting of the energy density distribution was proposed in [6] based on approximate analytical solutions of two-rod [6] and multi-rod [7] probes. In addition, the authors of [8] used numerical methods to solve the Laplace's equation in a two-dimensional model to determine the distribution of electric potentials and the energy density of TDR probes. The numerical method has been used to determine the region of porous material and fraction contributing to the TDR measurement that can be used to represent the sampling area [9]. It has also been used to determine some other probe configurations [10,11]. The authors of [12] used the two-dimensional weighting theory and the Geostudio software commonly used in hydrology to compute the apparent dielectric constant of the shaft probe. The total weighting factors consist of two parts: the weighting factor of the shaft and the weighting factor of the target medium. However, the two-dimensional weighting theory can only consider the computational plane and neglects the effects of the rest of probe and sample and so forth. For instance, the diameter variation of the probe rod and an irregular shape affect the response of the probe. The software used to directly compute the electric field utilizes an unambiguous theory to calculate the weight factors. Therefore, a new method based on three-dimensional (3D) weighting theory and software that can compute the electric field is proposed to simultaneously cover the effects of whole probe and sample on the response of the probe.
The investigation aims to assess the performance of the new method, combining three-dimensional weighting theory with numerical computation of electrostatic field, by means of comparing the computational values of the dielectric constant and the existing experimental data of TDR probes. In detail, the two-dimensional spatial weighting function was first extended to a three-dimensional one, followed by calculation of the spatial sensitivity of the shaft and rod probes that had been calibrated in the experimental studies. Besides, a comparison of the calculated and experimental results was conducted to validate the feasibility of the three-dimensional weighting theory.

Theory
As shown in Figure 1, a TDR device generally comprises a pulse generator, a sampler, an oscilloscope, a coaxial cable, and a probe. A step voltage pulse signal generated in the pulse generator propagates in an electromagnetic wave along the coaxial cable and is reflected at both the beginning and end of a probe that is normally in contact with the measured medium. The signal is sampled by the sampler and displayed by the oscilloscope. The travel time (∆t) can be measured from the reflected waveform. The apparent dielectric constant K a of the testing medium is a function as the travel time ∆t, the velocity of light c (3.0 × 10 8 m/s), and the length of probe in the testing medium. In addition, some planes and axes, such as planes transverse to the direction of electromagnetic wave propagation, the probe center plane, and the long axes of the rods, are plotted in the Figure 1 and will be involved in this work.
The spatial weighting concept proposed by [6] and the numerical approach of [8] are commonly used to determine the sample volume and distribution of energy density in the plane perpendicular to TDR rods that contributes significantly to the response of a TDR probe. The authors of [6,9] defined a spatial weighting factor (w(x, y)) at each point for the sensitivity of TDR probes, which can be expressed in Equation (1) as follows: w(x, y) = ∇Φ(x, y) 2 / where Φ(x, y) is the electrical potential in the heterogeneous field, Φ 0 (x, y) is the electrical potential in the equalized homogeneous field, Ω is the domain of integration, and dA is the area of the element [8]. The spatial weighting function concept underlying this analytical solution was extended to construct a numerical approach capable of predicting the response of TDR probes of any geometry to any spatial distribution of relative dielectric permittivity in the plane transverse to the direction of electromagnetic wave propagation. Ignoring end effects, the three-dimensional sample volume was defined as the projection of this two-dimensional sample area along the length of the rods [9]. However, the plane perpendicular to the long axis of the rods sometimes cannot represent the practical conditions. A three-dimensional spatial weighting factor (w(x, y, z)) at each point can be expressed in Equation (2) as follows: where Φ(x, y, z) is the electrical potential in the heterogeneous field, Φ 0 (x, y, z) is the electrical potential in the equalized homogeneous field, Σ is the domain of integration, and dV is the volume of the element. The predicted apparent relative dielectric permittivity (K a ) is then, where K(x, y, z) is the dielectric constant of a point (x, y, z) in the domain Σ.
Water 2020, 12, x 3 of 14 The spatial weighting concept proposed by [6] and the numerical approach of [8] are commonly used to determine the sample volume and distribution of energy density in the plane perpendicular to TDR rods that contributes significantly to the response of a TDR probe. The authors of [6,9] defined a spatial weighting factor (w(x, y)) at each point for the sensitivity of TDR probes, which can be expressed in Equation (1) as follows: where Φ(x, y) is the electrical potential in the heterogeneous field, Φ0(x, y) is the electrical potential in the equalized homogeneous field, Ω is the domain of integration, and dA is the area of the element [8]. The spatial weighting function concept underlying this analytical solution was extended to construct a numerical approach capable of predicting the response of TDR probes of any geometry The measurement sensitivity can be defined as the contribution level of a medium to the apparent relative dielectric permittivity. The spatial weighting function describes the sensitivity of the instrument to any given point in the domain. The sensitivity of the whole domain is defined in Equation (4) and is equal to 1. For uniform K a the electrostatic potential satisfies Laplace's Equation (5), which can be used to calculate electrical potential in the field,

Shaft Probes
Numerical analyses were conducted to investigate the spatial sensitivity of two types of TDR probes. According to [12], Type A has a shaft separating two conductors and Type B has four straight copper strips embedded on the nonconductor shaft surface with configurations similar to the probes reported by [13,14]. Both types of the probes have a diameter of 30 mm and the details are presented in [12]. The length of the two probes was 200 mm. Then, the electromagnetic field was calculated by a seepage analysis software and the effective apparent dielectric permittivity of the simulated medium was obtained based on the spatial sensitivity equations by [12].
In order to verify the reliability of the proposed methods including the three-dimensional spatial weighting function theory and numerical simulation technique, we computed the electric field distributions inside the sample using the finite element modeling (FEM) algorithm as implemented in the Comsol software [15], obtaining the distribution of electrostatic potential surrounding the probe, and analyzed the measurement sensitivity along the rods. Then, the sensitivity results calculated by the Comsol software were compared with published data using other numerical techniques reported by [12]. In order to verify the feasibility and reliability of the numerical simulation technique used in this work, two-dimensional modeling was constructed to compute the fields with the same model dimensions with that reported by [12], but it should be noted that the numerical simulation software and technique were different. Moreover, based on the three-dimensional spatial weighting factor (w(x, y, z)) and the measurement sensitivity equation for a TDR shaft probe, three-dimensional numerical analyses were conducted to investigate the spatial sensitivity for the two types of TDR shaft probes, and the modeling and computation details are presented below.
The Comsol multiphysics package has graphic user interfaces (GUIs) [15], which can develop complex geophysical models with simple geometries such as blocks and cylinders, and the meshes or grids required for representing the fields can be generated easily. Two representative shaft probes, type A and B probes, were selected to verify the proposed method. Using the 3D type B TDR shaft probe as an example, the transparent view of 3D numerical model is presented in Figure 2. The model was established in terms of the dimensions mentioned in Figure 2 including conductors and shaft. M1 and M2 respectively represent the testing targeted medium and shaft material, while M3 represents the four conductors. S1, S2, S3, and S4 represent the lateral boundaries of the four conductors and S5 is representative of the lateral boundaries of the targeted medium. S1 and S3 were set to a constant voltage value of +1 V, while S2 and S4 were set to a constant voltage value of −1 V. Meanwhile, S5 was set as a zero charge boundary that indicates that the boundary corresponds to air. The dielectric constants of the targeted medium varied from 10 to 80. The dielectric constants of the shaft in type A and B probes were 5.9699 and 3.5477, respectively. The details of configurations can be found in Figure 2. In the entire domains, we used free tetrahedral elements to discretize the models. The length ratio of the side of the mesh should be close to 1. After building and meshing the models, Comsol generated the partial differential equations and solved the equations on the mesh. We chose the damped Newton method to solve the nonlinear problem.
Since a number of probe designs incorporate nonmetallic probe components of known K a that lie within the sampling volume of the probe, it is important to know how these materials affect the spatial weighting function to compare the performance of probes and to optimize their designs [9,12]. A two-dimensional equation was proposed considering that the value of the effective apparent dielectric permittivity measured using the shaft probe is a weighted average of the apparent dielectric permittivity of the surrounding soil (i.e., targeted medium) and the apparent dielectric permittivity of Water 2020, 12, 545 5 of 14 the nonconductor shaft. According to Equation (3), the effective apparent dielectric permittivity K a,eff can be expressed in Equation (6) as follows: where K a,o and K a,t are the apparent dielectric permittivity of the shaft and targeted medium, respectively, and w o (x, y, z) and w t (x, y, z) are the weighting factors that are distributed in the shaft zone and the testing targeted medium zone, respectively.
Water 2020, 12, x 5 of 14 of the nonconductor shaft. According to Equation (3), the effective apparent dielectric permittivity Ka,eff can be expressed in Equation (6) as follows: K , = K , w x, y, z dV K , w x, y, z dV (6) where Ka,o and Ka,t are the apparent dielectric permittivity of the shaft and targeted medium, respectively, and wo(x, y, z) and wt(x, y, z) are the weighting factors that are distributed in the shaft zone and the testing targeted medium zone, respectively. After computation of the energy density field, the weighting factors of the shaft and the target medium could be calculated according to Equation (2). Then, the effective apparent dielectric permittivity dielectric (also named effective dielectric constant) was calculated by Equation (6) in terms of the known dielectric constant of the shaft and the target medium.

Sample Calibration
Generally, probes are calibrated by the soil medium under the condition that the probes are orthogonally vertical to the surface of the medium, but cylindrical samples have sometimes been used to calibrate probes that are penetrated into the arc lateral face [16]. Probe calibration was performed in an acrylic tubing (7.62 cm inner diameter by 8.89 cm outer diameter) and the probe was laterally extended into the medium placed in the tube. In order to measure the variation of volumetric water content of soil cylinder for determining its unsaturated water conductivity in the laboratory, the authors of [17] used a small cylindrical stiff clay sample with a diameter of 38 mm, a height of 40 mm, and a dry density of 2.0 g/cm 3 to calibrate a miniature probe, and the three-rod probe was radially penetrated into the sample with the probe plane paralleling with the axis of the cylindrical sample. Verification of the proposed method was performed based on 3D weighting theory and the Comsol modeling technique using the shaft probe. In this section, the effect of sample shape on the relationship of dielectric constant Ka was quantitatively computed and analyzed based on the measurement sensitivity of the probes developed in the present work. In detail, probes with different For the type B probe, M1 and M2 represent the testing medium and shaft material, and M3 involves the four conductors. S1, S2, S3, and S4 represent the lateral boundaries of the four conductors (M3), respectively. S5 represents the lateral boundaries of the testing medium. The unit of the dimensions of the model is mm.
After computation of the energy density field, the weighting factors of the shaft and the target medium could be calculated according to Equation (2). Then, the effective apparent dielectric permittivity dielectric (also named effective dielectric constant) was calculated by Equation (6) in terms of the known dielectric constant of the shaft and the target medium.

Sample Calibration
Generally, probes are calibrated by the soil medium under the condition that the probes are orthogonally vertical to the surface of the medium, but cylindrical samples have sometimes been used to calibrate probes that are penetrated into the arc lateral face [16]. Probe calibration was performed in an acrylic tubing (7.62 cm inner diameter by 8.89 cm outer diameter) and the probe was laterally extended into the medium placed in the tube. In order to measure the variation of volumetric water content of soil cylinder for determining its unsaturated water conductivity in the laboratory, the authors of [17] used a small cylindrical stiff clay sample with a diameter of 38 mm, a height of 40 mm, and a dry density of 2.0 g/cm 3 to calibrate a miniature probe, and the three-rod probe was radially penetrated into the sample with the probe plane paralleling with the axis of the cylindrical sample. Verification of the proposed method was performed based on 3D weighting theory and the Comsol modeling technique using the shaft probe. In this section, the effect of sample shape on the relationship of dielectric constant K a was quantitatively computed and analyzed based on the measurement sensitivity of the probes developed in the present work. In detail, probes with different diameters (0.55 mm, 2.55 mm, and 4.55 mm), a net separation of 4.5 mm, and a length of 30 mm were used. The diameter of the sample varied from 30 to 80 mm. Figure 3 shows the transparent view of the 3D numerical model for computing the electrical potential field surrounding the TDR probe based on experimental studies reported by the authors of [17]. The model plotted in Figure 3 comprises five domains: one of them is of soil and the rest of domains are of air. Upon simulation, the model was segmented into the meshes that were constituted by free tetrahedral elements. The length ratio of the side of the mesh should be close to 1. The boundary condition and rod boundary voltage were set to be the same as in the type B model in Figure 2. The dielectric constant of air was equal to 1 and the dielectric constant of the sample ranged from 10 to 80 in the computation. After computation, the apparent dielectric constants were calculated from the dielectric constant of the air and the sample based on the Equation (3). In other words, the 3D weighting function w(x, y, z) was first calculated based on the numerical simulation technique mentioned above. Then, the integration was performed in the domain of air and sample in terms of the Equation (3) to obtain the apparent dielectric constants. There were five groups of models and they were numbered in terms of sample and probe rod diameter. For example, for DS30R0.55 the diameters of the sample and probe rod in the model were 30 mm and 0.55 mm, respectively.
Water 2020, 12, x 6 of 14 diameters (0.55 mm, 2.55 mm, and 4.55 mm), a net separation of 4.5 mm, and a length of 30 mm were used. The diameter of the sample varied from 30 to 80 mm. Figure 3 shows the transparent view of the 3D numerical model for computing the electrical potential field surrounding the TDR probe based on experimental studies reported by the authors of [17]. The model plotted in Figure 3 comprises five domains: one of them is of soil and the rest of domains are of air. Upon simulation, the model was segmented into the meshes that were constituted by free tetrahedral elements. The length ratio of the side of the mesh should be close to 1. The boundary condition and rod boundary voltage were set to be the same as in the type B model in Figure 2. The dielectric constant of air was equal to 1 and the dielectric constant of the sample ranged from 10 to 80 in the computation. After computation, the apparent dielectric constants were calculated from the dielectric constant of the air and the sample based on the Equation (3). In other words, the 3D weighting function w(x, y, z) was first calculated based on the numerical simulation technique mentioned above. Then, the integration was performed in the domain of air and sample in terms of the Equation (3) to obtain the apparent dielectric constants. There were five groups of models and they were numbered in terms of sample and probe rod diameter. For example, for DS30R0.55 the diameters of the sample and probe rod in the model were 30 mm and 0.55 mm, respectively.

Water Calibration
As shown in Figure 4 [18], a three-rod probe was calibrated with a sensing electrode length of 150 mm. The dimensions of the rods had a diameter of 6 mm and the spacing between the center of the electrodes was 30 mm. The probe was fitted to the center of the short side of a rectangular organic glass container with a height of 100 mm, a width of 100 mm, and a length of 200 mm. The plane that the three rods lay in was parallel to the horizontal plane. Initially, the container was full of water and the water surface had a vertical distance of 50 mm to the plane that the center of the probes lay in. The level of water surface was reduced 2 mm each time by draining water out of the container, and the dielectric constant of the TDR probe was recorded. This procedure was repeated until the level of the water surface was reduced to the level of the plane that the center of the probes lay in. During the simulation, the dielectric constants of water and air were 80 and 1, respectively.

Water Calibration
As shown in Figure 4 [18], a three-rod probe was calibrated with a sensing electrode length of 150 mm. The dimensions of the rods had a diameter of 6 mm and the spacing between the center of the electrodes was 30 mm. The probe was fitted to the center of the short side of a rectangular organic glass container with a height of 100 mm, a width of 100 mm, and a length of 200 mm. The plane that the three rods lay in was parallel to the horizontal plane. Initially, the container was full of water and the water surface had a vertical distance of 50 mm to the plane that the center of the probes lay in. The level of water surface was reduced 2 mm each time by draining water out of the container, and the dielectric constant of the TDR probe was recorded. This procedure was repeated until the level of the water surface was reduced to the level of the plane that the center of the probes lay in. During the simulation, the dielectric constants of water and air were 80 and 1, respectively.
In addition, the authors of [19] presented laboratory measurements of the effects of steep gradients in relative dielectric permittivity on the spatial sensitivity of TDR probes, in which three-rod probes were placed horizontally and vertically within an experimental box. The inner dimensions of the box were 0.3 × 0.3 × 0.3 m. Then, an air-liquid interface was raised upwards from below and past the probe. The effective length of the probe rods was approximately 285 mm. The three-rod probe with a diameter equal to 4.6 mm and outer rod spacing equal to 70.0 mm was inserted through the box walls.
In this work, the laboratory measurement obtained by the three-rod probes was used for comparison with the computational results. In addition, the authors of [19] presented laboratory measurements of the effects of steep gradients in relative dielectric permittivity on the spatial sensitivity of TDR probes, in which threerod probes were placed horizontally and vertically within an experimental box. The inner dimensions of the box were 0.3 × 0.3 × 0.3 m. Then, an air-liquid interface was raised upwards from below and past the probe. The effective length of the probe rods was approximately 285 mm. The three-rod probe with a diameter equal to 4.6 mm and outer rod spacing equal to 70.0 mm was inserted through the box walls. In this work, the laboratory measurement obtained by the three-rod probes was used for comparison with the computational results.
The computational models used in this work considering different computational zones at the initial state are shown in Figure 5. There were four models (N1, N2, N3, and N4) established to study the impact of water or air on the dielectric constant. In Figure 5a the model N1 was constructed considering water in the container and air surrounding the water, and Figure 5b shows the model N2 without considering the air around the container, while models N3 and N4 plotted in Figure 5c,d involved the end effect including electric response in the water medium near the probe end and sharp probe end with a length of 20 mm. Note that these four models are in the initial state. When the water surface level falls, the air medium will replace the initial water medium. Model N1 is taken as an example for introducing the model computed by software Comsol [15]. As shown in Figure 5a, M1 is representative of water zone with an initial height of 100 mm, width of 100 mm, and length of 200 mm. M2 represents air with an initial height of 200 mm and width of 200 mm that the organic glass wall of the container does not consider in the model. The boundary B1 of the rod in the middle of the probe was set with the voltage of +1 V, while the boundary B2 of the rods on either side of the probe was set with the voltage of −1 V. The outer boundary of the overall model was of zero charge, while the medium of M1 and M2 confirmed charge conservation. As the water surface fell, the water was replaced by the air until the level of the water surface level fell to the level of the plane that the center of the probes lay in, i.e., half of water in the container was replaced by the air medium. The model mesh, computation, and analysis method were the same as in the model used in the sample calibration section. The computational models used in this work considering different computational zones at the initial state are shown in Figure 5. There were four models (N1, N2, N3, and N4) established to study the impact of water or air on the dielectric constant. In Figure 5a the model N1 was constructed considering water in the container and air surrounding the water, and Figure 5b shows the model N2 without considering the air around the container, while models N3 and N4 plotted in Figure 5c,d involved the end effect including electric response in the water medium near the probe end and sharp probe end with a length of 20 mm. Note that these four models are in the initial state. When the water surface level falls, the air medium will replace the initial water medium. Model N1 is taken as an example for introducing the model computed by software Comsol [15]. As shown in Figure 5a, M1 is representative of water zone with an initial height of 100 mm, width of 100 mm, and length of 200 mm. M2 represents air with an initial height of 200 mm and width of 200 mm that the organic glass wall of the container does not consider in the model. The boundary B1 of the rod in the middle of the probe was set with the voltage of +1 V, while the boundary B2 of the rods on either side of the probe was set with the voltage of −1 V. The outer boundary of the overall model was of zero charge, while the medium of M1 and M2 confirmed charge conservation. As the water surface fell, the water was replaced by the air until the level of the water surface level fell to the level of the plane that the center of the probes lay in, i.e., half of water in the container was replaced by the air medium. The model mesh, computation, and analysis method were the same as in the model used in the sample calibration section. Figure 6 shows the model N5 for simulation of calibration of vertical probe in terms of the report by [19]. The probe had a diameter of 4.6 mm and an outer separation of 70 mm. The length of the model was 285 mm, which was the same as that of the probe. The height and width of the model were 300 mm. In order to simulate water rising from the bottom of the container to progressively fill the container, the water surface elevated 20 mm every time and the dielectric constant was calculated. The medium parameters, boundary conditions, and probe voltage were set to be the same as in model N1. The model mesh, computation, and analysis method were also the same as in model N1. Water 2020, 12, x 8 of 14  Figure 6 shows the model N5 for simulation of calibration of vertical probe in terms of the report by [19]. The probe had a diameter of 4.6 mm and an outer separation of 70 mm. The length of the model was 285 mm, which was the same as that of the probe. The height and width of the model were 300 mm. In order to simulate water rising from the bottom of the container to progressively fill the container, the water surface elevated 20 mm every time and the dielectric constant was calculated. The medium parameters, boundary conditions, and probe voltage were set to be the same as in model N1. The model mesh, computation, and analysis method were also the same as in model N1.    Figure 6 shows the model N5 for simulation of calibration of vertical probe in terms of the report by [19]. The probe had a diameter of 4.6 mm and an outer separation of 70 mm. The length of the model was 285 mm, which was the same as that of the probe. The height and width of the model were 300 mm. In order to simulate water rising from the bottom of the container to progressively fill the container, the water surface elevated 20 mm every time and the dielectric constant was calculated. The medium parameters, boundary conditions, and probe voltage were set to be the same as in model N1. The model mesh, computation, and analysis method were also the same as in model N1.

Shaft Probes
After building the models using the probe parameters and computation using a finite element technique, the three-and two-dimensional weighting factors were obtained and the apparent dielectric Water 2020, 12, 545 9 of 14 constant was computed according to the theories mentioned above. Then, the results and discussions of the shaft probe were obtained. Figure 7 presents a comparison of the measurement sensitivity results of for the two types of TDR shaft probes determined using numerical and experimental investigations. Note that NS_A and ER_A represent data obtained from numerical simulation and experimental results, respectively using the type A probe reported by [12]. The experimental data were obtained by measuring a mixture of ethanol and deionized water set as the targeted medium. Furthermore, 2D_A and 3D_A represent data separately obtained from two-and three-dimensional numerical simulations using the type A probe in this work. It was observed that the effective dielectric constants obtained by the Comsol model were close to those derived by the seepage analysis software reported by [12]. This indicates that the numerical simulation technique performed by Comsol software was feasible and reliable. In addition, the effective dielectric constant obtained by the 3D Comsol model was slightly larger than that obtained by the 2D numerical results with the same software. In Table 1, the linear regression of measurement sensitivity equation is expressed and all the determination coefficients are 1, which indicates that for a given shaft dielectric constant the effective dielectric constant is a ideally linear function of the dielectric constant of the target medium. The 2D numerical results indicated that the measurement sensitivity values for the targeted medium for the A, B probes were 0.2669 and 0.4710, respectively, which is almost consistent with the values of 0.2664 and 0.4760 obtained in [12]. Meanwhile, the measurement sensitivity values for the targeted medium for the 3D numerical results of the A, B probes were 0.2842 and 0.4815, respectively, slightly higher values than those found in the 2D results. Overall, the Comsol modeling results match the numerical solutions reported in [12] very well, which proves that the 3D weighting theory and the Comsol numerical analyses are suitable for computing measurement sensitivity of the shaft probes. More importantly, the root mean square error (RMSE) values between the experimental data of the probe A and data computed by 3D or 2D models were 0.329 and 0.526, respectively, and the RMSE values between the experimental data of the probe B and data computed by 3D or 2D models were 1.609 and 2.115, respectively, which indicates that data obtained by the 3D models were more representative of the experimental data than those derived by the 2D models. This may be because the 3D method can encompass the spatial effect, which does not occur in the 2D method. In this regard, the 3D method has an advantages over the 2D method. discussions of the shaft probe were obtained. Figure 7 presents a comparison of the measurement sensitivity results of for the two types of TDR shaft probes determined using numerical and experimental investigations. Note that NS_A and ER_A represent data obtained from numerical simulation and experimental results, respectively using the type A probe reported by [12]. The experimental data were obtained by measuring a mixture of ethanol and deionized water set as the targeted medium. Furthermore, 2D_A and 3D_A represent data separately obtained from two-and three-dimensional numerical simulations using the type A probe in this work. It was observed that the effective dielectric constants obtained by the Comsol model were close to those derived by the seepage analysis software reported by [12]. This indicates that the numerical simulation technique performed by Comsol software was feasible and reliable. In addition, the effective dielectric constant obtained by the 3D Comsol model was slightly larger than that obtained by the 2D numerical results with the same software. In Table 1, the linear regression of measurement sensitivity equation is expressed and all the determination coefficients are 1, which indicates that for a given shaft dielectric constant the effective dielectric constant is a ideally linear function of the dielectric constant of the target medium. The 2D numerical results indicated that the measurement sensitivity values for the targeted medium for the A, B probes were 0.2669 and 0.4710, respectively, which is almost consistent with the values of 0.2664 and 0.4760 obtained in [12]. Meanwhile, the measurement sensitivity values for the targeted medium for the 3D numerical results of the A, B probes were 0.2842 and 0.4815, respectively, slightly higher values than those found in the 2D results. Overall, the Comsol modeling results match the numerical solutions reported in [12] very well, which proves that the 3D weighting theory and the Comsol numerical analyses are suitable for computing measurement sensitivity of the shaft probes. More importantly, the root mean square error (RMSE) values between the experimental data of the probe A and data computed by 3D or 2D models were 0.329 and 0.526, respectively, and the RMSE values between the experimental data of the probe B and data computed by 3D or 2D models were 1.609 and 2.115, respectively, which indicates that data obtained by the 3D models were more representative of the experimental data than those derived by the 2D models. This may be because the 3D method can encompass the spatial effect, which does not occur in the 2D method. In this regard, the 3D method has an advantages over the 2D method. ER_A [12] NS_B [12] ER_B [12] 2D_A 3D_A 2D_B 3D_B Figure 7. Comparison of measurement sensitivity of the TDR probes obtained from numerical simulation (NS) and experimental results (ER). NS_A, ER_A, NS_B, and NS_B represent data obtained from numerical simulation and experimental results using type A and B probes, respectively, as reported by Zhan et al. (2015). Furthermore, 2D_A, 3D_A, 2D_B, and 3D_B represent data separately obtained from two-and three-dimensional numerical simulations using type A and B probes, respectively, in this work.

Sample Calibration
In this section, different sample and probe rod diameters were considered in the case of modeling. According to the simulation scheme and the new computation method, the results and discussions based on sample calibration could be obtained. Figure 8 shows the relationship between apparent dielectric constant and dielectric constant of the sample. It could be observed that for a given probe with a diameter of 0.55 mm, the apparent dielectric constant was closer to the dielectric constant of the sample as the sample diameter increased from 30 mm to 80 mm. It was also observed that the difference of the apparent dielectric constant became larger as the dielectric constant of the sample increased. In addition, for a given sample diameter (the sample diameter was 30 mm in this work) the apparent dielectric constant was closer to the dielectric constant of the sample as the rod diameter decreased from 4.55 mm to 0.55 mm. This means that the smaller the probe size or the larger the sample diameter, the smaller the difference in the apparent dielectric constant with respect to the dielectric constant of the sample. Moreover, it was also observed that the influence of the change of sample and rod diameter on the apparent dielectric constant increased when the dielectric constant of the sample increased. Nevertheless, there was a limited impact of the cylindrical sample and rod diameter on apparent dielectric constant.
Water 2020, 12, x 10 of 14 obtained from two-and three-dimensional numerical simulations using type A and B probes, respectively, in this work. In this section, different sample and probe rod diameters were considered in the case of modeling. According to the simulation scheme and the new computation method, the results and discussions based on sample calibration could be obtained. Figure 8 shows the relationship between apparent dielectric constant and dielectric constant of the sample. It could be observed that for a given probe with a diameter of 0.55 mm, the apparent dielectric constant was closer to the dielectric constant of the sample as the sample diameter increased from 30 mm to 80 mm. It was also observed that the difference of the apparent dielectric constant became larger as the dielectric constant of the sample increased. In addition, for a given sample diameter (the sample diameter was 30 mm in this work) the apparent dielectric constant was closer to the dielectric constant of the sample as the rod diameter decreased from 4.55 mm to 0.55 mm. This means that the smaller the probe size or the larger the sample diameter, the smaller the difference in the apparent dielectric constant with respect to the dielectric constant of the sample. Moreover, it was also observed that the influence of the change of sample and rod diameter on the apparent dielectric constant increased when the dielectric constant of the sample increased. Nevertheless, there was a limited impact of the cylindrical sample and rod diameter on apparent dielectric constant.

Water Calibration
In this section, different model boundary size, probe end effect, and probe orientation were considered in the process of modeling. According to the simulation scheme and the new computation method, the results and discussions based on water calibration could be obtained. Figure 9 shows contour surfaces of 90% and 50% sample volume. The 90% level is to represent the sample volume and the 50% sample volume is to demonstrate the distribution of probe sensitivity

Water Calibration
In this section, different model boundary size, probe end effect, and probe orientation were considered in the process of modeling. According to the simulation scheme and the new computation method, the results and discussions based on water calibration could be obtained. Figure 9 shows contour surfaces of 90% and 50% sample volume. The 90% level is to represent the sample volume and the 50% sample volume is to demonstrate the distribution of probe sensitivity within the sample volume. It was observed that sample volume in the cross section view was in accordance with that obtained by the Geostudio software commonly used in hydrology [18]. Meanwhile, it was observed that the shape of sample volume at the end of the probe was different from that at the main part of the probe and the distribution zone exceeded the length of the probe into the deeper medium, which obviously influenced the computational results. In order to determine this difference, the calculating zone only employed the zone that the probe located, as shown in Figure 9c. It was observed that the contour surfaces of the 90% and 50% sample volume were slightly influenced by the end of the probe. Moreover, the sharp end that exists in the practical application is involved in this work, as shown in Figure 9d. The shape of the sample volume progressively decreased as the rod of the probe became slimmer, which also had an impact on the computational dielectric constant.
Water 2020, 12, x 11 of 14 within the sample volume. It was observed that sample volume in the cross section view was in accordance with that obtained by the Geostudio software commonly used in hydrology [18]. Meanwhile, it was observed that the shape of sample volume at the end of the probe was different from that at the main part of the probe and the distribution zone exceeded the length of the probe into the deeper medium, which obviously influenced the computational results. In order to determine this difference, the calculating zone only employed the zone that the probe located, as shown in Figure 8c. It was observed that the contour surfaces of the 90% and 50% sample volume were slightly influenced by the end of the probe. Moreover, the sharp end that exists in the practical application is involved in this work, as shown in Figure 8d. The shape of the sample volume progressively decreased as the rod of the probe became slimmer, which also had an impact on the computational dielectric constant.  Figure 10 shows a comparison of computational and experimental results. Note that the water surface level was above the plane of the probe in Figure 10. The water surface level fell in N1, N2, N3, N4, and [18], whereas the water surface was elevated in [19]. It was observed that as the distance of water surface from the probe center plane was large, the apparent dielectric constant was high and was equal to the dielectric constant of water. As water surface fell, the apparent dielectric constant was simultaneously influenced by the water and air medium. As the dielectric constant of the air is much less than that of the water, the apparent dielectric constant decreases as the distance of the water surface from the probe center plane decreases. The decrease becomes obvious as the water surface approaches to the probe center plane because the probe has a limited influence scope. In Figure 10a, the apparent dielectric constant Ka obtained by model N1 is slightly different from that derived from the rest of models; this difference may be attributed to the dimension ratio of the different media and their locations relative to the probe. Meanwhile, the data obtained from model N2, N3, and N4 were almost identical, indicating that the end effect, including the sharp end and probe end electric field distribution, was slight and may be negligible. In addition, when comparing the computational apparent dielectric constant and the experimental ones at a given distance from the water surface to the plane where the center of the probe was located, the apparent dielectric constant derived from both methods had the same variation trend in terms of the distance from the water surface to the probe center. There was less of a difference when the distance was large and more of a difference as the distance became smaller. In other words, when the distance from the water  Figure 10 shows a comparison of computational and experimental results. Note that the water surface level was above the plane of the probe in Figure 10. The water surface level fell in N1, N2, N3, N4, and [18], whereas the water surface was elevated in [19]. It was observed that as the distance of water surface from the probe center plane was large, the apparent dielectric constant was high and was equal to the dielectric constant of water. As water surface fell, the apparent dielectric constant was simultaneously influenced by the water and air medium. As the dielectric constant of the air is much less than that of the water, the apparent dielectric constant decreases as the distance of the water surface from the probe center plane decreases. The decrease becomes obvious as the water surface approaches to the probe center plane because the probe has a limited influence scope. In Figure 10a, the apparent dielectric constant K a obtained by model N1 is slightly different from that derived from the rest of models; this difference may be attributed to the dimension ratio of the different media and their locations relative to the probe. Meanwhile, the data obtained from model N2, N3, and N4 were almost identical, indicating that the end effect, including the sharp end and probe end electric field distribution, was slight and may be negligible. In addition, when comparing the computational apparent dielectric constant and the experimental ones at a given distance from the water surface to the plane where the center of the probe was located, the apparent dielectric constant derived from both methods had the same variation trend in terms of the distance from the water surface to the probe center. There was less of a difference when the distance was large and more of a difference as the distance became smaller. In other words, when the distance from the water surface level to the plane where the center of the probe is located becomes smaller, the values of dielectric constant obtained from the experimental results reported by [19] are slightly greater than those computed by numerical simulation. As the distance from the water surface to the probe center plane becomes zero, half of the probe is immersed in the water, while the other half is in the air. In Figure 10b, it can be observed that the ratio of apparent dielectric constant K a to the initial apparent dielectric constant K 0 obtained by models in the initial water level state is almost identical for all of models, which means the model dimension has less impact on the K a /K 0 ratio. Moreover, there is also a difference as the distance of water surface to the plane where the center of the probe is located decreases, as shown in Figure 10a. This difference may be attributed to the size errors of the probe or boundary conditions or frequency characteristics, and further studies are needed.
Water 2020, 12, x 12 of 14 surface level to the plane where the center of the probe is located becomes smaller, the values of dielectric constant obtained from the experimental results reported by [19] are slightly greater than those computed by numerical simulation. As the distance from the water surface to the probe center plane becomes zero, half of the probe is immersed in the water, while the other half is in the air. In Figure 10b, it can be observed that the ratio of apparent dielectric constant Ka to the initial apparent dielectric constant K0 obtained by models in the initial water level state is almost identical for all of models, which means the model dimension has less impact on the Ka/K0 ratio. Moreover, there is also a difference as the distance of water surface to the plane where the center of the probe is located decreases, as shown in Figure 10a. This difference may be attributed to the size errors of the probe or boundary conditions or frequency characteristics, and further studies are needed.  Figure 11 shows the calculated relative dielectric permittivity by model N5 using the method proposed in this work and the measured values reported by [19] as a function of the distance of the air-water interface to the probe midpoint plane. Note that the water surface level is initially below the plane of the probe and progressively elevates until it is above the plane of the probe. It was observed that the three-rod vertical probes were sensitive to the advance of the water surface as the surface approached the middle and upper rods. There was a smaller difference in the computational and experimental results at the region further away the middle rod, while a greater difference was observed as the water surface approached the middle rod. With a 2D method, the authors of [19] also found that the three-rod vertical probe design resulted in poor quality waveforms when the interface was located between the central and upper rods. Figure 11. Calculated relative dielectric permittivity by model N5 and the measured values reported by [19] as a function of the distance of the air-water interface to the probe midpoint plane.  Figure 11 shows the calculated relative dielectric permittivity by model N5 using the method proposed in this work and the measured values reported by [19] as a function of the distance of the air-water interface to the probe midpoint plane. Note that the water surface level is initially below the plane of the probe and progressively elevates until it is above the plane of the probe. It was observed that the three-rod vertical probes were sensitive to the advance of the water surface as the surface approached the middle and upper rods. There was a smaller difference in the computational and experimental results at the region further away the middle rod, while a greater difference was observed as the water surface approached the middle rod. With a 2D method, the authors of [19] also found that the three-rod vertical probe design resulted in poor quality waveforms when the interface was located between the central and upper rods.
Water 2020, 12, x 12 of 14 surface level to the plane where the center of the probe is located becomes smaller, the values of dielectric constant obtained from the experimental results reported by [19] are slightly greater than those computed by numerical simulation. As the distance from the water surface to the probe center plane becomes zero, half of the probe is immersed in the water, while the other half is in the air. In Figure 10b, it can be observed that the ratio of apparent dielectric constant Ka to the initial apparent dielectric constant K0 obtained by models in the initial water level state is almost identical for all of models, which means the model dimension has less impact on the Ka/K0 ratio. Moreover, there is also a difference as the distance of water surface to the plane where the center of the probe is located decreases, as shown in Figure 10a. This difference may be attributed to the size errors of the probe or boundary conditions or frequency characteristics, and further studies are needed.  Figure 11 shows the calculated relative dielectric permittivity by model N5 using the method proposed in this work and the measured values reported by [19] as a function of the distance of the air-water interface to the probe midpoint plane. Note that the water surface level is initially below the plane of the probe and progressively elevates until it is above the plane of the probe. It was observed that the three-rod vertical probes were sensitive to the advance of the water surface as the surface approached the middle and upper rods. There was a smaller difference in the computational and experimental results at the region further away the middle rod, while a greater difference was observed as the water surface approached the middle rod. With a 2D method, the authors of [19] also found that the three-rod vertical probe design resulted in poor quality waveforms when the interface was located between the central and upper rods. Figure 11. Calculated relative dielectric permittivity by model N5 and the measured values reported by [19] as a function of the distance of the air-water interface to the probe midpoint plane. Figure 11. Calculated relative dielectric permittivity by model N5 and the measured values reported by [19] as a function of the distance of the air-water interface to the probe midpoint plane.

Conclusions
The TDR technique is an effective method to measure the volume of soil water content. Many types of TDR probes have been developed and can be used according to the condition of soils, and numerical simulation has been used to design and optimize the probes. However, the complex probe pattern and surrounding conditions limit the application of the two-dimensional weighting factor theory. Therefore, a new method including a three-dimensional weighting theory and simulation technique by the finite element Comsol software based on Laplace's equation was proposed and an extensive verification performed with shaft probes and rod probes in the present work. The results of the simulation were compared with those of the experiments, and conclusions were reached as follows.
For the shaft probes, the combination of the three-dimensional weighting theory and numerical simulation based on Laplace's equation could be used to calculate the spatial sensitivity of the probes well through comparing the experimental and calculated results. The values of the dielectric constant obtained by the three-dimensional method were slightly larger than those derived by two-dimensional method, which indicates that the 3D method encompasses more spatial factors. In addition, the 3D simulation values were more representative of the experimental data, which presents an advantage over the 2D method for the shaft probe.
For rod probe computation for samples, the smaller the probe size or the larger the sample diameter, the smaller the difference between apparent dielectric constant and the dielectric constant of the sample, but there was a limited impact of cylindrical sample shape on the apparent dielectric constant. In addition, for rod probe computations for water, as the distance between the water surface and the probe center plane increased, the difference between the computational and experimental results decreased. There was a smaller difference when the water surface was located at the probe center for the horizonal probe and the vertical probe. However, the variation trend of the computational and experimental data was identical. Moreover, there was less of an impact of the end effect including the sharp end and probe end electric field on the calculation of the dielectric constant and the ratio of the K a to K 0 , while model size had a certain effect on the dielectric constant but not on the ratio of them. The new method can provide a basis to preliminarily design the probe, especially for complex spatial shape of probes and samples.