## 2 Materials and Methods for Coupled Thermal-Electrical Modeling of the Grounding Electrode
### The Finite Differences Method for Heat Equation (4)
It can be easily verified that, for the soils considered in this work, made up of materials equivalent to clay, the variation of thermal conductivity with temperature is insignificant and, therefore, it will not be considered in the temperature profile calculation. For any choice of the model, the variability in electrical resistivities and thermal conductivities can affect the values of the longitudinal current \(l_{c}\) and the transversal current \(I_{s}\). Therefore, these two currents must be treated as new variables of the problem, with their own equations associated with the electrical circuit to which the conductor belongs, which will give the values of \(I_{c}(t)\) and \(I_{s}(t)\). Therefore, in strict terms, in (7) and (9), the power densities of the sources \(\dot{\epsilon}_{1}\) and \(\dot{\epsilon}_{2}(r_{i})\), should be replaced by \(\dot{\epsilon}_{1}^{n}(T_{i}^{n})\) and \(\dot{\epsilon}_{2}^{n}(T_{i}^{n},r_{i})\), respectively, where it is made explicit that there is a time dependence for the currents \(I_{c}\) and \(I_{s}\). In practice, the time loop for the implementation of (9) must also include the algorithm by which these currents are obtained at each time step. The numerical scheme would therefore proceed as follows:  
\[ \begin{array}{l}T_{i}^{n+1}=T_{i}^{n}(1-2C)+T_{i-1}^{n}\cdot C(1-\frac{\Delta r }{2r_{i}})+T_{i+1}^{n}\cdot C(1+\frac{\Delta r}{2r_{i}})+\frac{C\Delta r^{2}}{ k_{2}}\dot{\epsilon}_{2}(T_{i}^{n},r_{i})\quad(1<i<N)\\ T_{1}^{n+1}=\left[\frac{k_{2}}{2\Delta r\zeta(1-\frac{\Delta}{2R_{1}})}+\frac {k_{1}R_{1}}{2\Delta t\Delta t}\right]^{-1}\cdot\left\{\begin{array}{l}T_{1 }^{n}\bigg{[}\frac{k_{1}R_{1}}{2\Delta t\Delta t}+\frac{k_{2}(1-2C)}{2\Delta r \zeta(1-\frac{\Delta}{2R_{1}})}\bigg{]}+\\ T_{2}^{n}\bigg{[}\frac{k_{2}}{2\Delta r\zeta}+\frac{k_{2}(1+\frac{\Delta}{2R_{1 }})}{2\Delta(1-\frac{\Delta}{2R_{1}})}\bigg{]}+\\ \dot{\epsilon}_{2}(T_{1}^{n},R_{1})\bigg{[}\frac{\Delta r}{2(1-\frac{\Delta r }{2R_{1}})}\bigg{]}+\frac{R_{1}}{2}\dot{\epsilon}_{1}\big{(}T_{1}^{n}\big{)} \end{array}\right\}\\ T_{N}^{n+1}=T_{s}\end{array}  \tag{17}\]  
with the heat power densities defined by  
\[\dot{\epsilon}_{1}(T_{i}^{n})=\frac{\rho_{0}\exp(\alpha_{1T}(T_{i}^{n}-T_{0}) )\cdot I_{c}^{2}}{\left(\pi R_{1}\right)^{2}} \tag{18}\]  
for the sources in the conductor and  
\[ \begin{array}{l}\dot{\epsilon}_{2}(T_{i}^{n},r_{i})=\frac{\rho_{0}\exp(- \alpha_{2T}(T_{i}^{n}-T_{0}))\cdot I_{c}^{2}}{\left(2\pi L\right)^{2}}\frac{ 1}{r_{i}^{2}}\\ \dot{\epsilon}_{2}(T_{i}^{n},r_{i})=\frac{\rho_{0}\left(\alpha_{3T}+\beta_{3} \exp(-\frac{T_{3T}^{n}}{2\Delta t\Delta t})\cdot I_{c}^{2}\right)}{\left(2 \pi L\right)^{2}}\frac{1}{r_{i}^{2}}\\ \dot{\epsilon}_{2}(T_{i}^{n},r_{i})=\frac{\frac{\rho_{0}}{1+\alpha_{4T}(T_{i} ^{n}-T_{0})}\cdot I_{c}^{2}}{\left(2\pi L\right)^{2}}\frac{1}{r_{i}^{2}} \end{array}  \tag{19}\]  
for the heat sources in the conductive soil. In (18) and (19), \(\rho_{0}\) is the resistivity at the reference temperature \(T_{0}\), and the coefficients \(\alpha_{1T}\), \(\alpha_{2T}\), \(\alpha_{3T}\), \(\alpha_{4T}\), \(\beta_{3T}\), and \(\gamma_{3T}\) are associated with the variation of resistivity with temperature for the materials. In general, for the common metals used in conductors, the temperature coefficient is positive (for copper, \(\alpha_{1T}=\alpha_{Cu}=0.0039\ ^{\circ}\mathrm{C}^{-1}\)), meaning that the resistance increases with temperature. For a conductive soil, in the absence of specific additional information, the coefficients will take the typical values mentioned earlier: \(\alpha_{2T}=0.001\ ^{\circ}\mathrm{C}^{-1}\), \(\alpha_{3T}=0.447\), \(\beta_{3T}=1.4034\), \(\gamma_{3T}=26.815\ ^{\circ}\mathrm{C}\), and \(\alpha_{4T}=0.019\ ^{\circ}\mathrm{C}^{-1}\) with \(T_{0}=25\ ^{\circ}\mathrm{C}\) for models in (19). Given that determining the currents \(I_{c}\) and \(I_{s}\) from the material resistivities in the system involves a slow calculation, whereas solving (9) requires a time loop with several thousand steps, an approximation is proposed.