## 2 Modeling Grounding Electrodes with Internal Resistance: Backgrounds  
The calculation of the grounding resistance of an electrode needs the model for the electrical conductivity of the soil to be known. From a general point of view, it would be necessary to solve the equation \(\overset{\rightarrow}{\
abla}\cdot(\sigma(\overset{\rightarrow}{r})\overset{ \rightarrow}{\
abla}\phi(\overset{\rightarrow}{r}))=0\) together with some boundary conditions, where \(\phi(\overset{\rightarrow}{r})\) is the absolute potential at the point \(\overset{\rightarrow}{r}\) and \(\sigma(\overset{\rightarrow}{r})\) represents the conductivity of the ground. For a multilayered model it is possible to find acceptable semi-analytical solutions for the potential. This model assumes that conductivity is just a non-continuous piecewise constant function. The soil is composed of horizontal layers of infinite extension and specific thickness in which the electrical conductivity takes a constant and different value in each layer. Thus, the actual conductivity is replaced by a step function although in practice, due to the difficulty of the calculations, only a few layers are considered. The soil conductivity profiles of interest for this work comprise layers of almost any conductivity with thicknesses ranging from a few centimeters to tens of meters. Although for applications in grounding design it is common to deal with models in which the top layer is a few meters thick in which to place the grounding electrode. Choosing a constant conductivity in each layer allows calculating the potential from the equation for a generic conductivity. The connection between the solutions corresponding to each layer is made by imposing the continuity of the potential at the different interfaces as well as the continuity of the normal component of the current density through interfaces. Due to the supposed cylindrical symmetry of the potential, which is associated with conductivity dependent only on the depth as the z coordinate, the procedure to solve the equation for the potential relies on using the method of separation of variables in a cylindrical coordinate system. Thus, the potential created in layer \(j\) by an electric point current source of strength \(I\) in layer \(i\) could be written as [12]  
\[\phi_{ij}=\frac{\rho_{i}I}{4\pi\left|\overrightarrow{r}_{i}-\overrightarrow{r} _{j}\right|}\delta_{ij}+\int\limits_{0}^{\infty}{(f_{ij}(\lambda)e^{-\lambda( z_{i}-z_{j})}+g_{ij}(\lambda)e^{\lambda(z_{i}-z_{j})})J_{0}(\lambda r)d\lambda} \tag{1}\]  
where \(\delta_{ij}\) is the Kronecker delta and the unknown functions \(f_{ij}\) and \(g_{ij}\) need to be calculated by imposing the boundary conditions at each interface, namely \(\overrightarrow{\
abla}\phi_{i1}\cdot\overrightarrow{n}_{z}|_{z=0}=0\) which stands for the null current flux through the soil surface and  
\[ \begin{split}\phi_{ij}(\overrightarrow{r}_{j,j+1})=\phi_{ij+1}( \overrightarrow{r}_{j,j+1})\\ \sigma_{j}\cdot\overrightarrow{\
abla}\phi_{ij}\cdot\overrightarrow {n}\bigg{|}_{\overrightarrow{r}_{j,j+1}}=\sigma_{j+1}\cdot\overrightarrow {\
abla}\phi_{ij+1}\cdot\overrightarrow{n}\bigg{|}_{\overrightarrow{r}_{j,j +1}}\end{split}  \tag{2}\]  
what guarantees continuity of potential and normal current flux conservation through the interface separating the layers \(j\) and \(j+1\). In (1) and (2), \(\rho_{i}\), is the constant resistivity of layer \(i\), \(r=\sqrt{\left(x_{i}-x_{j}\right)^{2}+\left(y_{i}-y_{j}\right)^{2}}\) and \(J_{0}\) is the zero-order Bessel function. In (1) the functions\(f_{ij}(\lambda)\) and \(g_{ij}(\lambda)\) must be calculated by (2) at the above-mentioned interfaces. The calculation of the functions \(f\) and \(g\) can be difficult and the evaluation of the integral in (1) is only approximate due to the oscillatory character of the \(J_{0}\) Bessel function, in addition to the own singularities of the \(f\) and \(g\) functions [13].