## 2 Wenner's Method and Description of Horizontally Layered Soil Model
### Description of Horizontally Layered Soil Model  
The results obtained with the Wenner method are only specific resistivities as a function of distance and do not provide information about the soil structure independently. They need to be linked to horizontally multi-layered analytical soil models developed in the past [15; 16; 17; 18; 19]. It is necessary to determine the parameters of analytical soil models so that the calculated resistivities fit the measured data best. An inverse problem is described, which can be solved successfully using optimization methods. \(N\) is the number of soil model layers, which can be two or more. In the rest of this article, a three-layered model is used, so Figure 2 shows a general \(N\)-layered and a three-layered soil model that was used. Each soil layer has a specific soil resistance, marked from \(\rho_{1}\) to \(\rho_{N}\), and a layer thickness marked from \(h_{1}\) to \(h_{N-1}\). The deepest, the \(N\)th layer, has an infinite depth. The soil is homogeneous within a single layer. The analytic expression of apparent resistance depends on distance \(d\) (shown in Figure 1), and it is written in (3) [29]. \[\rho=\rho_{1}\bigg{(}1+2d\int_{0}^{\infty}f(\lambda)[J_{0}(\lambda d)-J_{0}(2 \lambda d)]\,d\lambda\bigg{)} \tag{3}\]  
Figure 2: (a) \(N\)- and (b) three-layered soil model. \(\rho_{1}\) is the specific resistivity of the first layer, \(J_{0}\) is Bessel's function of the first kind, and \(d\) is the distance between the electrodes. Equation (3) contains the function \(f(\lambda)\), which is calculated using (4). \[f(\lambda)=\alpha_{1}(\lambda)-1 \tag{4}\]  
\(\alpha_{1}\) for the \(N\)-layered model is calculated with the sequence of equations written in (5). \[ \begin{array}{c}K_{1}(\lambda)=\frac{\rho_{2}\alpha_{2}(\lambda)-\rho_{1}}{ \rho_{2}\alpha_{2}(\lambda)+\rho_{1}};\qquad\alpha_{1}(\lambda)=1+\frac{2K_{1 }e^{-2\lambda b_{1}}}{1-K_{1}e^{-2\lambda b_{1}}}\\ \qquad\qquad\cdots\\ K_{N-2}(\lambda)=\frac{\rho_{N-1}\alpha_{N-1}(\lambda)-\rho_{N-2}}{\rho_{N-1} \alpha_{N-1}(\lambda)+\rho_{N-2}};\ \alpha_{N-2}(\lambda)=1+\frac{2K_{N-2}e^{-2\lambda b_{N-2}}}{1-K_{N-2}e^{-2 \lambda b_{N-2}}}\\ K_{N-1}(\lambda)=\frac{\rho_{N}-\rho_{N-1}}{\rho_{N}+\rho_{N-1}};\qquad\alpha_ {N-1}(\lambda)=1+\frac{2K_{N-1}e^{-2\lambda b_{N-1}}}{1-K_{N-1}e^{-2\lambda b_ {N-1}}}\end{array}  \tag{5}\]  
In the continuation of the paper, a three-layered model is used and, based on (3) for the three-layered model, \(\alpha_{1}\) is calculated with the equations presented in (6). \[\begin{array}{c}K_{1}(\lambda)=\frac{\rho_{2}\alpha_{2}(\lambda)-\rho_{1}}{ \rho_{2}\alpha_{2}(\lambda)+\rho_{1}};\qquad\alpha_{1}(\lambda)=1+\frac{2K_{1 }e^{-2\lambda b_{1}}}{1  
The approach of using two or more directions of measurements is shown in Figure 3a. The average values of measurements are calculated for all directions and then are used in the whole area of interest. The disadvantage of this approach is that the soil resistivity is too low in some areas and, in some, too high. With that, the grounding system in the areas with excess resistivity is over-dimensioned, and in the areas with insufficient resistivity, it is dimensioned inadequately. The approach of dividing the area into smaller parts is used in Figure 3b. In this case, the average values of measurements L3 and L4 are calculated and used in the area of crossing lines L3 and L4. The same approach is also used for the areas with lines L5 and L6, L7 and L8, and L9 and L10. The strong disadvantage of this approach is that the current does not flow deep into the soil, and, with that, the resistivity of deeper soil is not considered. We could increase the area of interest presented in Figure 3b and also extend lines L3-L10 beyond the field of interest; however, in many cases, in a real environment, it is not possible to make measurements outside the presented area of interest, due to factors such as buildings.