New Detailed Modeling of GICs in the Spanish Power Transmission Grid

The threat of Geomagnetically Induced Currents (GICs) driven by severe Space Weather looms over technological systems such as power grids. Assessing their vulnerability is thus vital to avoid damages or even disruption of the electrical power supply. This endeavor, however, entails an interdisciplinary approach, ranging from the characterization of the geoelectrical structure of the Earth beneath and around the area of interest, or the modeling of the power network from its parameters and topology, and including the validation of the modeling process by means of (direct or indirect) GIC flow measurements. In this paper, we summarize our current achievements focused on mainland Spain, concentrating on the improvements reached after going from a homogeneous Earth's resistivity to an alternative 3D electrical resistivity distribution approach to geoelectric field computation, which is still in progress because new empirical impedance tensors are needed, mainly at sites in the west of the Iberian Peninsula. The second major achievement has come from the addition of the 220 kV level to the network model. The overall improvement has been validated against real GIC data in one area of the country. The new vulnerability maps show that in some nodes the predicted GIC has been substantially reduced by the sum of both effects. The assessment has been carried out down to the level of the individual windings of each transformer, and examples of the estimated GIC flow are given for substations with numerous power transmission lines converging to them at diverse orientations.

• Improvements are achieved by using a 3D inversion of magnetotelluric (MT) data to derive the electric field and adding the 220 kV level to the network model • More accurate geomagnetically induced current (GIC) estimates for all power transmission lines, nodes and transformers are provided • The system experiences a reduction in GICs flowing to ground, but at the transformer winding level they exhibit, in general, an augmentation and/or increased corrosion of pipeline steel, which in turn may have an impact on the critical infrastructures that provide services to the population such as electric power, gas or oil delivery and train systems.
The induced geoelectric fields have periods typically ranging from tens of seconds to several hours (or frequencies in the range of 10 −1 -10 −4 Hz). These electric fields are conveyed to conductors running on the surface such as power lines, whose transformers are generally grounded working in an alternating current (AC) regime (50)(51)(52)(53)(54)(55)(56)(57)(58)(59)(60). Because the GICs are quasi-direct currents (DC) flowing in response to the geoelectric field, they may disrupt or permanently damage these transformers.
These currents are responsible for both dramatic system-level failures, such as the famous Hydro-Québec 1989 outage that affected millions of people, and less severe but far more frequent effects such as pipeline corrosion (Gummow & Eng, 2002;Pulkkinen et al., 2001), railway support system failures (Boteler, 2021) and power transformer wear. Although traditionally viewed as a high-latitude phenomenon, during the Halloween storms of 2003 reports indicated important GIC impacts in South Africa (Gaunt & Coetzee, 2007), where geomagnetic latitudes range between 25°S and 35°S. After those reports and our first analysis in north-eastern Spain (Torta et al., 2012), vulnerability assessments of the GIC hazards in mid-to-low latitude countries have proliferated all over the world. Kelbert (2020) provides a complete review of the state of the GIC studies around the world.
GIC hazard assessment in a system does not require precise knowledge of GIC values, but rough estimates of its magnitude are sufficient (Pirjola, 2008). However, this assertion must not be misunderstood and convey to lax approaches when constructing approximate models. Only using the highest voltage level tends to overestimate the predicted currents, and this could favor falling into certain catastrophism to be avoided. In the same way, ocean-land interfaces or the juxtaposition of provinces with notably different conductive structure could produce geographically distributed differences up to a factor of 1,000 in the derived geoelectric fields (e.g., Bedrosian & Love, 2015), that can be masked when using solely homogeneous or 1D resistivity Earth structures.
With the above considerations in mind, after a decade of effort and continuous improvements, in this paper we present our latest results on the modeling of GICs in Spain. One of our main achievements has been the replacement of the homogeneous Earth's resistivity by a 3D lithospheric-scale geoelectric model in the Iberian Peninsula. Next, our efforts have been directed at determining how the derived geoelectric fields interact with the power transmission grid infrastructure by including the 220 kV voltage level. Both achievements have been validated by comparing model results with actual GIC measurements under particular power transmission lines.

New 3D Lithospheric Resistivity Model
Aimed at improving our previous analyses on the threat posed by GICs to the Spanish grid, we have undertaken a critical review of the existing magnetotelluric (MT) data and models and created a new 3D Electrical Resistivity Model of the Iberian Lithosphere (ERMIL). Several local and regional 2D and 3D crustal and lithospheric electrical resistivity models had been obtained in the last 25 years (Alves Ribeiro et al., 2017;Campanyà et al., 2012Campanyà et al., , 2018Carbonell et al., 1998;Garcia et al., 2015;Ledo et al., 2000;Marti et al., 2009;Monteiro Santos et al., 2002;Pous et al., 1995Pous et al., , 19971999, 2004Rosell et al., 2011). However, a lithospheric scale 3D model of the whole peninsula was still lacking. The available magnetotelluric sites ( Figure 1) have an irregular distribution, being the NE and SE part of the peninsula the areas with a denser MT sites distribution. The western part of the Iberian Peninsula does not contain enough MT sites to obtain a reliable electrical resistivity model at lithospheric scale. New long period data are being collected, and existing long period data already acquired by other authors are being retrieved to complement the model presented in this work in the near future.
The initial model was set for four different units, namely sediments, upper crust, lower crust, lithospheric mantle and asthenosphere, as well as the seawater around the peninsula. The limits between different units were obtained from the EUcrust07 model (Tesauro et al., 2008) and Fullea et al. (2021) for the lithosphere-asthenosphere boundary. It includes the bathymetry and the topography for the area. The resistivity values used for each of those units were determined from previous 2D and local 3D MT models (see references above). Figure 2 shows cross sections of the initial 3D model. For the sediments, a value of 50 ohm·m was chosen, 5,000 ohm·m for the crust, 1,000 ohm·m for the lithospheric mantle, 30 ohm·m for the asthenosphere (not shown in the figure) and 0.3 ohm·m for the seawater.
The inversion of the 56 sites of Figure 1 was done using the ModEM 3D inversion code (Egbert & Kelbert, 2012;Kelbert et al., 2014) inverting the four components of the impedance tensor at 19 logarithmic scale equispaced periods from 1 to 30,000 s. Due to noise issues some of these values were rejected, giving rise to a total number of 2,603 values used. An error floor of 10% and 15% of · xy yx Z Z was chosen for the anti-diagonal and diagonal components of the impedance tensor, respectively. A metric to characterize the global fit of the model responses to the measured data is the normalized root mean square error, defined as: where d obs,i and d calc,i refer to the observed and calculated model responses, respectively, of the different components of the impedance tensor, and e i is the aforementioned error floor of the components and N is the total number of data used. Mesh dimensions of the model were 103 × 122 × 119, the horizontal size of the mesh of the whole peninsula was 11 × 11 km and the vertical elements above and below the bathymetry/topography changes were 100 m. The initial normalized RMS error was 40 and the final one after 150 TORTA ET AL.

10.1029/2021SW002805
3 of 17 iterations was 2.1. Comparisons between observed data and model responses can be found in Figure S1. Figure 3 shows cross sections of the final ERMIL.
The final model keeps the sharp limit between the sediments and the upper crust although the resistivity value of the crust and mantle structures shows some differences. The final model eliminates the electrical difference between the upper and lower crust, resulting in an overall moderately resistive crust. Low resistivity structures appear at the Pyrenean (C1 in Figure 3) and the Betics orogens (C2 in Figure 3) as well as below the Iberian Chain (C3 in Figure 3). Those low resistivity structures have been reported previously by local studies (Campanyà et al., 2018;Rosell et al., 2011).
Another difference with the initial model is the contact between the electrical lithosphere-asthenosphere boundary (e-LAB) that is a bit shallower than the initial one based on other geophysical properties (Fullea et al., 2021) although the lack of data precludes a sufficiently accurate characterization.
Magnetotellurics is most sensitive to vertically integrated conductivity or conductance. For a uniform layer, the conductance is the product of conductivity and thickness. Figure 4 shows the comparison between the conductance down to 120 km depth for the initial and final 3D resistivity models (Figures 4a and 4b, respectively). We have chosen this depth, which roughly coincides with the whole lithosphere, because this is the resolution limit of the available data. The conductance of the western part of the final model is like the original one, and below the Pyrenean and Betics orogens enhances because the electrical resistivity of the lower crust below those orogens was overestimated in the initial model.

Constructing a Model of the Spanish Power Network
The Spanish power network used in Torta et al. (2012Torta et al. ( , 2014    lengths, we assume that they would not add much difference to the modeled GICs. The effect of shield wires has been neglected. The schematic power network can be seen in Figure 5. From Figure 5, although it seems that the improvement does not represent an enormous increase of substations and transmission lines, in practice we went from a total of 173 substations, 375 transformers and 300 lines to 629, 1,401, and 912, respectively (these include a series of newly created 400 kV substations and lines). The counterintuitive appearance comes from the huge concentration of substations in densely populated areas, especially around the metropolitan regions of Madrid and Barcelona, which appear superposed in the figure. As in the previous studies, we have only included the border foreign stations of France, Portugal and Morocco.

New GIC Model Results
Once the information on the extended power network was gathered and properly arranged, the new DC-equivalent network model was created, from which the GIC values at the different parts of the system were calculated by subjecting the model network to unit geoelectric fields in the North and East directions. The procedure for solving the GICs was carried out using a standalone code written in Matlab for the generation of a DC-equivalent network and GIC calculation which works for power grids with multiple voltage levels. The equivalent network thus constructed consists of as many nodes per substation as the maximum number of voltages encountered in one single substation plus one. The former nodes play the role of buses of the equivalent substations and the latter plays the role of the neutral point. Our scheme, whose detailed description will be presented in a subsequent paper, is an alternative scheme to the mathematically equivalent (Boteler and Pirjola, 2014) methods of Lehtinen and Pirjola (1985) (LP), which solves for the earthing currents at each node, and to the Nodal Admittance Matrix (NAM) method, which calculates first the nodal voltages, and the GICs in a subsequent step. Alternatively, our method solves the circuit laws for the current flowing between the bus nodes and the neutral point which, unlike LP, requires no (infinite-resistance) earth connection for the buses. Note that all of these methods require subsequent steps to calculate the TORTA ET AL.
10.1029/2021SW002805 6 of 17 GICs in the different parts of interest in the system (windings, transmission lines, and grounding current). However, the number of equations/unknowns is reduced in our method compared to the traditional ones, thus optimizing the computational cost (an advantage that becomes important for large networks) while keeping a good condition of the matrix to be inverted. In any case, we have checked that all three methods provide exactly the same GIC results.
This test allowed comparing the new results with those obtained by Torta et al. (2014) when applying uniform electric fields of 1 V/km oriented North and East and looking at the resulting GICs calculated at each of the nodes and transformer windings. Figure 6 shows the absolute values in red for the 400 kV voltage level and in blue when both the 400 and 220 kV levels are considered.
Comparing both results, we can see a general decrease in the GIC amplitude which, as expected, is mainly due to the fact that currents not only flow through the neutral points of the 400 kV circuit but are also distributed over the 220 kV lines. A very clear example can be seen in Figure 6 (bottom) in Galicia (Northwest of Spain).
To assess the maximum expected GIC in each transformer because of extreme geomagnetic storms, postevent analysis of data from the Spanish geomagnetic observatories during the October 29-31, 2003 ("Halloween storm") were performed, although other episodes coincident with very abrupt storm onsets, which have proven to be more hazardous at these mid-latitudes, were analyzed as well. The resistivity model described in section 2 allowed obtaining the impedance tensor at each point of a dense, regularly distributed grid, which was then convolved with the geomagnetic field interpolated at the same points of that grid from the field measured at the Spanish and some surrounding observatories using the spherical elementary current systems (SECS) technique as shown in Torta et al. (2017), to provide the geoelectric field. Figure S2 shows the locations of all these observatories. The interpolated geoelectric fields from that grid were integrated along the power lines to obtain the necessary line voltages, from which the GICs were ultimately calculated in the earthing points, in the power lines and in each separate winding of the substation transformers.
Movie S1 (top) shows the modeled GICs flowing to ground at each substation during the Halloween storm, with circles whose diameters are proportional to the flowing current. It is important to emphasize that we assumed the network configuration and elements currently in operation, which is possibly different from the one at the time of this event. Red/blue circles denote GICs flowing to/from the earth, respectively, while red and blue lines correspond to the 400 and 220 kV voltage levels, respectively. At the bottom, the time evolution is shown for the Ascó substation, where the model predicts the highest GIC throughout the whole time series.
Since it is the flow of the GICs through the winding of transformers that produces DC flow distorting the alternate AC flow, causing the saturation of the transformer core during the one-half cycle, it is more significant to seek the maximum expected GIC flow throughout the windings of each transformer in order to look for its dangerousness. In connection with this, it is worth mentioning that different types of transformers are not equally prone to GIC consequences, but they depend on their design (e.g., Rajput et al., 2020). Figure 7 (top) shows the peak absolute winding current at each transformer of the network during the Halloween storm, that is, the winding with the highest GIC is displayed for each transformer (either the primary or the secondary for full-winding transformers, or the series or common for autotransformers). These results contrast with those obtained with just the 400 kV grid (Figure 7, bottom). In both figures, null values correspond either to junctions (i.e., substations without transformers) or transformers at the first substation of neighboring (i.e., foreign) systems, where results are inaccurate due to the omission of the rest of the systems. Figures S3 and S4 correspond to the same Figure 7, but for the cases of using horizontal geoelectric fields of 1 V/km in the north and east directions, respectively. When adding the 220 kV level, even though the system experiences an overall reduction in GIC flowing from the network to ground or vice versa, as reflected in Figure 6, the latter results show that at the transformer winding level they exhibit, in general, an augmentation. This is a consequence of the increase in the number of transmission lines that converge at substations, sometimes with sharp changes in line orientation, causing them to exhibit a large GIC. For example, when a 220 kV line continues in the same direction as a 400 kV line that formerly ended in an edge substation when only the higher voltage was considered, it TORTA ET AL. tends to increase the GICs through its windings. This is somehow illustrated in Zheng et al. (2013) and it is the case of Muruarte station (see Figure 8), the location of which is indicated by a red arrow in Figure 6, bottom. Another example of GIC increase is seen at Mudarra substation (Figure 9), the location of which is indicated by a red arrow in Figure 6, top. Although by shifting from a system model with only 400 kV voltage level to a 400 and 220 kV model the substation earthing current diminishes, the maximum winding GIC increases. The corresponding GIC values to/from ground can be scrutinized by zooming in on those locations on the map in Figure 6, and are clearly detailed in Figures 8 and 9, along with the GICs in the series windings of the corresponding transformers, and in each of the power lines that converge on those substations. Again, our analysis here focusses on the effects of horizontal geoelectric fields of 1 V/km in the north and east directions, but, as indicated by Bernabeu (2013), the extent of an actual impact depends on the magnitude and orientation of the geoelectric field, both of which vary stochastically over the course of an actual storm.
TORTA ET AL.
10.1029/2021SW002805 9 of 17   located far enough from the line to exclude the signals derived from the GICs (as in, e.g., Hübert et al., 2020). The East-to-West line indicated by a red arrow at the northwest of EBR Observatory, for example, is a good candidate to perform such new measurements.
Validation is the latest fundamental step of the modeling process, as it allows assessing the validity of the assumptions made by comparing its output with real GIC observations. Using our own methodology for achieving GIC data in the way outlined above, we have started measuring the GICs flowing in a few power lines of the grid in regions where the ERMIL resistivity model was constructed with a dense distribution of MT sites. According to our experience, only induced currents above about 1 A give magnetic signatures that exceed the noise threshold in most deployments. As we started measuring during the solar minimum and Spain is a mid-latitude country, the latter fact limited the significance of the available recorded data, but we can already report the results for a number of minor geomagnetic storms. In Figure 11 we show them for an event in August 2020 at the power line where we obtained the best fits. It corresponds to a 400 kV line located in the southeasternmost corner of the Iberian Peninsula (indicated by the green arrow in Figure 10). Figure 11 suggests that we are progressing correctly in the model improvement, as the measured GICs are significantly better predicted using the geoelectric field produced by the final ERMIL than using the initial model. When we calculate typical skill scores for that occasion, such as the correlation coefficient, r, or TORTA ET AL.

10.1029/2021SW002805
11 of 17 Figure 9. Same as Figure 8 for Mudarra substation. Note that the geomagnetically induced current (GIC) flow at the substation is spread between two autotransformers built in parallel at the substation. Note also that the 400 kV circuit line through which the greatest current flows (connecting Mudarra substation to the North-West), in fact, includes three power lines. the performance parameter, P′ (see Marsal & Torta, 2019), the former goes from 0.83 to an almost perfect correlation of 0.97; while the latter goes from 0.44 to 0.70.
The not yet perfect match is a consequence of the model's underestimation of the measured signal, mainly originating from the numerous assumptions still made concerning the resistances of the network elements.

Discussion and Conclusion
During the last decade, we have followed a rather natural and common way to progressively approach the problem of assessing the hazard from GICs in a country. In our first steps, we established a modeling procedure to derive GIC estimates in Spain that assumed a uniformly resistive Earth model with a conductivity of 0.001 S/m, and just used the elements of the 400 kV voltage level as the load of the circuit to which the geoelectrical source of the GICs was connected (Torta et al., 2014). Using the current network configuration TORTA ET AL.

10.1029/2021SW002805
13 of 17 and all elements in operation yielded maximum transformer neutral GIC amplitudes of 70 A, and more than 100 A for the major storms of "Halloween" and March 24, 1991 (the one with the most abrupt onset ever recorded), respectively. The transmission lines in Spain cover a considerably large area, with edge-toedge lengths of the order of 1,000 km for several paths. Note that GIC flows not just in individual lines, but through the whole system, and the important length is not the "individual line length" but the "system length" (Zheng et al., 2013). However, those historical maximum amplitudes of the GICs certainly seem too large compared to those recorded or estimated for other countries with comparable geomagnetic latitudes. Unfortunately, GIC measurements were not available in Spain for the times of those historical events to validate those predictions. In any case, they apparently produced no harmful effects for the period 2000-2010, as shown by a correlation analysis between the power utility's failure logs and the geomagnetic activity (Cid et al., 2013). The performance of such initial modeling could only be analyzed by assessing its agreement with actual measurements recorded in one substation (Vandellòs) during 2011 and 2012, and we, as many others, realized that the use of empirical magnetotelluric impedance tensors is essential to improve the match between model predictions and actual observations.
Geoelectric fields are rarely measured directly but can be estimated from ground observations of the magnetic field and knowledge of the surface impedance. Torta et al. (2017) acquired MT data at one site in the proximity of the Vandellòs substation, which is located on the Mediterranean coast in the NE of Spain and considered the measured impedances as the representative response of the regional geoelectrical structure for long periods. This pilot experience of using customized values of impedance instead of generic ones proved to be effective because the agreement between GIC model predictions and measurements improved significantly in that area. This single MT site has been supplemented now by other 55 stations, mainly distributed in the NE and SE areas of the Iberian Peninsula, which allowed to perform a 3D inversion that amended an initial resistivity model. The initial model was constructed by using the available information from bathymetry, topography, sediment thickness and depth to the upper and lower crust obtained from seismic and potential field methods and assuming typical resistivity values for each layer.
The ERMIL model still needs to be complemented by the inversion of more MT measurements filling the remaining gaps, especially in the western part of the Peninsula (including Portugal), but it is starting to provide realistic results, according to the validation already performed with the GIC measurements under a power line in the SE.
Regardless of the Earth's resistivity structure used, the retrospectively predicted peak amplitudes for major historical events have in general been reduced at node level by adding the 220 kV level to the network model (see Figure 6). At the transformer previously considered as the most vulnerable (in Manzanares substation, indicated by the green arrows in that figure), the GIC drop is especially significant because that single autotransformer in the substation accumulated all the current flowing through a node that was wrongly considered as an edge node, as the 220 kV lines connected to it were omitted. This fact, along with a substantial change in the geoelectric field brought about by the replacement of an unrealistic highly resistive homogenous Earth for that region with one computed using the 3D approach, yields an extraordinary reduction in the earthing current predicted in that transformer again for the Halloween storm: now, it hardly exceeds a value of 10 A compared to the previous value of 70 A. Since GIC measurements at the neutral of that transformer are available, validation of the modeling is possible there (Figure 12), which gives satisfactory results.
The decreases have not been as dramatic everywhere, so that the maximum absolute value of GICs in the whole network to ground at node level for that event is reduced from the 78.2 A given by Torta et al. (2014) to 56.5 A. The value for this maximum absolute GIC flow to ground, which can be taken as a proxy for the vulnerability of the network, amounted up to 70.7 A on the occasion of the most abrupt storm commencement ever recorded (on March 24, 1991). As it is difficult to obtain geomagnetic data at the 1-min or lower cadence for that epoch, we just used magnetic field data from EBR (which had to be digitized from analog records using the technique described in Curto et al., 1996) to derive the geolectric field in that occasion (given the mid-latitudinal location of the Iberian Peninsula, this is a reasonable simplification). This is shown in Figure S5, while for comparison with Figure 7 (top), Figure S6 shows the peak absolute winding current at each transformer of the network during that storm. An additional relevant result of this study is that, although the GIC flowing to ground is reduced in general after inclusion of the 220 kV level, this is not necessarily the case when applied to individual transformer windings, which is arguably the most critical quantity for assessing the vulnerability of the grid to Space Weather events.
Modeling GICs in complex power transmission networks involves facing important challenges regarding the true heterogeneity of the Earth's geoelectrical structure, the true values of the network resistances, and all the different voltage levels of the network. Our results confirm preconceived notions that the simplifying assumptions that often must be made to overcome these difficulties can lead to significant errors when producing vulnerability assessments of power grid assets. There exists the possibility of approaching the problem by developing transfer functions that directly map geomagnetic field measurements to GICs in such a way that the errors by having assumed geoelectrical simple structures or simplified network characteristics are absorbed by the uncertainties of the derived parameters (Heyns et al., 2020).
However, empirical approaches need GIC data, but not all transformers in a network can be instrumented. The "traditional" modeling, such as the one outlined in this work, provide GIC estimates for all power transmission lines, nodes and transformers, although they can become inaccurate if all the necessary elements in the modeling chain are not exactly known or have substantial uncertainties. Therefore, it is recommended to progressively enhance the model by introducing more MT data and network characteristics as they become available. In any case, the analytical network modeling formulation and the assumptions made must be validated through significant GIC measurements. With the current available data, we have validated the modeling of GICs by measuring the magnetic effect of the induced currents on several 400 kV lines of the Spanish power grid. Results revealed good correlations with the GICs modeled in three of them, and especially in the Tabernas-Litoral line (as shown in Figure 11). This suggests that the 3D lithospheric resistivity model works quite well in areas where it enjoyed a dense distribution of MT data to be inverted (see Figure 1). However, the mismatch in the GIC amplitude in three of the four tested power lines could be due to the existence of still substantial uncertainties in some assumed values of transformer winding resistances and/or their grounding, or due to the effect of having ignored power lines at voltages below 220 kV. The research that led to these results was in part carried out using funds from "la Caixa" Foundation. The authors would like to thank Red Eléctrica de España (REE) for supporting this study. We are deeply thankful to J. Rupérez, M. Quintana, and L. Martínez, from REE, for providing the Spanish electrical transmission network's DC characteristics to us and attending to many questions concerning them. It should be noted that such information was provided with caveats and restrictions. We are also grateful to two anonymous reviewers for their useful comments, which helped to improve the manuscript. Some of the results presented in this paper rely on data collected at geomagnetic observatories; we thank the national institutes that support them, and INTERMAGNET for promoting high standards of geomagnetic observatory practice.