## 1 Introduction  
The safety of electric devices connected to a network depends on the appropriate dimensioning of the grounding systems. In the case of a technical failure or lightning strike, it is the grounding system that ensures the safety of people and animals. This means that the step voltage and touch voltage should not increase beyond the limit permitted by the standard. Larger grounding systems are often designed using the finite element method (FEM) [1; 2; 3; 4; 5; 6; 7; 8; 9; 10; 11; 12]. With the FEM model, the potential distribution in and on the soil's surface can be determined, and the step voltage and touch voltage are determined based on the calculated potentials. An FEM model can only be made if the data of the soil structure are known in advance. This information is usually obtained through measurements. Various measurement methods are used, the most common of which is the Wenner method [13; 14; 15; 16; 17]. A generalized method can also be found in the literature [14], and the Schlumberger method is also mentioned in [17]. This article focuses on the Wenner method. Determining the parameters of the soil is a difficult task because it requires an inverse approach. Analytical horizontal multi-layer soil models [15; 16; 17; 18; 19] are of great help in this task. The soil data are obtained using the Wenner method, in which the apparent resistivity is determined for different electrode spacings. It is necessary to determine such parameters of analytical models of soil with two or more layers [15; 16; 17; 18; 19] so that the calculated resistances coincide with the measured ones. Due to the complexity of the problem, various approaches are used, such as the steepest descent method, the Levenberg-Marquardt method, the Newton method, the generalized inverse method, and the quasi-Newton method, as shown in [20]. The use of the second-order gradient technique, the use of electrostatic images, and the integrated-circuit-based fast relaxed vector fitting approach are shown in [21; 22; 23]. We decided to solve this problem with modern metaheuristics [18; 24; 25; 26; 27; 28], because they are known for their usefulness in engineering problems. A problem arises when the soil is not homogeneous and the horizontally layered soil model poorly represents the actual conditions in the soil. This article shows an approach that partially avoids this weakness. As testing the presented approach requires detailed soil structure information, a basic FEM model of inhomogeneous soil was created, on which [PERSON]'s measurements were simulated. The final FEM model was also made, which was obtained on the basis of the presented procedure. [PERSON]'s measurements were also simulated on this and compared with those from the basic model. In order for the method to be successful, reliable determination of the model parameters is necessary. In order to find the most suitable one, we used different metaheuristics and compared their performance with each other. The methods used were Differential Evolution (DE) with three different strategies of DE/rand/1/exp, DE/rand/2/exp, and DE/best/1/bin; Teaching-Learning-Based Optimization (TLBO); and Artificial Bee Colony (ABC). Our contributions in this work are as follows:  
* This article presents an approach that enables the use of an analytical horizontally layered soil model for inhomogeneous soil. In the approach, the observation area is not easily divided into smaller parts, because the deeper layers of the soil would not be considered according to the smaller distances between the electrodes when measuring with the Wenner method. The presented approach was not found in the literature. * Three evolutionary methods were compared, DE (using three strategies), TLBO, and ABC, with the aim of obtaining the most suitable one. The quality of the evolutionary methods was also analyzed in terms of convergence and calculation time. * Memory usage was tested in order to speed up the calculations.