Investigation on soil resistivity of two-layer soil structures using ﬁnite element analysis method

Effective grounding is the most important aspect in every electrical installation because it can ensure efﬁcient equipment operation in the event of electrical faults. Measurement data of soil resistivity is used to assess the suitability of the soil for grounding purpose and also to evaluate the structures of the soil. However, the limitation of soil resistivity measurement is various investigations on two-layer soil structures cannot be performed due to the variation of natural soil structures is limited. Therefore, in this work, modelling of soil resistivity of two-layer soil structures was proposed using ﬁnite element analysis method to estimate the depth of upper layer of non-homogeneous soil based on the measurement results from Wenner and Schlumberger methods. Also, using the proposed method, the effect of soil resistivity of each layer on the overall two-layer soil resistivity was investigated. To validate the results obtained from the proposed method, comparison between the FEA results and measured ﬁeld data was made. It was found that the average error between both results is < 3%. From the model that has been developed, the soil properties can be varied, which enhances the information that can be obtained from two-layer soil structures.


INTRODUCTION
Good grounding system is very important in ensuring a proper equipment operation in the case of faults caused by thunderstorm, accidental short circuits or insulation failure within the electrical system [1][2][3]. Grounding system depends on the soil resistivity, which is affected by various factors such as the soil type, temperature, moisture content, mineral content and compactness. For an effective grounding system, it should be designed such that it can withstand the worst possible conditions. By measuring the soil resistivity at a particular area and considering the effect of various factors, a valuable insight on how the desired resistance value can be achieved and maintained over the life of the installation with minimum cost and effort, can be obtained. Therefore, it is important to perform soil resistivity measurement under various conditions. Moisture content is one of the major factors that greatly influences the soil resistivity. Soil with higher moisture content will This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. have a lower soil resistivity. However, the soil resistivity is less affected once the moisture content exceeds approximately 20% [4][5][6][7][8][9]. In area with high soil resistivity, the application of salt in a water solution can reduce the soil resistance. However, this will enhance the metal corrosion rate, which can damage the grounding rod [5,6,10]. For temperature above freezing point, its effect on soil resistivity is small but for temperature below the freezing point, the water in the soil starts to freeze and soil resistivity increases greatly [4,5,[11][12][13].
Several works have been performed on soil resistivity measurement in different number of soil layers [14][15][16][17][18]. The effect of buried metallic structure and nearby pipelines on the soil resistivity measurement depends on the measurement method, location of the measurement and the soil structure [19,20]. The effect of different test frequency on soil resistivity measurement is not significant especially on soil with high moisture content as most of the gap holes in the soil are filled by water. However, for soil with low moisture content, frequency variation of the current source will affect the soil resistivity measurement [4,[21][22][23]. Different types of soil have different resistivity due to their compositions [24,25].
Factors such as salt and water content in the soil play an important role on the soil resistivity, which affects the grounding system design of a facility [17,23,[26][27][28]. If the soil structure is not suitable for grounding of a facility, the grounding point has to be located far away from the facility, increasing the overall installation cost. Thus, treatment of soil such as modifying the salt and water content is performed to decrease the soil resistivity value for better grounding condition.
Several mathematical models related to soil resistivity have also been proposed since the past [29][30][31][32][33]. These include modelling of tree-layer stratified soil and the electrostatic equipotential surfaces to analyse the influence of concrete buried structures in grounding systems [33]. A simplified expression for uniform soil resistivity calculation that can represent a twolayer stratified was also proposed [30]. Equations for estimating the apparent resistivity with different depth and layer thickness in a multi-layer earth structure were proposed for geotechnical investigations [31]. However, estimating the depth of upper layer of non-homogeneous soil and investigation on the effect of soil layer on soil resistivity using finite element analysis method are less likely to be found in literature.
Comparison of earth resistance measurement based on the IEEE standard-81 between a non-operational utility substation and imitative earthing systems built to emulate real earthing system was investigated under lightning impulse in [34]. The investigation found that the post-ionisation earth resistance of the imitative earthing system was found to be similar under lightning response and low magnitude currents.
Numerous soil-related works using finite element methods have been performed in the past. These include analysis of behaviour of soils under cyclic loading [35], a finite element model for the coupled flow of heat and moisture in soils under atmospheric forcing [36], analysis of integral pavement/soilwall structures in soft soil [37] and equivalent dynamic infinite element for soil-structure interaction [38]. However, these works do not focus on two-layer soil structures.
Earthing system was also investigated using finite element method (FEM) in the past. The measurement results of earthing systems at steady-state condition are similar to that of simulation from FEM [39]. However, the results are not applicable for soil resistivity measurement. The earthing systems under transients were found to be non-linear. Another work reported that the calculated earth resistance of a gas insulated sub-station of a 275/132 kV system using mathematical formulas was compared with the field measurement results [40]. The results were within reasonable agreement with each other.
Soil resistivity measurement can be used to evaluate the structures of the soil, such as the upper and lower layer of soil resistivity and the thickness of the upper layer. However, the limitation is various investigations on two-layer soil structures cannot be performed via soil resistivity measurement due to the natural soil structures are limited. Therefore, in this work, modelling of soil resistivity of two-layer soil structures was proposed using finite element analysis (FEA) method to estimate the depth of upper layer of non-homogeneous soil based on the measurement results from Wenner and Schlumberger methods. Also, using the proposed method, the effect of soil resistivity of each layer on the overall two-layer soil resistivity was investigated. From the model that has been developed, the soil properties can be varied, which enhances the information that can be obtained from two-layer soil structures.

MEASUREMENT TECHNIQUES
In selecting the measuring technique for soil resistivity, factors such as the maximum electrode depth, length of cables required, ease of the measuring technique, cost and ease of interpretation of the data need to be considered [41]. There are many techniques available for soil resistivity measurement. However, in this work, Wenner and Schlumberger methods have been chosen because they are the most commonly used methods. Measurements of soil resistivity were performed at three different locations, sandy loamy soil (Location 1), muddy soil (Location 2) and loamy soil (Location 3).

Wenner method
In Wenner method (Figure 1), four electrodes of equal size are driven into the soil surface in a straight line at equal spacing a [15]. A current I is injected between the outer pair of the electrodes and the potential difference V is measured between the two inner pair of the electrodes [41]. Using Ohm's Law, the resistance between the two potential electrodes can be determined. This value is used to calculate the resistivity of the soil. If the spacing a between the four electrodes is considerably small, the measured soil resistivity is the result of the resistivity of the upper layer of the soil. To accurately obtain the resistivity of this layer, more measurements with small increment of a are required. If the spacing a is increased, the current will penetrate deeper into the soil. The measured soil resistivity will be an indication of the average soil resistivity over a larger area [41]. If the penetration depth b is very small compared to a and R is the resistance, the resistivity ρ can be calculated using [16,25]

Schlumberger method
In Schlumberger method (Figure 2), the spacing between four electrodes is not equal. There are two cases for the electrodes spacing; fixing the spacing between the potential electrodes d or fixing the spacing between current and potential electrodes c. The former case requires only two current electrodes to be moved. The latter case requires all electrodes to be moved. The soil resistivity for Schlumberger method is calculated using [25] where ρ s is the soil resistivity, c is the spacing between the current and potential electrodes, d is the spacing between two potential electrodes and R is the measured resistance. Figure 3 shows the measurement results of soil resistivity using Wenner and Schlumberger methods for different distances between the current electrode spacing on sandy loamy soil. From both methods, it is noticed that the soil resistivity increases initially up to a certain level and then decreases when the spacing between the electrodes is increased. The soil in the area where measurement was conducted is likely to be nonhomogeneous because the soil resistivity measured varies widely with the position of the electrodes [42]. The results using Wenner method can be explained by treating the soil as a three-layer model, with a higher resistivity layer between two lower resistivity layers. For a small electrode spacing a, almost the entire applied current will flow within the upper layer. Thus, the measured soil resistivity will be closer to the resistivity of the upper layer. When a is increased, the current will penetrate deeper into the soil, which is the second layer of the soil. As the current-flow lines reach the interface between the two layers, they are distorted by the presence of the subsurface boundary and flow towards the least resistive path, in this case towards the first layer. Since current is flowing in the second layer, which has a higher resistivity, it steadily increases the measured soil resistivity as the average resistivity of the two layers. The resistivity keeps increasing with larger a, until the maximum value, where most of the current flows in the second layer. As a is increased further, the current will penetrate deeper into the third layer, causing it to be refracted towards the lower resistivity layer, which is the third layer. Thus, the measured soil resistivity starts to decrease steadily.

Location 1 (sandy loamy soil)
From the measurement results of soil resistivity using Schlumberger method by fixing the potential electrode spacing d with increasing the electrode distance c, the measured soil resistivity initially increases as c is larger. Then, it reaches the maximum value and decreases as the value of c is increased further. For small spacing of c, the current flows mainly in the lower resistivity layer (first layer). For medium spacing of c, the current flows mainly in the higher resistivity layer (second layer) and for larger spacing of c, the current flows mainly in the lower resistivity layer (third layer).
From Figure 3, both methods yield the same response for the same distance of the current electrodes. However, it is noticed that the soil resistivity measured by Schlumberger method is slightly lower than the Wenner method for current electrodes spacing greater than 15 m. This is due to when the current electrode spacing exceeds 15 m, for example 30 m, the current penetration for both methods is the same. However, in Schlumberger method, the two potential electrodes are closer to each other than in Wenner method. This results in a lower current density around two potential electrodes because most of the current flows into deeper layer of the soil. Hence, the measured soil resistivity using Schlumberger method could give a better indication of the soil resistivity of the lower layer. However, if the spacing is too large, the potential difference between the potential electrodes will be very small. Hence, a more sensitive measuring meter is required.

Location 2 (muddy soil)
The measurement results of soil resistivity using Wenner and Schlumberger methods for different electrode spacing on muddy soil are shown in Figure 4. From the Wenner method results, the resistivity of the soil increases as the distance a between the electrodes is increased. The result can be explained by treating the soil as a two-layer soil model, with a lower resistivity soil layer overlaying a higher resistivity soil layer [42]. For small distance of a, in this case a < 1 m, almost the entire applied current will flow within the upper layer. As a result, the measured soil resistivity will be closer to the resistivity of that layer.
As a is increased, the current will penetrate deeper into the soil, which is the second layer. Since the second layer has a higher soil resistivity than the first layer, the current-flow lines will be distorted at the sub-surface boundary and flow towards the least resistive path, which is the first layer. However, there is a current flow in the second layer, which has a higher soil resistivity. Hence, the measured resistivity of the soil will start to increase from the resistivity of the first layer. Since current tends to flow in the first layer, it is expected that the measured resistivity of the soil approaches the resistivity of the second layer more slowly. For Schlumberger method by fixing the potential electrode spacing d but changing the electrode distance c, the measured soil resistivity steadily increases as the electrode spacing c is increased, similar to the results obtained using Wenner method. However, for current electrode spacing greater than 15 m, the soil resistivity measured by Schlumberger method is slightly higher than the Wenner method. This could indicate that the resistivity of the deeper layer is measured by Schlumberger method with larger spacing of the current electrodes. Figure 5 shows the measurement results of soil resistivity using Wenner and Schlumberger methods at loamy soil. For Wenner method, the measured resistivity of the soil decreases as the distance a between the electrodes is increased. The result can be explained by treating the soil as a two-layer soil model with a higher resistivity layer overlaying a lower resistivity layer [42]. For small a, most of the current will flow within the upper layer. This results in the measured soil resistivity to be closer to the resistivity of the upper layer. As a is increased, the current will penetrate deeper into the soil. At the boundary between the two layers, the current will be distorted and flow towards the least resistive path, which is the bottom layer of the soil. Hence, the measured soil resistivity will start to decrease and approach the resistivity of the bottom layer as the value of a is increased. Since current tends to flow in the bottom layer, it is expected that the measured resistivity of the soil will approach the resistivity of the bottom layer faster.

Location 3 (loamy soil)
For measurement results of soil resistivity using Schlumberger method by fixing the potential electrode spacing d but varying the electrode spacing c, the measured soil resistivity decreases as c is increased. The results obtained using Schlumberger method is within good agreement with the results obtained using Wenner method. It is also noticed that for larger current electrode spacing, Schlumberger method gives a lower soil resistivity than Wenner method.

FEA MODELLING OF SOIL STRUCTURES
Homogeneous and two-layer soil models were developed to reproduce the soil resistivity values that were obtained from the measurement. The proposed model geometry of the soil resistivity measurement was built using finite element analysis software (FEA) software, COMSOL Multiphysics in a 3D space. The reason of using FEA in this work is for extension into determination of the effect of soil treatment, buried metallic structures, impurities and voids in soil under investigation and for multiple soil layers on soil resistivity value for future work. These investigations could not be performed using analytic method.
FEA is the most widely used for solving engineering and mathematical model problems. FEA method is used to solve partial differential equations in two or three space variables. The FEA subdivides a large system into smaller elements called finite elements to solve a problem by discretising a particular space into the space dimensions. A mesh of the object is constructed, which is the numerical domain for the solution, with a finite number of points. Then, the FEA method formulates the boundary value problem to yield a system of algebraic equations by approximating the unknown function over the domain. The equations that model the finite elements are assembled to form a larger equation system, which models the whole problem.
In order to solve the problem in this work by FEA, the governing partial differential equation is derived from the where E is the electric field, ρ is the resistivity, σ is the conductivity and ε is the permittivity. The relation between the electric potential V and E is defined by Inserting Equation (2) into Equation (1), Rewriting Equation (3) in differential equation, the equation used to solve the electric potential in the model is [43,44] where V is the electric potential, ε is the permittivity and σ is the conductivity. Referring to Equation (4), the equation is solved within each of the mesh elements to determine the electric potential V according the value of σ and ε that are assigned in the FEA model. Figure 6 shows the flowchart of the FEA method and how simulation was done using FEA method in this work. The simulation was done by first, the model geometry was constructed in FEA software. Then, the material properties and boundary conditions were assigned. After that, the model was meshed. Finally, the model was solved by Equation (4) to obtain the electric potential and field distribution in the geometry. A model geometry of homogenous soil, which consists of four points to represent four electrodes and a rectangular soil For the homogenous soil model, the soil resistivity was set as 100 Ω m to test the validity of the proposed calculation method of soil resistivity. If the voltage determined from this model is divided by the injected current and is converted to resistivity using Equation (1) or (2) equals to 100 Ω m, the proposed model is considered acceptable. Once the model is found to be acceptable, the calculation is extended into two-layer soil structure model.

Two-layer soil model in FEA
A homogeneous soil model can only be used when the variation in the measured soil resistivity is small. However, this case is less likely to be found in real situation as soil is rarely homogeneous. If there is a large variation in the measured soil resistivity, the use of homogeneous soil model is likely to yield inaccurate results. Hence, a two-layer soil model can be used. Figure 8 shows a complete model geometry of the proposed two-layer soil structure in FEA, which consists of an upper layer with a finite depth h, a lower layer with infinite depth and four points representing four electrodes. The boundary conditions assigned are similar to that of homogeneous soil model. The parameters assigned for every material in two-layer soil model are shown in Table 1. The calculation of the soil resistivity is similar to that of homogenous soil FEA model.   Figure 9 shows the electric potential distribution between the two point sources from the simulation. The red colour point is where the current is injected into the soil while the blue point is where the current sinks. The potential electrodes are somewhere in between the two current electrodes. From this distribution, the electric field in the soil further from the current electrodes is uniform. The simulation results from this model were validated using the soil resistivity equation of Wenner and Schlumberger methods. Table 2 shows the results of soil resistivity in 100 Ω m homogeneous soil FEA model using Wenner and Schlumberger methods. The soil resistance R was obtained from the FEA model and this value was used in Equations (1) and (2). From this table, the calculated soil resistivity values are found to be almost similar to the assigned resistivity value (100 Ω m) of the soil in the FEA model. The average error between the soil resistivity using Equations (1) and (2) and the assigned resistivity x-axis (m) y-axis (m)

FIGURE 9
Electric potential distribution for homogeneous soil model value (100 Ω m) is 0.1% and 0.04% for Wenner and Schlumberger methods respectively. Hence, the homogeneous soil FEA model that has been developed can be considered acceptable.

Two-layer soil FEA model
The concept of soil resistivity calculation for homogenous soil FEA model was adopted for two-layer soil FEA model. In this study, four cases were simulated; (1) with the lower layer fixed at ρ 2 = 300 Ωm and upper layer with varying resistivity ρ 1 ; (2) with the upper layer fixed at ρ 1 = 100 Ω m and lower layer with varying resistivity ρ 2 , (3) the thickness of the upper layer is varied, and; (4) the current electrode distance is varied. Case (1): Figure 10a shows the simulated overall soil resistivity vs soil resistivity of the upper layer ρ 1 for two-layer soil FEA model, where the soil resistivity of the lower layer ρ 2 is fixed at 300 Ω m, the depth of the upper layer h is 1 m and the distance between the current electrodes is 18 m. In this case, when ρ 1 is increased, the amount of current that flows into the lower layer of the soil increases because higher ρ 1 reduces the current flow through it. Hence, the soil resistivity value detected at the measuring point becomes nearer to the value of the soil resistivity at lower layer ρ 2 , which is 300 Ω m.
Case (2): Figure 10b shows the simulation results of soil resistivity vs ρ 2 for two-layer soil FEA model, where ρ 1 is fixed at 100 Ωm, h at 1 m and the distance between the current electrodes at 18 m. When ρ 2 increases, the overall soil resistivity becomes higher. Although higher soil resistivity at certain layer will limit the current flow through it, which in turn reduces the overall soil resistivity, this behaviour is not shown in Figure 9b. Hence, it can be deduced that when h = 1 m and the distance between the current electrodes is 18 m, the overall soil resistivity is for the lower layer of the soil. Hence, when ρ 2 is increased, the detected soil resistivity also increases.
Case (3): The overall soil resistivity vs depth of the upper layer h from the FEA model is shown in Figure 11 for two cases: (a) ρ 1 (300 Ω m) > ρ 2 (100 Ω m) and (b) ρ 2 (300 Ω m) > ρ 1 (100 Ω m). The current electrode distance was fixed at 18 m. For both cases, the overall soil resistivity seems to be equal to the soil resistivity of the lower layer. This can be seen from ρ 1 (300 Ω m) > ρ 2 (100 Ω m), where the overall soil resistivity equals to ≈100 Ω m, which is the soil resistivity of the lower layer. For ρ 2 (300 Ω m) > ρ 1 (100 Ω m), the overall soil resistivity equals to ≈300 Ω m, which is the soil resistivity of the lower layer. For ρ 1 (300 Ω m) > ρ 2 (100 Ω m), when h is increased, the overall soil resistivity increases. As the depth of the upper layer of the soil increases, more current is occupying the space in the upper layer. Since the soil resistivity of the upper layer of the soil is 300 Ω m, the overall soil resistivity approaches 300 Ω m when the depth of the upper layer is increased further. However, for ρ 2 (300 Ω m) > ρ 1 (100 Ω m), when h is increased, more current is flowing in the upper layer of the soil because the resistivity is lower. Hence, the overall soil resistivity decreases.
Case (4): Figure 12 shows the overall soil resistivity vs current electrode distance for ρ 1 > ρ 2 and ρ 2 > ρ 1 from the FEA model, where the depth of the upper layer of the soil h was kept  Simulated soil resistivity vs upper layer soil depth h for ρ 1 > ρ 2 and ρ 2 > ρ 1 at 1 m. For ρ 1 (300 Ω m) > ρ 2 (100 Ω m), when the current electrode distance is increased, the overall soil resistivity was found to be lower. This is due to the current penetrates deeper into the soil layer. Hence, the overall soil resistivity is approaching the soil resistivity of the lower layer ρ 2 (100 Ω m) when the current electrode distance is increased. However, for ρ 2 (300 Ω m) > ρ 1 (100 Ω m), the overall soil resistivity increases with larger current electrode distance, approaching the soil resistivity of the lower layer ρ 2 (300 Ω m). Again, this is due to the current penetrates deeper into the soil layer.

Comparison between measurement data and results from FEA model
From the measurement results that have been obtained in this work using Wenner method, from Figure 4 (where the soil  resistivity measurement was performed on muddy soil), the soil structure can be treated as a two-layer soil, which a lower resistivity layer is overlaying a higher resistivity layer, ρ 2 > ρ 1 [42,45]. The measured soil resistivity is approaching 195 Ω m when the distance between the current electrodes is further. Hence, the soil resistivity of the lower layer ρ 2 is set as 195 Ω m in the twolayer soil FEA model. The measured soil resistivity is around 105 Ω m when the distance between the current electrodes are very near to each other. Hence, the soil resistivity of the upper layer ρ 1 is set as 105 Ω m in the two-layer soil FEA model.
Once the values of ρ 1 = 105 Ω m and ρ 2 = 195 Ω m were set in the FEA model, the thickness of the upper layer h was obtained by tuning its value until the mean square error (MSE) between the measurement and simulation results is the lowest. Table 3 shows the MSE for different h from the simulation results for this case. The value of h was increased from 1 m with a step of 0.2 m until a local minima was found. The tuning of h was refined between 1.2 and 1.6 m because the lowest MSE was found to be within this range. It was found that the value of h that yields the lowest MSE is 1.5 m. Figure 13 shows the measurement and simulation results of soil resistivity vs distance between the current electrodes for the two-layer soil model for ρ 2 > ρ 1 or muddy soil. Referring to this figure, both measurement and simulation results are within good agreement with each other. The errors between the measurement data and results from the FEA model are 24.60 (MSE) and 1.35% (average).  The same method used for Figure 4 was used for the measurement results shown in Figure 5. From Figure 5, where the soil resistivity measurement was performed on loamy soil, the soil structure can also be treated as two-layer soil, where a higher resistivity layer is overlaying a lower resistivity layer, ρ 1 > ρ 2 [42]. From the measurement results using Wenner method, the measured soil resistivity is approaching 108 Ω m when the distance between the current electrodes is further. Hence, the soil resistivity of the lower layer ρ 2 is set as 108 Ω m in the two-layer soil FEA model. The measured soil resistivity is around 295 Ω m when the distance between the current electrodes is very near to each other. Hence, the soil resistivity of the upper layer ρ 1 is set as 295 Ω m in the two-layer soil FEA model.
Once the values of ρ 1 = 295 Ω m and ρ 2 = 108 Ω m were set for the FEA model, the thickness of the upper layer h was obtained by tuning its value until the mean square error (MSE) between the measurement and simulation results is the lowest, similar to the steps taken for Table 3. Table 4 shows the MSE for different h from the simulation results for this case. The value of h was increased from 1 m with a step of 0.2 m until a local minima was found. The tuning of h was refined between 1.2 and 1.6 m because the lowest MSE was found to be within this range. It was found that the value of h that yields the lowest MSE is 1.4 m. Figure 14 shows the measurement and simulation results of soil resistivity vs distance between the current electrodes for the two-layer soil structure for ρ 1 > ρ 2 or loamy soil. Again, both measurement and simulation results are within good agreement with each other. The errors between the measurement data and results from the FEA model are 75.17 (MSE) and 2.56% (average).

CONCLUSIONS
Homogenous and two-layer soil structure models have been successfully developed using FEA method. The models were able to estimate the depth of upper layer of non-homogeneous soil based on the measurement results from Wenner and Schlumberger methods. Also, using the proposed method, the effect of soil resistivity of each layer on the overall two-layer soil resistivity was successfully investigated. It was found that the overall soil resistivity increases linearly with the resistivity of the lower layer of the soil. However, the overall resistivity increases logarithmically with the resistivity of the upper layer of the soil. From the comparison between the measurement field data of soil resistivity and the results obtained from the FEA model, the average error between both results was found to be <3%. Therefore, the proposed FEA model developed in this work can be considered acceptable for investigation on soil resistivity of two-layer soil structures. From these models, soil properties can be varied, enhancing the information that can be obtained from two-layer soil structures. The proposed models can also be extended into determination of the effect of soil treatment, buried metallic structures, impurities and voids in soil under investigation and for multiple soil layers on soil resistivity value.