## 2 Materials and Methods for Coupled Thermal-Electrical Modeling of the Grounding Electrode
### The Finite Differences Method for Heat Equation (4)  
To numerically solve Equation (4), this can be achieved using a finite difference scheme in cylindrical coordinates, with the radial variable as the only operative spatial coordinate. Figure 2: Reduction of the initial system to a simpler geometry. In this way, a cross-section of the pipe is considered, where \(R_{1}<r<R_{2}\), dividing this domain into N \(-1\) segments of size \(\Delta r=\frac{R_{2}-R_{1}}{N-1}\), with a total of N nodes \(r_{i}\), \(i=1\ldots\) N, \(r_{1}=R_{1}\), and \(r_{N}=R_{2}\). Similarly, the chosen time interval [\(t_{i}\),\(t_{j}\)] is divided into M intervals of size \(\Delta t=\frac{t_{f}-t_{i}}{M-1}\), with \(t_{1}=t_{i}\) and \(t_{M}=t_{f}\). Equation (4) is discretized as follows:  
\[ \begin{array}{l}T_{i}^{n+1}-T_{i}^{n}=\frac{\alpha_{3}\Delta t}{\Delta r^{2}} (T_{i-1}^{n}-2T_{i}^{n}+T_{i+1}^{n})+\frac{\alpha_{2}\Delta t}{\Delta r^{2} \Delta r}(T_{i+1}^{n}-T_{i-1}^{n})+\frac{\alpha_{3}\Delta t}{\Delta r}\dot{ \varepsilon}_{2}(r_{i})\quad\text{for $R_{1}<r<R_{2}$}\\ -\frac{k_{2}}{\Delta r}(T_{2}^{n}-T_{1}^{n})_{r=R_{1}}=-\frac{k_{1}R_{1}}{2R_ {1}}\frac{(T_{1}^{n+1}-T_{1}^{n})_{r=R_{1}}}{\Delta t}+\frac{R_{1}}{2}\dot{ \varepsilon}_{1}\quad\text{for $r=R_{1}$}\\ T_{N}^{n+1}=T_{s}\quad\text{for $r=R_{2}$}\end{array}  \tag{6}\]  
The convergence is determined by the condition \(\frac{\alpha_{3}\Delta t}{\Delta r^{2}}\leq 0.5\), which limits the value of the time step \(\Delta t\) taken for the evolution of the system. From Equation (6), the first one applies to the internal nodes. For the nodes adjacent to the boundary, when the temperatures at the boundary appear in the calculation, the second equation representing the Neumann boundary condition must be applied. The third equation, Equation (6), represents a Dirichlet boundary condition indicating that the temperature of the boundary node is always Ts and, thus, \(T_{N}^{n}=T_{s}\). This boundary condition can be considered valid as long as the evolution of the system is studied for times that are sufficiently short so that the boundary \(r=R_{2}\) does not show significant changes from its initial state. Alternatively, for the node adjacent to the boundary \(R_{1}\), that is, node 2, node 1 is located at the left boundary and will satisfy the second Equation (6), so  
\[T_{1}^{n+1}=\frac{2\alpha_{1}k_{2}\Delta t}{k_{1}R_{1}\Delta r}(T_{2}^{n}+( \frac{k_{1}R_{1}\Delta r}{2\alpha_{1}k_{2}\Delta t}-1)T_{1}^{n}+\frac{R_{1} \Delta r}{2k_{2}}\dot{\varepsilon}_{1}) \tag{7}\]  
With these considerations, the equation for the internal nodes becomes  
\[T_{i}^{n+1}=T_{i}^{n}(1-\frac{2\alpha_{2}\Delta t}{\Delta r^{2}})+T_{i-1}^{n} (\frac{\alpha_{2}\Delta t}{\Delta r^{2}}-\frac{\alpha_{2}\Delta t}{2r_{i} \Delta r})+T_{i+1}^{n}(\frac{\alpha_{2}\Delta t}{\Delta r^{2}}+\frac{\alpha_{ 2}\Delta t}{2r_{i}\Delta r})+\frac{\alpha_{2}\Delta t}{k_{2}}\dot{\varepsilon }_{2}(r_{i}) \tag{8}\]  
which can be simplified somewhat if \(C=\frac{\alpha_{2}\Delta t}{\Delta r^{2}}\) is introduced, so that (8) becomes  
\[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{ \varepsilon}_{2}(r_{i}) \tag{9}\]  
An important comment must be made at this point regarding (6). The Neuman boundary condition at \(r=R_{1}\) has been discretized at a low order (\(o(\Delta r)\)). We can achieve greater accuracy (\(o(\Delta r^{2})\)) by using the so-called ghost nodes method outside the domain.