New Method for Modelling Seasonal Variation in Resistance and Performance of Earthing Systems

Abstract

The current techniques utilized for estimating seasonal fluctuations in earthing system resistance, including artificial neural networks (ANNs) and correlation/correction factors, rely on resistance records, soil resistivity measurements, and meteorological data collected across broad areas. However, they frequently fail to consider the impact of soil conditions and properties at the actual earthing location. As a solution, this research introduces a new method that models atmospheric conditions as soil suction and incorporates hydraulic soil properties (soil water retention characteristics and hydraulic conductivity) to estimate the seasonal changes in earthing resistance and performance. To illustrate this approach, this study constructs geometric models of vertical earthing rods for three homogeneous soil textures (clayey, silty, and sandy) utilizing COMSOL Multiphysics software. By coupling the differential equations governing electric current and water flow using Archie’s formula and solving numerically with the finite element method (FEM) for various soil suctions, this research reveals that soil water retention and resistivity variations are notably influenced by soil texture. Sandy soil displays higher variability, silt soil demonstrates moderate changes, while clayey soil exhibits lower fluctuation. By linking soil resistivity changes to soil suction and hydraulic properties, this innovative method predicts seasonal trends in earthing resistance and performance.

1. Introduction

The effectiveness of earthing electrodes in dispersing fault currents and surges is contingent upon multiple factors such as electrode type, size, composition, and the resistivity of the soil surrounding the electrode [1,2,3,4,5]. To ensure that power system components function reliably, it is imperative that the resistance of earthing electrodes within the ground remains low and within the specified limits for the application. Earthing resistance is determined by the resistivity of the soil and geometric factors (shape and size) of the earth electrode. Soil resistivity is the ability of soil to impede electric current flow. Because soil resistivity depends on salinity, compaction/porosity, temperature, and water content [1,2,4], which are factors susceptible to atmospheric conditions, soil resistivity is subject to seasonal variability. Sometimes, resistivity variation attains high values that can inhibit the reliable operation of earthing systems [6,7,8,9]. To limit the impact of high seasonal resistivity, the IEEE Std. 81 recommends regularly scheduled inspections of earthing installations to identify sites with high resistances to implement corrective measures [10]. Due to inadequate budget and manpower for planning and executing inspection and resistance measurement of thousands of earthing sites, most electric utility companies cannot perform periodic earthing inspections. This means, that earthing sites with deteriorating conditions may be undetected in preventing further deterioration that can jeopardize the performance of the earthing system.

The current methods for estimating the seasonal variation in earthing resistance use artificial neural networks (ANNs) or correlation/correction factors that are formulated with meteorological data and records of earthing resistance and/or soil resistivity measurements [1,11,12]. By associating soil resistivity with meteorological data, these methods ignore the influence of soil properties and site features on the seasonal variation in soil resistivity. These properties and features regulate the rate and extent of water loss from the subsoil to the atmosphere. For instance, while the atmosphere determines the amount of rainfall, it is the topography of the ground, a site feature, that controls water runoff and the rate of infiltration of water in the soil, while the quantity of water retained in the soil is determined by the soil property. Therefore, incorporating site features and soil properties in the modelling of seasonal variation in soil resistivity should enhance the accuracy of the prediction of seasonal variation in earthing resistance and performance.

The model of soil resistivity currently used in the computational analysis of earthing systems is generic. Such models lack the capacity for evaluating the seasonal variation in earthing resistance as resistance values are fixed and not dynamic. Formulating soil resistivity models as functions of soil properties makes the evaluation and prediction of the variability in soil resistivity and resistance reduction capacity of earthing enhancement materials (EEMs) in different soil conditions, seasons, and climates possible. The effect of atmospheric and environmental conditions on the energy state of soil water determines the quantity of water flowing through, retained in, or released from the soil matrix. This implies that hydraulic conductivity, water retention characteristics, and diffusivity of soil are the properties of soil that indicate the response of the soil to these external influences. Under certain assumptions, these influences can be modelled in the form of soil suction or pressure head.

Consequently, a new method for modelling soil resistivity as a function of hydraulic conductivity and water retention characteristics is proposed. It requires the multiphysics coupling of the differential equation for hydraulic and electric flows in the ground and is solved numerically for different values of soil suction or pressure head.

1.1. Theoretical Fundamentals of the Proposed Method

Fundamental descriptions from the domain of soil physics are required in developing the proposed multiphysics earthing prediction model, subsequently called the hydro-geoelectric model of earthing systems. These descriptions are requisite for understanding the hydraulic properties of soil and their roles in the seasonal variability in soil resistivity in response to external influences. This subsection discusses fundamental concepts such as soil composition and texture, water content of the soil, forms of energy state of the soil, and soil water retention characteristics. It concludes with the modelling of external influences such as relative humidity and its correlation to the energy state of soil water.

1.1.1. Soil Composition and Texture Classification

The subsoil or ground consists of stratified layers, each containing solid parent materials of different physical, morphological, and chemical properties [13]. These amorphous solids also called particles or grains have different sizes with spaces (voids or pores) between them. Depending on the state of the soil, the voids are filled with either air (gas), water (liquid), or both. The distribution of particle sizes determines the size and number of pores, which determine the water retention and release capacity of the soil under external influences. Soil grains of fine sizes such as clay and silt have numerous small-sized pores that can retain more water for longer periods unlike coarse aggregates (sand) with fewer large-sized pores that quickly drain and retain little water. Depending on size, the grains of soils are classified as clay (<0.002 mm), silt (0.002 mm–0.05 mm), and sand (0.05 mm–2 mm). The percentage of clay, silt, and sand in soil samples determines its textural classification as clay, silty clay, sandy clay, silty clay loam, clay loam, sandy clay loam, silt, silt loam, loam, sandy loam, loamy sand, or sand.

1.1.2. Soil Water Content

From the description of soil composition, a typical volume of soil consists of volumes of amorphous solids, liquid (water), and air (gas) within the spaces (voids). The total volume of soil, $V_t$, is the sum of the volumes of solid aggregate, $V_s$, the volume of liquid or water, $V_w,$ and the volume of air, $V_a$. The volume of void, $V_v$, is the sum of $V_s$ and $V_s$. The amount of water in the soil or soil water content can be described by the volumetric water content, $\mathsf{\theta }$, degree of saturation, $S_d$, or porosity $\varphi$.

Volumetric water content is the volume fraction of soil that is water. Porosity represents the ratio of the volume of void to the volume of soil, while the degree of saturation of soil is the ratio of the volume of water to the volume of void. Mathematically, $\mathsf{\theta }$, $S_d$, and $\varphi$ are expressed by Equations (1)–(3).

$$\mathsf{\theta }=V_w/V_t$$

(1)

$$\varphi =V_V/V_t$$

(2)

$$S_d=V_W/V_V=\mathsf{\theta }/\varphi$$

(3)

Depending on the amount of water within the soil pores, the degree of saturation of soil can be classified as saturated, partially saturated, or unsaturated. The pores of saturated soils are filled with water, whereas air fills the pores of unsaturated soil. Partial saturation is when air and water coexist within the soil pores. The value of the degree of saturation of soil lies between 0 and 1, which are the extreme values representing saturated and unsaturated soils.

1.1.3. Energy State of Soil Water

Under hydro-static conditions, the movement of water in the soil is entirely due to mechanical potentials. The interdependence shared between the energy potential of soil water and water content ensures that changes in one produce changes in the other. The energy state of soil water is described by the soil water potential (SWP), which is the sum of the osmotic, electric, chemical, gravitational, pressure, and matric potentials. In groundwater flow, the energy state of soil water is expressed in the form of a hydraulic head and defined as SWP per unit weight of water. It is measured in the unit of length and expressed in Equation (4), as the sum of pressure head, $h_P=P/\mathsf{\gamma }$, and elevation head, $h_z$.

$$H=P/\mathsf{\gamma }+h_z$$

(4)

where $\mathsf{\gamma }=\mathrm{g}{\rho }_w$ is the unit weight of water, $P$ [Pa] is the pressure due to the weight of the water column above the point, $\mathrm{g}$ [ms⁻²] is the acceleration due to gravity, and ${\rho }_w$ is the density of water [kgm⁻³].

Using the ground surface as a reference datum and the soil region above the water table as unsaturated, the elevation head component can be ignored such that $H\approx h_P$.

The energy state of soil water can be described by soil suction. Also called matric or capillary suction, ${\psi }_s$, represents the strength of attraction of water molecules to a medium, which allows water to be held tightly within the soil pore under high tension or negative pressure [14]. ${\psi }_s$ is measured in the unit of pressure (Pa). Although describing different parameters, ${\psi }_s$ and $h_P$ are related as in Equation (5).

$$h_P={-\psi }_s/\mathrm{g}{\rho }_w$$

(5)

1.1.4. Relative Humidity and Soil Suction

The ground surface is the interface for the exchanges between the atmosphere and the ground. The cumulative effect of these atmospheric conditions on the subsoil can be described by soil suction at the ground surface in terms of relative humidity, $H_R$ [15,16]. The magnitude of ground surface suction depends on atmospheric pressure and humidity, which are temperature-dependent. Relative humidity represents the measure of the amount of water vapour in the air relative to the saturation capacity of air at a specific pressure and temperature.

Equation (6) is the Magnus–Teten’s equation expressed in terms of atmospheric temperature, $T$, and dew point temperature, $T_d$, which are both measured in °C and are used to determine $H_R$ in percentage.

$$H_R=100\exp \left(\frac{17.625T_d}{243.04+T_d}-\frac{17.625T}{243.04+T}\right)$$

(6)

Equation (7) [16] relates the water potential of the soil in Mega-Pascal (MPa) to the relative humidity, $H_R$.

$$\psi ={\psi }_0-(RT/M)\ln \left(0.01H_R\right)$$

(7)

where $R$ is the universal gas constant (8.314 J·K⁻¹·mol⁻¹), $T$ is the temperature in K, $M$ is the molecular weight of water (0.01802 kg/mol), and ${\psi }_0$ is the potential of pure water under standard temperature and pressure (STP), 0 MPa. Equation (7) can be used to evaluate the suction within the soil under equilibrium conditions [15].

To illustrate how Equation (7) is applied, on a particular day the atmospheric and dew point temperatures were 7 °C and 0 °C for Johannesburg and 37 °C and 17 °C for Jeddah, respectively. Using Equation (6), the relative humidity is computed as 61.05% for Johannesburg and 30.83% for Jeddah. Consequently, if the ground surface and air are at equilibrium, then using Equation (7), the magnitude of water potential at the ground surface is 63.78 kPa in Johannesburg and 168.37 kPa in Jeddah.

1.1.5. Soil Water Retention Characteristics

Soils of different textures have distinct water retention and release because they contain different percentages of clay, silt, and sand. As a parameter that describes the distribution of soil grain size, soil texture has been shown to influence the characteristics of soil [17]. Consequently, since layers of the subsoil have different soil textures, water content is different in each layer and discontinuous across adjacent soil layers [18].

The soil water characteristics curve (SWCC) is the curve used to describe the water retention and release attribute of the soil [19,20,21]. It is a graphical relationship between the soil volumetric water content, $\mathsf{\theta }$, (or degree of saturation) and the soil water potential, $\psi$, which is represented by soil suction or pressure head. Due to hysteresis behaviour between the saturation and desaturation cycles, SWCC typically has two curves for the saturation and desaturation processes. Desaturation processes impact soil resistivity more than saturation; therefore, rather than the SWCC, the soil water retention curve (SWRC) is used to describe water retention during water loss.

Figure 1 is a typical SWRC showing the critical regions of saturation, capillary, and residual states of the soil. Saturation and residual regions represent the soil states where water and air predominantly fill the soil pores. The capillary region is a transitional region, where in the presence of air, water is held by capillary forces. Within this mid-region, a slight change in soil suction produces a substantial change in the volumetric water content and resistivity of soil.

Effective saturation, $S_e$, is sometimes preferred to the degree of saturation, $S_d$, in describing the hydraulic state of a soil medium. $S_e$ is similar to the normalized degree of saturation and expressed in Equation (8) as a function of volumetric water content, $\mathsf{\theta }$, residual volumetric water content and saturated volumetric water content of the medium, ${\mathsf{\theta }}_r$ and ${\mathsf{\theta }}_s$.

$$S_e=\frac{\mathsf{\theta }-{\mathsf{\theta }}_r}{{\mathsf{\theta }}_s-{\mathsf{\theta }}_r}$$

(8)

Constitutive models proposed by Brooks & Corey [22], van Genuchten [23], and Fredlund & Xing [24], etc., are used to fit the SWRC. These models use fitting parameters to mathematically model the effective degree of saturation, $S_e$, or volumetric water content, $\mathsf{\theta }$, in terms of soil suction or pressure head, $h_p$. The Mualem–Genuchten (MG) model, a modified version of the van Genuchten model, describes the SWRC as follows:

At saturation, $h_p\ge 0;$ $S_e\left(h_p\right)=1,$ and $\mathsf{\theta }\left(h_p\right)={\mathsf{\theta }}_s$.

During unsaturation, $h_p<0;$ then, $S_e\left(h_p\right)$ and $\mathsf{\theta }\left(h_p\right)$ are given in Equations (9) and (10) as

$$S_e\left(h_p\right)={\left(\frac{1}{1+{\left({\alpha }_h\left|h_p\right|\right)}^n}\right)}^m$$

(9)

$$\mathsf{\theta }\left(h_p\right)={\mathsf{\theta }}_r+({\mathsf{\theta }}_s-{\mathsf{\theta }}_r){\left(\frac{1}{1+{\left({\alpha }_h\left|h_p\right|\right)}^n}\right)}^m$$

(10)

where ${\alpha }_h$ [m⁻¹] is a fitting parameter equal to the inverse of the bubbling pressure or suction/pressure head at air entry and $n$ is a dimensionless fitting parameter related to the particle size distribution of the medium. At $n>1$, the relationship between $n$ and $m$ is $m=1-1/n$ [25]. Suction or pressure head at air entry is a unique parameter of the medium that can be read off from the SWRC and represents the pressure where the largest pore begins to drain [26].

1.1.6. Soil Resistivity Characteristics Curve

The relationship between the pairs of soil resistivity, suction, and water content or degree of saturation has been investigated and established in the literature [27,28]. Laboratory investigation on residual soil observed that the electrical resistivity–water characteristic curve (RWCC) and soil–water characteristic curve (SWCC) have similar shapes and have a nonlinear relationship with the volumetric water content [28]. Like SWRC, the RWCC of the soil can be described by the van Genuchten model using appropriate fitting parameters. The linear proportionality between the fitting parameters enables the prediction of the RWCC from the SWCC.

The tripartite relationship between water content, soil suction, and resistivity of soil is gaining profound usefulness in engineering applications. For instance, doctoral-level research utilized it to model, analyse, and predict corrosion processes for underground metallic pipes subjected to differential aeration [29].

This paper extends the application of the tripartite relationship for evaluating and predicting the resistance of earthing systems under variable saturation of the ground. It is demonstrated by multiphysics simulation modelling. Atmospheric conditions are modelled as soil suction and are used as variables to estimate water flow and retention. A form of Archie’s equation is subsequently used to evaluate the electrical resistivity of soil corresponding to the soil suction and retained water content. Earthing resistance is evaluated from post-numerical analysis of electric current conservation in the soil. This paper proposes and demonstrates a simulation method for estimating the variation in earthing resistance due to changes in atmospheric conditions.

1.2. Formulation of the Hydro-Geoelectric Earthing Model

Although atmospheric conditions regulate the availability of water to the soil and the rate, frequency, and direction of water movement within the subsoil [20], it is soil texture, and the thickness and depth of each layer relative to the ground surface and water table, that determines the quantity and energy state of soil water in the subsoil. Since climate and environmental conditions control the mechanism of water retention and release in the subsoil, they describe and regulate the energy state of soil water [7,30].

As water flows in the soil into or out of the soil, the soil water content influences the electrical conductivity of the soil. The water content of soil affects soil resistivity by influencing the mechanisms of electrical conduction in the soil, which is either electrolytic or surface conductance [4]. Soil water decreases and increases soil salinity, which promotes surface conductance, while the increase in soil water content decreases salinity and promotes electrolytic conduction. A significant decrease in soil water content produces high resistivity values of the soil [28].

To describe electric current flow in variably saturated soil, the relationship between electrical conductivity and water content of the soil under changing soil conditions is required. The following subsection contains a discussion of water flow and electric potential distribution in the soil.

1.2.1. Water Flow in Porous Medium

The flow of water in a porous medium is fundamentally described by Darcy’s law and expressed mathematically in Equation (11). Darcy’s law states that the flux rate of liquid through a saturated porous medium is proportional to the gradient of the hydraulic head, $\nabla H$ [20].

$$q=-K_s∆H/∆L$$

(11)

where $q$ [L3L-2T-1] is the volumetric flow rate, $K_s$ [LT-1] is the saturated hydraulic conductivity and a characteristic of the medium, $∆H$ [L] is the difference in the total hydraulic head between two points in the saturated porous medium separated by a distance of $∆d$ [L], then, $∆H/∆d=\nabla H$. The negative sign indicates that water flows in the direction of decreasing hydraulic head.

The vadose zone (i.e., the unsaturated region of soil that lies between the ground surface and water table) undergoes variable degrees of saturation. It is saturated during and immediately after heavy rainfall and becomes unsaturated hours to days after. As the soil desaturates, water is displaced by air in the soil pores, which creates the discontinuity of water molecules in the soil. This discontinuity increases the drag between the solid and fluid phases, which prevents the free flow of water through the soil. Consequently, hydraulic conductivity is considerably higher at saturation than during desaturation.

Soil texture also influences the hydraulic conductivity of soil. Coarse grain soils have hydraulic conductivity of several magnitudes higher than fine-grained soils at saturation. However, due to their large pores that quickly drain water, coarse soils have significantly smaller values of hydraulic conductivity under desaturation conditions than fine-textured soils.

Whereas Equation (11) describes water flow for saturated soil, Buckingham suggested that by replacing $K_s$ with $K$, the hydraulic conductivity function (HCF), Darcy’s law is modified and becomes suitable flow describing water flow in unsaturated porous medium [31]. HCF is expressed either as a function of water content, $\mathsf{\theta }$, as $K\left(\mathsf{\theta }\right)$ or by pressure head, $h_p$, as $K\left(h_p\right)$ [31,32]. Therefore, the Darcy–Buckingham description of unsaturated water flow in the soil is expressed in Equation (12) as

$$q=-K\left(h_p\right)\nabla \mathrm{H}$$

(12)

Applying the conservation of mass ensures that the amount of water entering or leaving a volume of the porous medium is consistent with the rate of change in volumetric water content. This is expressed mathematically by the Continuity equation or Equation (13).

$$\frac{\partial \mathsf{\theta }}{\partial t}=-\nabla ·q$$

(13)

where $\mathsf{\theta }$ [L3L-3] is the volumetric water content, $t$ [T] is time, and $q_x$, $q_y$ and $q_z$ are the components of volumetric flow rate in the x-, y-, and z-axis, respectively.

Equation (14) is a partial differential known as Richards’ equation, which is formed by substituting $q$ of Darcy–Buckingham law for $q$ in the Continuity equation.

$$\frac{\partial \mathsf{\theta }}{\partial t}=\nabla ·\left[K\left(h_p\right)\nabla H\right]$$

(14)

If lateral flows in the x and y directions are insignificant or ignored, then a 1D form of Richards’ equation, representing the vertical flow along the z-axis (negative downward) can be derived. Equation (15) is the 1D partial differential equation (PDE) expressed as a function of $h_p$ and $h_z$, given that $H$ is the sum of $h_p$ and $h_z$.

$$\frac{\partial \mathsf{\theta }}{\partial t}=\frac{\partial }{\partial z}\left[K\left(h_p\right)\left(\frac{\partial h_p}{\partial z}+1\right)\right]$$

(15)

Equation (15) describes the vertical flow of water in a porous medium at different degrees of saturation. The first and second terms account for the effects of capillarity gravitational flux, respectively.

Equation (16) is formed by expressing the left-hand side of Equation (15) in terms of $C_m\left(h_p\right)$ [L-1], the specific water capacity of the soil.

$$C_m\left(h_p\right)\frac{\partial h_p}{\partial t}-\frac{\partial }{\partial z}\left[K\left(h_p\right)\left(\frac{\partial h_p}{\partial z}+1\right)\right]=0$$

(16)

Equation (17) is a representation of Richards’ equation in terms of the specific storage, $S_s$ [L-1], and $S_e$ of the soil. This form is equivalent to Equation (18), which is the form of Richards’ equation in COMSOL Multiphysics software that gives an average error of less than 10% between the observed and simulated data [33,34].

$$\left[S_s{S}_e\left(h_p\right)+C\left(h_p\right)\right]\frac{\partial h_p}{\partial t}-\frac{\partial }{\partial z}\left[K\left(h_p\right)\left(\frac{\partial h_p}{\partial z}+1\right)\right]=0$$

(17)

$${\rho }_w\left(\frac{C_m}{{\rho }_{w}g}+S_{e}S\right)\frac{\partial P}{\partial t}+\nabla ·{\rho }_w\left(-\frac{K_s{K}_r}{\mu }\left(\nabla P+{\rho }_{w}g\nabla z\right)\right)=Q_m$$

(18)

where $K_s$ [LT-1] and $K_r$[-] are saturated and relative hydraulic conductivity, respectively, μ is fluid dynamic viscosity, $K$ [LT-1] is the hydraulic conductivity, which relates to $K_s$ as $K_s=(g{\rho }_w/\mu )K$, $z$ is elevation, and $Q_m$ [ML-2 T-1] represents the fluid source term or mass deposit in/out of the porous medium.

The values of $C_m$, $S_e$, and $K$ are determined as functions of $h_p$ and can be evaluated using the water retention curve. For instance, $C_m$ is the first derivative of the water retention curve,

$$C_m\left(h_p\right)=d\mathsf{\theta }/d{h}_p$$

(19)

$K$ is determined from the values of $K_s$ and $K_r$, since $K=K_r{K}_s$. $K_s$ can be evaluated in the laboratory from the falling or constant head tests of fine and coarse-textured soils, respectively. Alternatively, $K_s$ [m/day] can also be estimated using Equation (20) in terms of soil porosity $\varphi$ and ${\alpha }_h$ [35].

$$K_s=4.65\times {10}^{-2}\varphi {\left({\alpha }_h\right)}^2$$

(20)

The value of $K_r$ can be evaluated using Equation (21) [23,36].

$$K_r\left(h_p\right)=\frac{{\left[1-{\left({\alpha }_h{h}_p\right)}^{n-1}{\left[1+{\left({\alpha }_h{h}_p\right)}^n\right]}^{1/m}\right]}^2}{{\left[1-{\left({\alpha }_h{h}_p\right)}^n\right]}^{\frac{m}{2}}}$$

(21)

1.2.2. Electric Current Flow in Conductive Medium

The small finite size of the earthing electrode relative to the infinite soil makes modelling the electrode as a current point source suitable. It is the electrical conductivity, $\sigma$, property of the soil (inverse of resistivity) that permits the flow of electric current. The dispersion of current $I_f$ from the surface of the earthing electrode into the soil is governed by Ohms’ law. Equation ($22$) is a form of Ohms’ law expressed as a function of current density, $\mathrm{J}$ (A/m²), electric field intensity, $E$ (V/m), and electrical conductivity, $\mathsf{\sigma }$ (S/m).

$$\mathrm{J}=\mathsf{\sigma }E$$

(22)

The gradient operator relates electric potential to the electric field as in Equation (22). Equation ($23$) is used to determine the electric potential $V_p$ at any point due to the dispersed current

$$E=-\nabla V_p$$

(23)

Substituting for $E$ in Equation $\left(23\right)$, applying divergence on both sides of the resulting equation and simplifying gives Equation ($24$).

$$\nabla ·\mathrm{J}=-\sigma {\nabla }^2{V}_p$$

(24)

Equation ($24$) is the differential equation that relates the Laplacian of the electric potential to the divergence of current density, which is solved by using appropriate initial and boundary conditions.

1.2.3. Electrical Conductivity and Water Content

Archie’s empirical equation provides the relationship between the degree of saturation (water content) and the electrical conductivity of the porous medium. Equation (25) is Archie’s for variably unsaturated medium expressed in terms of the surface conductivity, ${\sigma }_s$, of fine-grained soils, and the electric conductivity of the liquid, ${\sigma }_w,$ in the porous medium [37].

$$\sigma ={\sigma }_w{\varphi }^{\mathsf{ɱ}}{\left(S_d\right)}^{\beta }+{\sigma }_s$$

(25)

where $\mathsf{ɱ}$ and $\beta$ are the cementation factor and saturation index respectively and $\beta \gg \mathsf{ɱ}$.

1.2.4. Coupling of Electric and Water Flows in Conductive Porous Medium

The following assumptions were adopted in developing the hydro-geoelectric model of earthing terminals in the ground.

(1)

Water in the medium is treated as incompressible.

(2)

No dynamic change in liquid mass in the soil volume at any test pressure head, ${∆Q}_m=0$.

(3)

Soil electrical conductivity is isotropic and determined by the volumetric water content function of the pressure head.

(4)

Surface electrical conductivity of soil is neglected.

(5)

Ground surface suction is equilibrated within the soil domain.

(6)

Steady-state condition; $\partial h_p/\partial t=0$.

If ${\sigma }_s$ is neglected, and the product term ${\sigma }_w{\varphi }^{\mathsf{ɱ}}$ (a property of the porous medium) is assigned as $A_{\sigma }$, then Equation (25) takes the form of Equation $\left(26\right)$.

$$\sigma \left(S_d\right)=A_{\sigma }{\left(S_d\right)}^{\beta }$$

(26)

At saturation, the volumes of void and water are equal; therefore, $\varphi ={\mathsf{\theta }}_s$. Equation (26) can then be expressed as a function of $\mathsf{\theta }$ in Equation (27).

$$\sigma \left(\mathsf{\theta }\right)=A_{\sigma }{\left(\frac{\mathsf{\theta }}{{\mathsf{\theta }}_s}\right)}^{\beta }$$

(27)

To account for the actual saturation or desaturation, $S_e$ is used rather than $S_d$, such that Equation (26) becomes

$$\sigma \left(S_e\right)=A_{\sigma }{\left(S_e\right)}^{\beta }$$

(28)

Substituting for $S_e$ using Equation (8) or Equation (9) expresses Equation (28) as Equation (29) as a function of $\mathsf{\theta }$ or as Equation (30) as a function of $h_p.$

$$\sigma (S_e)=A_{\sigma }{\left(\frac{\mathsf{\theta }-{\mathsf{\theta }}_r}{{\mathsf{\theta }}_s-{\mathsf{\theta }}_r}\right)}^{\beta }$$

(29)

$$\sigma (h_p)=A_{\sigma }{\left[\frac{1}{1+{\left(\alpha \left|h_p\right|\right)}^n}\right]}^{m·\beta }$$

(30)

Consequently, substituting for $\sigma$ in Equation (24) yields Equations (31) and (32), respectively.

$$\nabla ·J=-\left[A_{\sigma }{\left(\frac{\mathsf{\theta }-{\mathsf{\theta }}_r}{{\mathsf{\theta }}_s-{\mathsf{\theta }}_r}\right)}^{\beta }\right]·\left({\nabla }^2{V}_p\right)$$

(31)

$$\nabla ·J=-\left\{A_{\sigma }{\left[\frac{1}{1+{\left(\alpha \left|h_p\right|\right)}^n}\right]}^{m·\beta }\right\}·\left({\nabla }^2{V}_p\right)$$

(32)

Equations (31) and (32) are the coupled equations that describe the conservation of water and electric current flow in variably saturated conductive porous medium expressed in terms of water content and pressure head, respectively.

2. Materials and Methods

2.1. Methodology

To demonstrate the capacity of the proposed hydro-geoelectric method in predicting the resistance and performance of earthing systems under variable conditions of soil, a simulation model of vertical earth rods in homogeneous soil was built in COMSOL Multiphysics software. Within the model builder, 1D–3D geometric models of an object can be created, discretized/meshed, and parameters and numerical solvers assigned for solving the governing equations over the discretized space. The software also enables probes and post-simulation analysis that may be presented in tables and annotated plots.

Actual dimensions and properties of solid cylindrical earth rods were used, while the hydraulic and electrical parameters of the three homogeneous soils (sandy, clayey, and silty soils) were sourced from references [35,38,39]. Using appropriate initial and boundary conditions, the partial differential equations governing the physics of electric and water flows were solved numerically using the finite element method (FEM) for different degrees of soil saturation. The influence of atmospheric conditions on the soil was modelled as soil suction on the ground surface. For the range of soil suctions, the soil water content or degree of saturation and corresponding resistivity of soil are determined at equilibrated suction values. Earthing resistance is subsequently computed from post-simulation analysis using the electric potentials at any soil boundary relative to the electric potential at the surface of the electrode. To validate the simulation model, resistance evaluated from post-simulation analysis is compared to resistance computed using existing analytical formulas in the literature.

2.1.1. Geometric Modelling

The geometric model of the vertical earth rod in homogeneous soil is built using the 2D axisymmetric space of the Model builder of the software. Leveraging on the vertical radial symmetry of cylindrical rod and ground enables the geometric model of the earth rod to be built with the 2D axisymmetric space of the software [40]. Unlike the 3D model equivalent, the 2D results in reduced meshing and computational costs. For a cylindrical rod of length 3 m and radius 9.525 mm, the width and depth of the homogeneous soil domain are 25 m and 20 m, respectively, as shown in Figure 2. The model shows a critical space around the earth electrode called the region of influence (ROI) of the earth electrode. Also referred to as the sphere of influence, the ROI is the region around the electrode where the energy dissipated from the earth electrode is significant. The volume of the ROI depends on the physical shape and size of the earthing electrode [41]. For a vertical earth rod of length $L$, the ROI is bounded by the space 2.5 L from the surface of the rod [42].

2.1.2. Interface of Flow Physics in COMSOL Multiphysics Software

The equations governing the flow of electric current and fluid in a conductive porous medium are accessed in the Electric Currents (EC) interface and Richards’ Equation (RE) interface, respectively. EC interface is hosted in the Alternating Current and Direct Current (AC/DC) module under the Electric Fields and Currents submodule, while RE interface is hosted with the Fluid Flow module under the Porous Media and Subsurface Flow submodule. The EC interface uses electric potential, $V_p$, as the dependent variable to solve the current conservation based on Ohm’s law [34]. The RE interface uses pressure as the dependent variable to solve fluid flow in variably saturated porous media with permeability or hydraulic conductivity of the medium and any of the chosen or user-defined water retention models.

EC and RE as with other interfaces of COMSOL Multiphysics software have default and user-assigned nodes for specifying parameters and initial and boundary conditions. In this application, five EC interface nodes are used. The main node is current conservation, a default node that contains the governing equation (Equation (24)) and provides the field for specifying the value or expression of the electrical conductivity of the medium. An electric insulation node is another default for specifying the boundaries inhibiting electric current flow. The third and fourth default nodes are the initial value and ground nodes, which are used for specifying the initial value of electric potential and assigning the boundary where the electric potential is zero, respectively. To specify the initial condition of the surface of the electrode, a user-assigned node is used, which is either the floating potential node to specify electric potential or the terminal node to specify current.

The nodes of the RE interface required for the hydro-geoelectric earthing model are the No Flow node, a default node for specifying the boundary impermeable to fluid flow, Initial Value node, a default for specifying the initial value of pressure, pressure head, or hydraulic head within the porous medium, either the Pressure node, Pressure Head node, or Hydraulic Head node for specifying value or representing the initial value of pressure, pressure head, or hydraulic head at the ground surface. Richards’ Equation node is the main and default node of the RE interface which contains the governing equation (Equation (18)). This node also provides the interface for specifying the properties of the fluid and porous media such as density, viscosity, and hydraulic conductivity, specifying the storage model, saturated and residual volumetric water contents, and parameters of water retention models.

An essential aspect of the model builder is the Definition branch. This branch enables the definition of variables, parameters, and functions and assigns infinite element domain (IED) to boundaries. It also allows other tasks like geometric level selection, integration, and variable probes for domains, boundaries, or points.

A domain probe average with variable name $W_c$ is performed for the soil domain or region of influence of the earth electrode to determine the volumetric water content or effective saturation. This probe evaluates the average volumetric water content or degree of saturation of soil for different equilibrated values of soil suction. This is required to determine the bulk electrical conductivity of the soil medium.

The EC interface requires the electrical conductivity of the soil medium, as an expression or value. Equation (28) or Equation (29), which represent the electrical conductivity of the medium and provide coupling for the two physics interfaces, can be expressed as an analytic function of $W_c$ with Econ as a function name. A user-defined electrical conductivity Econ ($W_c$) may be specified in the isotropic electrical conductivity subfield of the Constitutive Jc-E relation field that is under the Current conservation node. Econ ($W_c$) estimates the electrical conductivity of the soil or medium at any degree of soil saturation.

2.1.3. Initial and Boundary Conditions

The initial and boundary conditions used in these nodes are described as follows:

• Initial condition of soil domain
  • The soil is initially assumed to be saturated, i.e., ${\mathsf{\theta }=\mathsf{\theta }}_s$. This implies that the initial value of pressure, the hydraulic head, or the pressure head in the soil domain is zero. This condition is specified at the Initial Value node of the REI.
  • Before the dissipation of electric current, the initial value of the electric potential within the soil domain is zero, i.e., $V_p=0V$. This condition is specified at the Initial Value node of the EC interface.
• Boundary conditions
  • Boundary at infinity
    The boundary of the soil domain impermeable to water flux is defined by the governing equation $-\mathbf{n}·{\rho }_{w}u=0$; where $\mathbf{n}$ is the vector normal to the boundary. For water flow, the No Flow node of the REI is used to specify the no flow boundaries such as the right and bottom sides of Figure 3. At boundaries near infinity, the effect of electric current is substantially diminished and electric potential is assumed zero, i.e., $V_p=0V$. The Ground node of the EC interface is used to specify such boundaries.
  • Ground surface
    The top boundary of Figure 3 is the ground surface interfacing with the ground soil and atmosphere. For electric current flow, the ground surface acts as an insulator impermeable to electric current flow. This means that the current density normal to this boundary is zero, i.e., $\mathbf{n}·\mathit{J}=0$, where $\mathbf{n}$ is normal to the ground surface. This boundary is assigned in the Electric Insulation node of the EC interface.
    The pressure exerted on the ground surface by the atmosphere produces a change in the energy potential of soil water that becomes constant within the unsaturated soil volume. The energy state of soil water at this boundary is specified using any of the following nodes in the RE interface: Pressure, Hydraulic Head, or Pressure Head nodes. For time-dependent analysis, the value of pressure and hydraulic/pressure head assigned must be consistent with the corresponding initial value assigned in the Initial Value node of the RE interface.
  • Soil–electrode interface
    The soil–electrode interface defines the boundary where the total electric current density $(\mathit{J}$) normal to the surface of the buried electrode is equal to the current, $I_f$, dispersed from the electrode surface, $S$, to the soil domain. The governing equation at the soil–electrode surface is $\int (-\mathbf{n}·\mathit{J})dS=I_f$: where $\mathbf{n}$ is normal to the electrode surface. This boundary condition is specified as terminal current or current $I_f$ in the Terminal or Floating potential nodes, respectively. Since earthing resistance is independent of the current dispersed from the surface of the buried electrode, any value may be assigned to $I_f$. A value of $I_f=1A$ was used in this modelling analysis.

2.1.4. Geometric Meshing/Discretization

The numerical solution of differential equations requires the discretization or meshing of the space. This process generates small-sized triangular, quadrilateral, or tetrahedron elements of the space dimension, which enables a linearized solution of the differential equation at their vertices. For 2D space, discretization results in triangular and/or quadrilateral elements. The radius of the soil domain and rod is 25 m and 9.525 mm, respectively; therefore, the dimension of the soil is 2624 times more than the rod. Due to this significant dimensional difference, two levels of discretization are required to achieve a good numerical solution.

Free triangular meshes were formed for the rod surface and soil domain resulting in minimum and maximum mesh element (free triangles) sizes of 3.175 mm and 9.525 mm, respectively. Figure 4 represents the discretization of the 20 m $\times$ 25 m or 500 m² rectangular soil space, which produced 4715 free triangles, 585 edge elements, 34 vertex elements, and 2531 mesh vertices at an element ratio of 2.947 $\times {10}^{-5}$. The minimum and average element quality, defined in terms of skewness, maximum angle, volume versus circumradius, volume versus length, growth rate, and curved skewness, the minimum and average element quality, are (0.5278, 0.8257), (0.6745, 0.907), (0.6395, 0.9306), (0.7903, 0.9531), (0.4507, 0.8126), and (0.5278, 0.8256). The discretization is good since the values of the minimum and average element qualities are all >0.1.

Figure 4 is the meshing plot of the discretization process of the 2D axisymmetric in COMSOL Multiphysics. It shows the relative size of mesh elements for the rod–soil domain. COMSOL Multiphysics assigns colour codes and values to distinguish between good (green or 1) and bad (red or 0) mesh elements. For the mesh plot, the least mesh element has a value of 0.58.

2.1.5. Study Mode and Solvers

The numerical solution of the differential equations at steady state is achieved with the Stationary study node of COMSOL Multiphysics. A preferred solver can also be selected and set within the node. Multifrontal Massively Parallel Sparse (MUMPS) direct solver is the default solver, but the Parallel Sparse Direct Solver (PARDIS) converges more often and faster than MUMPS. Within the setting window of the stationary sub-node, parametric and/or auxiliary sweeps can be specified to determine the solution for a range of values of geometric parameters or variables of domain property without modifying the model. The auxiliary sweep was performed for a range of pressure heads spanning from saturation to an unsaturated state.

2.1.6. Post Simulation Analysis

The major post-simulation analysis involves the evaluation of the resistance of the earthing rod for the range of soil suction at equilibrium. The formula for this computation is the energetic formula proposed by [43] and expressed in Equation (33).

$$R_{ei}=\frac{{\left(V_{e,i}-V_{r,i}\right)}^2}{P_i}$$

(33)

where $R_{ei}$ is the earthing resistance of an earthing electrode in $i$ layers, $i$ represents the layer index of the soil, $V_{e,i}$ is the electric potential at the electrode surface within the $i$th layer of soil, $V_{r,i}$is the electric potential at lateral distance $r$ from the electrode surface in the $i$th soil layer, $P_i$ is the power dissipated by the electrode in the $i$th soil layer. As in Equation (34), $P_i$ is equal to the volume integral of the product of ${\sigma }_i$, the electrical conductivity of $i$th soil layer, and the square of electric field strength, $E_{ir}$, at distance lateral $r$ from the electrode in the $i$th layer.

$$P_i={\sigma }_i\int {\left|E_{ir}\right|}^{2}dV$$

(34)

Since homogeneous soils are single-layered, the layer index $i$ remains invariant. Given that the electrical conductivity of the soil domain is assumed isotropic and a function of $S_e$, the degree of saturation of the soil, then Equation (33) becomes

$$R_e={\left(V_e-V_r\right)}^2{\left[\sigma \left(S_e\right)\int {\left|E_r\right|}^{2}dV\right]}^{-1}$$

(35)

The resistance computed from the simulation model can be compared with the existing analytical resistance formula. Equation (36) is Tagg’s formula for computing the resistance of vertical earth rods in homogeneous soil of isotropic resistivity, $\rho$ [44]. The only modification to the original formulation of Tagg’s equation is that $\rho ={\sigma }^{-1}$ is replaced with $\sigma \left(S_e\right)$.

$$R=\frac{1}{2\pi L}\left(\ln \frac{4r}{L}-1\right)\frac{1}{\sigma \left(S_e\right)}$$

(36)

3. Results

3.1. Effect of Soil Hydraulic Property on Earthing Performance

Vertical earth rods in three homogeneous soil textures (clayey, silty, and sandy loam) are built to illustrate the effect of the difference in hydraulic and electrical properties of the soils on earthing resistance and performance of identical earthing systems for a range of soil suction. For the same soil suction, the resistance of the vertical rod and ground surface potential due to current dispersion from the electrode in each of the homogeneous soils is evaluated and compared for corresponding soil suction. The hydraulic and electrical parameters for the clayey, silty, and sandy soils are listed in

Table 1. Electrical and hydraulic properties of clayey, silty, and sandy loam soils.

			Soil Texture
Property	Soil Parameters	Clayey †	Silty *	Sandy Loam *
	${\mathsf{\theta }}_s$$[\mathrm{m}{-}^3·$m−3]	0.5100	0.4600	0.4100
	${\mathsf{\theta }}_r$$[\mathrm{m}{-}^3·$m−3]	0.1020	0.0340	0.0650
Hydraulic	$n$ [-]	1.0900	1.3700	1.8900
	${\alpha }_h$ [m−1]	2.1000	1.6000	7.5000
	$K_s$$[\mathrm{m}·$s−1]	$5.56\times {10}^{-7}$	$6.94\times {10}^{-7}$	$1.2278\times {10}^{-8}$
	$A_{\sigma }$ [S/m]	0.1456	0.0648	0.007386
Electrical **	$\beta$ [-]	1.7530	1.1320	2.2470

${\mathsf{\theta }}_s$ is the volumetric water content at saturation, ${\mathsf{\theta }}_r$ is residual volumetric water content, $n$ is a van Genuchten fitting parameter for the water retention curve, ${\alpha }_h$ is the inverse of pressure head at the point of air entry, $K_s$ is hydraulic conductivity at saturation, and $A_{\mathsf{\sigma }}$ and $\beta$ are the coefficient and exponent of a form of Archie’s expression for electrical conductivity. The rows marked with ** representing electrical conductivity are sourced from reference [35], while the columns marked with ^† and * are sourced from reference [38] and [39], respectively.

. As indicated, these parameters were sourced from references [35,38,39].

Figure 5a is the plot of the electric conductivity while Figure 5b is the soil water retention curve for clayey, silty, and sandy loam soils using the soil parameters listed in Table 1 and Equations (30) and (10), respectively. Figure 5b shows the SWRC for silty, clayey, and sandy loam soils for a complete range of the degree of saturation. Sandy soils are predominantly coarse textured, with large but few pores between relatively large particle sizes (0.05–2 mm) and have higher hydraulic conductivity, at least seventeen times more than clay and silt textures. These large pores make water retention in sandy soils difficult even for a small increase in soil suction. This allows sandy soils to transition quickly from the state of saturation to unsaturation over a small soil suction range. Therefore, sandy soils attain their residual states quicker than clay and silty soils. Their lack of good water retention capacity makes sandy-dominated soils difficult for earthing application. They quickly attain high resistivity resulting in excessively high earthing resistance.

Conversely, clay soils are finer with smaller soil particle sizes $(\le 2 \mathsf{\mu }\mathrm{m})$. These particles are closely packed, restricting the passage of liquid. Hence, clay soils have relatively lower hydraulic conductivity and can retain more water under high soil suction, which enables them to attain the residual state relatively slower than coarse soils. High water retention property makes the resistance of earthing systems buried in a clayey environment consistent for different seasons. This property also makes clay-rich compounds like bentonite clay suitable for enhancing earthing performance in sandy or rocky areas. Silty soils are also classified as fine-grained, with particle size bigger than clay but less than sandy soils. Their exhibited water retention characteristics are similar to clay soils.

3.1.1. Estimation of Soil Resistivity for a Range of Soil Suction

Figure 6a is the combined plot of volumetric water content and electrical resistivity of silty and clayey soils for ground surface suctions equilibrated within the soil matrix. Figure 6b is similar to Figure 6a but for sandy loam soil. These plots were created from the result of auxiliary sweeps performed for a range of soil suction (negative pressure) spanning between −0.981 and −98100 kPa. Some of the suctions within this range correspond to suction saturation, −1.0712 kPa (−0.109465 m) for silty soil and −1.758 kPa (−0.17984 m) for clayey soils, and at field capacity, −33 kPa (−3.365 m), and permanent wilting point, −1500 kPa (−152.76 m).

The graphs in Figure 6a,b are curves showing the variation in volumetric water content and soil electrical resistivity for a range of soil suction for clay and silty soils and sandy loam soil, respectively. Each plot indicates that as soil suction increases, water content decreases while resistivity increases. When any of these soils are saturated, soil resistivity is fairly constant as the curve is approximately parallel to the horizontal axis. As the soil begins to de-saturate, there is a sharp variation in water content and electrical resistivity beyond the air entry point. While this trend is evident in the capillary or transition region of each curve, the rate of change of soil resistivity for the same suction value differs for these soils. It may be inferred that like their SWRC, these curves are characteristic of soil textures.

Table 2. Electrical resistivity $\rho$(Ωm) of clayey, silty, and sandy loam soils and corresponding resistance $R\left(\mathsf{\Omega }\right)$ of earthing rod at selected soil suctions.

Soil Suction		Clay Soil		Silty Soil		Sandy Loam Soil
m	kPa	$\mathit{\rho }$ (Ωm)	$\mathit{R}\left(\mathsf{\Omega }\right)$	$\mathit{\rho }$ (Ωm)	$\mathit{R}$ (Ω)	$\mathit{\rho }$ (Ωm)	$\mathit{R}$ (Ω)
0.10	0.981	6.87	2.17	15.43	4.88	135.42	42.79
3.37	33.06	7.35	2.32	16.40	5.18	188.88	59.68
10.0	98.10	9.59	3.03	21.20	6.70	537.12	169.73
50.0	490.5	32.37	10.23	87.65	27.70	$1.2\times {10}^7$	$3.7\times {10}^6$
100.0	981.0	43.14	13.63	123.61	39.06	$6.1\times {10}^7$	$1.9\times {10}^7$
152.8	1500	50.72	16.03	150.05	47.42	$1.6\times {10}^8$	$4.9\times {10}^7$
300.0	2943	65.04	20.55	201.95	63.82	$6.4\times {10}^8$	$2.0\times {10}^8$
600.0	5886	83.42	26.36	271.93	85.93	$2.7\times {10}^9$	$8.4\times {10}^8$
1000	9810	100.02	31.61	337.76	106.73	$7.9\times {10}^9$	$2.4\times {10}^9$

shows the resistivity of the three homogeneous soils following the numerical solution of Richards’ equation and water retention models at different values of soil suction. Using Equation (35), the corresponding earthing resistance was also evaluated. The post-simulation analysis shows that as the soil desaturates, the magnitude and rate of change in resistance of the vertical rod in clay soil is relatively smaller than in silty soil and is excessively high in sandy soils. The excessive value of soil resistivity for sandy soil produced proportional values of earthing resistance that are beyond the acceptable limit of resistance for safe operation. Given that the degree of saturation at field capacity (FC) represents a partially saturated state of the soil, the corresponding resistivity of soil and resistance at FC can be used as a reference for saturated and unsaturated changes.

The results of Table 2 indicate that due to the tendency for sandy soil to drain fast under the slightest increase in soil suction, earthing resistance in sandy soil is consistently high. This is a significant problem for earthing applications involving sandy soil [30].

Table 3. Comparison of the percentage change in soil suction ($\%∆$ψ) and resistance ($\%∆R$) for vertical earth rod in homogeneous clayey and silty soils for a range of soil suctions ψ [kPa] relative to suction and resistance at field capacity (−33 kPa).

Soil Suction		Clay Soil		Silty Soil
ψ (kPa)	$\mathit{\%}∆$ψ	R	$\mathit{\%}∆\mathit{R}$	$\mathit{R}$	$\mathit{\%}∆\mathit{R}$
33.06	-	2.32	-	5.18	-
98.10	197	3.03	30.60	6.7	29.34
490.5	1384	10.23	340.95	27.7	434.75
981.0	2867	13.63	487.50	39.06	654.05
1500	4434	16.03	590.95	47.42	815.44
2943	8802	20.55	785.78	63.82	1132.05
5886	17704	26.36	1036.21	85.93	1558.88
9810	29574	31.61	1262.50	106.73	1960.42

contains the percentage resistance change relative to the resistance at field capacity for a range of soil suction: −33 kPa (−3.365 m) and −9806.81 kPa (−1000 m). By associating the change in soil suction with the change in resistance, the influence of texture on the response of soil to desaturation conditions can be estimated. For example, the 197% increment of suction form field capacity, 32.99 kPa (3.365 m) to 98.068 kPa, yields 30.6% and 29.34% increase in the resistance of the earthing rod in changing from 2.32 Ω to 3.03 Ω in clayey soils and 5.18 Ω to 6.7 Ω in silty soils, respectively.

For high orders of suction change in the homogeneous soil textures, there is a corresponding change in earthing resistance that is small in clayey soils, moderate in silty soils, and large in sandy soils.

3.1.2. Effect of Soil Condition on Ground Surface Electric Potential

The electric potential of the ground surface due to current dispersion from the earth electrode depends on the distance from the current dissipating earth electrode and soil resistivity.

Figure 2 represents the plots of electric potential at the ground surface as the electrode dissipates the current of 1A in silty soil under (a) low suction and (b) high suction. Both plots contain solid line and dashed line electric potential curves. The solid and dashed lines of Figure 7a represent the ground surface potential for suctions of −1.0712 kPa (−0.10947 m) and −33 kPa (−3.374 m), respectively. Similarly, the solid and dashed curves of Figure 7b represent the ground surface potential for high suction values of −1500 kPa (−152.76 m) and −9810 kPa (−1000 m), respectively.

The plots indicate that at any radial point away from the electrode surface, the electric potential of the ground surface is reduced by a consistent percentage of the maximum electric potential for that suction. For instance, if radial distances are measured from the surface of the earthing rod, then the electric potentials at radial points of distances 1 m, 5 m, 7 m, 10 m, 15 m, 20 m, and 25 m drops consistently by 72.45, 92.51, 95.04, 96.51, 98.41, 99.33, and 99.91% of the maximum electric potential regardless of the suction state of the soil.

During these low suction values, the soil resistivity is 15.43 Ωm and 16.40 Ωm, respectively, resulting in a maximum value of ground surface electric potentials of 4.88 V and 5.18 V, respectively. For high suction values of −1500 kPa and −9810 kPa, the resistivity of silty soil increased to values of 150.05 Ωm and 337.76 Ωm, resulting in a maximum ground surface electric potential of 47.42 V and 106.73 V, respectively. Therefore, as the soil under high suction states desaturates, the electric potential at any point on the ground surface increases proportionally to the increase in soil resistivity. This outcome is consistent with the existing analytical formula for electric potential at any point on the ground.

4. Discussion

From the result, variation in soil water content due to the relative changes in the condition of soil can be used to estimate the resistivity of soil, earthing resistance, and electric potential at any point on the ground. The volumetric water content and resistivity of the homogeneous silty soil are, respectively, 0.46 m³/m³ and 15.43 Ωm at saturation, 0.4376 m³/m³ and 16.41 Ωm at field capacity (FC), 0.0911 m³/m³ and 150.05 Ωm at permanent wilting point (PWP), and 0.06189 m³/m³ and 337.76 Ωm at suction of −9810 kPa. Through this volumetric water content, the percentage increase in soil resistivity relative to resistivity at saturation is 6.35% for FC, 872.46% for PWP, and 20,988.9% for −9810 kPa. Similarly, the electric potential of the ground surface at any point is 1.06, 9.72, and 21.88 times higher at FC, PWP, and −9810 kPa than the corresponding saturation value. The percentage increase in transitioning from saturation to any suction value produces the same percentage increment to soil resistivity and ground surface potential. Therefore, the resistance of an earthing electrode in homogeneous soil and ground surface potential can be estimated using this proportional percentage increment in suction at saturation.

A major limitation in the presentation of the proposed method is that the results are not compared or validated with actual measurements. However, the accuracy of the model of electric current conservation is determined by comparing the resistance computed using the energetic approach (Joules’ law) of the finite element method (FEM) with resistance evaluated using the existing analytical formula. In the case of the vertical rod in homogeneous isotropic soil, the percentage error between the energetic and theoretical resistance evaluations is approximately 3.06%. Based on similar analysis and comparison, and an error of less than 10%, the model describing electric current conservation in the soil is satisfactory [40,43,45].

The uniqueness of the post-processing analysis result for the three homogeneous soils, subjected to the same soil condition, demonstrates that the texture of the soil an earthing system is buried in influences the resistance and performance of earthing systems. Secondly, it indicates that modelling the soil with appropriate electrical and hydraulic properties (water retention characteristics and hydraulic conductivity) enables the modelling, estimation, and prediction of the variability in earthing resistance for different soil or ground conditions. Thirdly, if the condition of the soil is regulated by external environmental and atmospheric influences, then by modelling the impact of these influences as soil suction, the seasonal values and fluctuations in earthing resistance and performance can be estimated and predicted from the tripartite relationship of soil water content, the electrical resistivity of soil, and soil suction. Fourthly, a utility company with tens of thousands of earthing sites will save enormous resources and time by using this method to estimate the seasonal variation in soil resistivity and/or earthing resistance and to excessively high values rather than performing actual measurements on the earthing sites.

Consequently, the method is an effective means of evaluating the impact of soil texture, season, and climate on the resistance and performance of earthing systems and any deviation from actual measurement can be improved with model refinement and the use of accurate parameters of the soil properties.

Although lacking the validation of measurements from actual installations, the modelling efficacy of the proposed method was compared to similar FEM models that have been reported in the literature. Using the same parameters of soil (100 Ωm) and cylindrical earth rod (radius 0.05 m, length 5 m) as in [40], the resistance computed using the electrical part of the proposed model is 14.98 Ω, which represents an error of 5.66% and 1.47% for the resistance evaluated theoretically and reported in the reference literature, respectively. This indicates that the proposed model has appreciable accuracy that can still be improved upon.

5. Conclusions

By incorporating hydraulic properties of the soil and atmospheric data, this paper introduced a new hydro-geophysical method for modelling and evaluating the resistance and performance of earthing systems in different soils and soil conditions.

Under high soil suction, clay and silty soils retain more water, while sandy soil drains very quickly, retaining little or no water within its pores. For large changes in soil suction, clay and silty soils experience relatively little volumetric water content change, unlike sandy soils that have significant volumetric water content change. Therefore, fine-grain soils (clay and silt) have high water retention and low hydraulic conductivity, while coarse-grain soils (sandy soils) have high hydraulic conductivity and relatively low water retention.

Unlike the methods of ANN and correlation/correction factors, this method has the following advantages. First, it provides an effective means of estimating the performance of earthing systems in different climatic conditions. By incorporating the properties of the subsoil and features of the earthing site, this method is an alternative method for evaluating and predicting seasonal variation in earthing resistance and performance. Consequently, the method can provide a definite means of comparing the performance and resistance reduction capacity of earthing enhancement materials in different soils and soil conditions.