## 2 Materials and Methods for Coupled Thermal-Electrical Modeling of the Grounding Electrode
### The Finite Differences Method for Heat Equation (4)
If \(-k_{2}\frac{\partial T(r,t)}{\partial r}|_{r=R_{1}}\simeq-k_{2}(\frac{T_{1}^ {n}-T_{0}^{n}}{2\Delta r}+o(\Delta r^{2}))\), where \(T_{0}^{n}\) is the temperature of a hypothetical node to the left of the boundary node whose temperature is \(T_{1}^{n}\), then we only need to eliminate the presence of that node through the equation that now satisfies the internal node, which is (9), now being  
\[T_{1}^{n+1}=T_{1}^{n}(1-2C)+T_{0}^{n}\cdot C(1-\frac{\Delta r}{2  
\[T_{1}^{n+1}=\left[\frac{k_{2}}{2\Delta rC(1-\frac{\Delta r}{2k_{1}})}+\frac{k_{1}R_ {1}}{2\alpha_{1}\Delta t}\right]^{-1}\cdot\left\{ \begin{array}{l}T_{1}^{n} \left[\frac{k_{1}R_{1}}{2\alpha_{1}\Delta t}+\frac{k_{2}(1-2C)}{2\Delta rC(1- \frac{\Delta r}{2k_{1}})}\right]+\\ T_{2}^{n}\left[\frac{k_{2}}{2\Delta r}+\frac{k_{2}(1+\frac{\Delta r}{2k_{1}})}{2 \Delta r(1-\frac{\Delta r}{2k_{1}})}\right]+\dot{\varepsilon}_{2}(R_{1}) \left[\frac{\Delta r}{2(1-\frac{\Delta r}{2k_{1}})}\right]+\frac{R_{1}}{2} \dot{\varepsilon}_{1}\end{array} \right\} \tag{12}\]  
This procedure for refining the calculation at the boundary is equivalent to the energy balance method in elementary control volumes around the nodes, as found in [19]. It is important to note that, since we are in a transient regime with sources, the balance is a complex operation, requiring the equation to be formulated for the control volume:  
\[\int\limits_{\Delta V}\rho C_{p}\frac{\partial T}{\partial t}dV=-\oint \limits_{\Delta S}\overset{\rightarrow\rightarrow}{q}\overset{\rightarrow} {d}\overset{\rightarrow}{S}+\int\limits_{\Delta V}\dot{\varepsilon}dV \tag{13}\]  
The surface integral must be performed with the outward normal vector so that the net flux has the correct sign. At this point, it is necessary to complete the scheme by adding some other effects that may arise in practical problems. On the one hand, it is clear that the variation in the resistivities of the conductor and the soil with temperature must be taken into account, as they are found in \(\dot{\varepsilon}_{1}\) and \(\dot{\varepsilon}_{2}\). For many types of materials, including metals and small temperature ranges, this variation follows the well-known linear law \(\rho(T)=\rho(T_{0})(1+\alpha(T-T_{0}))\), where \(\alpha\) is the thermal variability coefficient of the material resistivity and \(T_{0}\) is a reference temperature, typically set at 20 \({}^{\circ}\)C or 25 \({}^{\circ}\)C. For larger temperature ranges, always below the material melting point, the above linear law is often replaced by an exponential law \(\rho(T)=\rho(T_{0})\exp(\alpha(T-T_{0}))\). Furthermore, for conductive materials such as soils, rocks, and similar materials, it is known that electrical resistivity decreases with temperature. In [20], the following law can be found:  
\[\rho(T)=\frac{\rho(T_{0})}{1+\alpha(T-T_{0})}\quad T>T_{0} \tag{14}\]  
with typical values of \(\alpha=0.019\)\({}^{\circ}\)C\({}^{-1}\) and \(T_{0}=25\)\({}^{\circ}\)C. In the same previous reference, the expression can also be found:  
\[\rho(T)=\rho(T_{0})(\alpha+\beta\exp(-\frac{T}{\gamma}))\quad T>T_{0} \tag{15}\]  
where \(\alpha=0.447\), \(\beta=1.4034\), \(\gamma=26.815\)\({}^{\circ}\)C\({}^{-1}\), and \(T_{0}=25\)\({}^{\circ}\)C. Both expressions provide acceptable temperature-dependent decreasing values, although (15) exhibits a more moderate decrease than (14). To complete the variability of the system parameters with temperature, the possible variation of the thermal conductivity of the materials with temperature could be included as a second effect. Here, the model proposed in [21] is cited:  
\[ \begin{array}{l}k(T)=\frac{k_{0}}{0.99+T(a-\frac{b}{b_{0}})}\\ k_{0}=0.53k(25)+0.5\sqrt{1.13k(25)^{2}-0.42k(25)}\end{array}  \tag{16}\]  
In (16), \(T\) is the absolute temperature, and the parameters have the values \(a=0.003\) K\({}^{-1}\) and \(b=0.0042\) W/mK\({}^{2}\).