Geomagnetically Induced Current Model in New Zealand Across Multiple Disturbances: Validation and Extension to Non‐Monitored Transformers

Geomagnetically induced currents (GICs) produced during geomagnetic disturbances pose a risk to the safe operation of electrical power networks. One route to determine the hazard of large and extreme geomagnetic disturbances to national electrical networks is a validated model to predict GIC across the entire network. In this study, we improve upon an earlier model for New Zealand, expanding the approach to cover transformers nationwide by making use of multiple storms to develop national scaling factors. We exploit GIC observations which have been made and archived by Transpower New Zealand Ltd, the national grid operator. For some transformers the GIC observations span nearly 2 decades, while for others only a few years. GICs can vary wildly between transformers, particularly due to differences in the electrical network characteristics, transformer properties, and ground conductivity. Modeling these individual transformers is required if an accurate representation of the GIC distribution throughout the network is to be produced. Here we model the GIC during 25 disturbed periods, ranging from large geomagnetic storms to weakly active periods. We adopt the approach of scaling model output using observed GIC power spectra, finding that it improves the correlations between the maximum model and observed GIC by between 10% and 40% depending on the transformer. The modeled GIC at the 73 transformers which have measured GIC are analyzed to create local and national scaling curves. These are used to allow modeling for transformers without in‐situ GIC. We present approaches to utilize this technique for future storms, including non‐monitored transformers.

• Modeling geomagnetically induced currents (GIC) in New Zealand for multiple past geomagnetic disturbances shows good agreement with measurements at South Island transformers • Scaling the modeled GIC in the frequency domain using past GIC measurements produces more accurate modeling results • Averaged frequency scaling can be applied to transformers without measurements enhancing the models value for future extreme storm studies

Supporting Information:
Supporting Information may be found in the online version of this article.
Modeling of GIC in New Zealand have primarily used the thin-sheet model (Divett et al., , 2018. These studies have focused on validating the model GIC in New Zealand's South Island against a large spatially varying set of GIC observations. Total Harmonic Distortion (THD), a strong indicator of GIC, show significant agreement with the model GIC output in the South Island and also highlight potential areas of concern in the North Island for locations where GIC observations are not currently taken . Ingham et al. (2017) used transfer functions created from MT measurements made in New Zealand to predict GIC at a number of transformers in the South Island of New Zealand during various geomagnetic storms. This included using a simulated storm of a similar magnitude to the 1859 Carrington event, showing that a transformer at the South Dunedin substation could experience upwards of 800 A. Rodger et al. (2017) estimated peak currents during extreme geomagnetic storms by extrapolating measurements from past events using the corresponding horizontal magnetic field rate of change. They found the currents during extreme storms were very large, at orders of magnitude in agreement with those presented by Ingham et al. (2017). Other studies have looked at calculating GIC in the North Island using MT measurements (Mukhtar et al., 2020), confirming high currents in transformers where THD increases were reported by Rodger et al. (2020).
In this study, we describe a technique to improve GIC modeling capability where local GIC observations do not exist, providing realistic values at the transformer level. Section 2 describes the steps involved in modeling the transformer-level GIC induced by a varying external electric field. We use a modified version of the thin sheet model described in Divett et al. (2020). Modeling the electric field requires a ground conductance model (Section 2.1) and geomagnetic field input (Section 2.2). The geoelectric field model is described in Section 2.3. Applying the electric field to a representation of the transmission network in New Zealand produces model GIC at the transformer-level (Section 2.4). In Section 3, we describe the New Zealand network of GIC observations and the selection of 25 geomagnetic disturbance events that are identifiable with the datasets available. These are used in Section 4 to relate the model GIC output to the extensive set of GIC observations available for the events. We show how scaling the modeled GIC in the frequency domain can allow for a more accurate comparison with the observed GIC. For each transformer where GIC measurements are available we determine a local frequency domain scaling factor from the 25 events. These are used to correct the GIC model output to more realistic values for each transformer. Finally we combine all local scaling factors into a national scaling factor which can be used to correct the modeled GIC output for transformers which do not have GIC observations. Using a test set of transformers we show the potential of this technique to improve GIC modeling capability even where local observations do not exist.

Geoelectric Field and GIC Modeling
Modeling GIC is a useful tool to understand the impacts of geomagnetic storms on power grids. In the current study, we focus entirely on New Zealand, although the fundamental steps undertaken here should be generally applicable to power networks across the globe. We should note, however that the approaches outlined in our study rely upon the availability of significant quantities of GIC observations. We believe that this is not a common 10.1029/2021SW002955 3 of 23 situation for most power networks around the world. The modeling process can be split up into a geophysical step where electric fields across New Zealand are calculated (Vasseur & Weidelt, 1977) and an engineering step which models GIC flowing through every earthed transformer winding and transmission line (Boteler & Pirjola, 2014) in the New Zealand high voltage transmission network.

Ground Conductance Model
A ground conductance model of New Zealand and the surrounding ocean is required to calculate the electric field. We utilize a thin-sheet conductance (TSC) model based on magnetotelluric studies, bathymetry, and geological maps. The TSC modeling technique was developed by Vasseur and Weidelt (1977). The TSC model used in our study to represent New Zealand has been explained before (Divett et al., , 2018 and currently consists of 168 × 168 square grid cells (28°N × 28°W) with a one sixth of a degree (roughly 20 km) grid spacing. Figure 1a shows the thin-sheet conductance map for New Zealand, zoomed in to show the land of New Zealand and nearby oceanic regions. The on-land conductance of each cell represents the integrated conductance of the upper 20 km of the crust. Figure 1b shows the underlying layered structure which consists of four layers, with resistivity of 1,000, 10,000, 100, 1 Ω/m at layer boundaries of 20, 60, and 320 km depths. This is a two-dimensional TSC model with a depth variation specific to the New Zealand geology. However, it should be noted that this representation vastly simplifies the region of active fault lines and tectonic variation.

Geomagnetic Field
A spatially uniform magnetic field equal to that measured by the Eyrewell magnetic observatory is applied to the region of New Zealand. Ideally, we would have liked to have also utilized two variometers, one in Middlemarch in central Otago and another in Te Wharau in the lower North Island. These are the same variometers used in Divett et al. (2020) when looking at the March 2015 geomagnetic disturbance. Unfortunately, due to the recent installation of these variometers, measurements from both sites are only available for two (March 2015 and September 2017) of the 6 significant geomagnetic events we are analyzing and 3 of the 19 smaller events (these 25 4 of 23 events are discussed in Section 3.3). We therefore only use Eyrewell observatory magnetic field observations for all of the geomagnetic disturbance events considered in this study, so we are consistent between all 25 events. Figures 2g-2j in Divett et al. (2020) show how the magnetic field varies across New Zealand during the March 2015 event. Here, we see that the strongest impact is in the lower South Island. Divett et al. (2020) concluded that modeled GICs are more sensitive to variations in the magnetic field than other factors such as land conductivity or network resistances. The various scaling techniques presented in this paper will attempt to account for this lack of magnetic field variation.
Magnetic field data of 1-min period at a resolution of 0.1 nT in the X (positive to geographic north), Y (positive toward east), and Z (positive vertically downward) is provided by the Eyrewell magnetic observatory. For the study presented here, the 1-min X and Y components are used and applied to every grid cell in the model domain.

Geoelectric Field
We calculate the geoelectric field using the thin-sheet model of Vasseur and Weidelt (1977) with the two-dimensional thin-sheet conductance model discussed earlier. Horizontal electric fields are induced at the surface of the Earth by the temporal variations in the northward and eastward magnetic field. As such the magnitude and direction of the magnetic field affects the magnitude and direction of the electric field.
This method is similar to that used in Divett et al. (2020) with one exception: in the current study the range of valid periods are not restricted. In Divett et al. (2020), the valid periods were restricted to those between 5 and 80.5 min by applying a Gaussian filter to those periods outside this range. For the purpose of that study, the period-limited approach was viewed as acceptable as the modeled GIC was being compared with observations that had also been filtered to this period range. In the current study we are trying to reproduce the true measured GIC downsampled to 1-min resolution (due to the 1-min magnetic field input) and so restricting the range of valid periods any more than this is not required. Hence, the range of valid periods for the current study span from 2-min (corresponding to the nyquist frequency) to 1,440 min (24 hr geomagnetic disturbance period modeled). This stretches the numerical limitations of the thin-sheet model but comparison of modeled and observed GIC justifies this approach.

GIC Model
The New Zealand high-voltage AC power network consists of transmission lines with three different voltage ranges: 50 or 66 kV, 110 kV, and 220 kV. Effectively the North and South Islands are isolated as they are only connected by the high voltage DC (HVDC) link, operating at ±350 kV DC which is not represented in the network used in this model. The HVDC link is not affected by GIC as there are no transformers in the circuit. The New Zealand network is made up of a large number of substations each with a varying number of transformers. Some of these transformers are earthed allowing GIC to flow to and from ground through the earthing point of the transformer. These substations are connected together by a network of transmission lines of varying voltages.
We have followed the approach of Lehtinen and Pirjola (1985) for substation-level and the modification in Boteler and Pirjola (2014) for transformer-level to represent New Zealand's power transmission network by voltage levels within each substation connected through transformers and transmission lines of varying resistance and earthed through the Earth grid resistance (EGR). We use this transformer-level network representation to calculate the current through each individual transformer winding and transmission line in the New Zealand network. This allows us to calculate the GIC flowing through a specific transformer winding along with the substation-level GIC through the EGR.
A detailed explanation of how to calculate transformer-level GIC is given in Sections 3.2-3.4 of Divett et al. (2018) and is also discussed in Section 3.2 of Divett et al. (2020). The basic calculation steps are summarized below.
First the current source J m,n , due to GIC is calculated along each transmission line with In Equation 1, we sum the electric field ⃖⃖ ⃗ along the path of the transmission line from node m to n in line segments ⃖⃖⃖ ⃗ . Each segment ⃖⃖⃖ ⃗ is the segment between transmission line pylons. Summing the current source due to all transmission lines connected to that node gives the total current at that node due to the Electric field. We call this the source vector J.
The network admittance matrix, Y, and the earthing impedance matrix, Z, describe the network resistance. The admittance connecting nodes m and n is given by Y m,n = −1/R . The diagonal elements of Y are the sum of the admittance of every transmission line or transformer connected to the node while the off-diagonal elements are the negative admittance of the connections between the two nodes in question. For Z the diagonal elements are the EGR of that node since we assume that all substations are far enough away that the impedance at one node does not effect another node. The substation GIC vector I sub flowing to Earth can be calculated by The GIC through each transformer winding I trans , is calculated with Ohm's law. From the substation GIC we can calculate the node voltage, relative to the local substation Earth as From Equation 3, we can calculate the voltage difference between each node as for all nodes within every substation. We then use the DC resistance of each transformer winding to calculate the GIC flowing through the transformer with = . (5) where m and n are nodes within a substation. As V and V are relative to local Earth, this calculation is only valid for nodes within the same substation. To calculate the GIC flowing through each transmission line, I line we need to add an extra calculation. This is the current source J m,n given in Equation 1, calculated along each transmission line between nodes of different substations. Therefore I line is given as where P line is the number of parallel lines between the two nodes (assuming resistance of each parallel line is the same).

Observations
Measurements of DC currents are available at a large number of substations for multiple transformers in each substation, primarily in locations across the lower and mid-South Island of New Zealand. An extensive description of the New Zealand DC observations are given in Mac  and Rodger et al. (2020) Clilverd et al. (2020). During fast changing measurements the time resolution is at its highest of 4 s for most transformers. Figure 2 shows the locations of substations with DC monitoring equipment and the location of the Eyrewell magnetometer as given by the blue hexagon.

Network Configuration
A transformer level network representation of the New Zealand power transmission network was developed earlier (explained extensively in Divett et al., 2017;Divett et al., 2018Divett et al., , 2020. For the current study, we have extended the network to include the North Island and are using the real transmission line paths for all calculations. The calculations required to model the GIC were previously discussed in Section 2.4. The New Zealand GIC model requires an understanding of four key network characteristics. These include;  Table 1 shows the number of substations, transformers, and transmission lines throughout the New Zealand electric network on the high voltage side utilized in the GIC model (66 kV or greater). The term "upper phase" transformers refers to the series winding on auto transformers between two nonzero kV buses. For example, the winding between the 110 and 220 kV buses on transformer #4 at Halfway Bush (HWBT4) is termed the "upper phase," while the common winding from 0 to 110 kV is termed an "earthed" transformer. It is also worth noting that those transformers referred to as "unearthed," are only unearthed on the high voltage side (66 kV or greater). These transformers are earthed on the low voltage (local) network. Past observations have shown that GIC impact from transformers earthed on the low voltage side is minimal (NERC, 2013). Low voltage networks have high line resistances and short line lengths when compared to the high voltage transmission lines, leading to less GIC and less stress on these transformers. It is worth noting that although only 40% of the New Zealand transmission lines are at 220 kV, they make up 80% of the GIC flowing through transmission lines.

Selection of Geomagnetic Disturbance Events
The New Zealand DC observations of GIC span about 20 years. Throughout this period there have been multiple geomagnetic events that have generated significant levels of GIC measured at various transformers. For our current study we have selected 6 significant disturbance events, each covering a 24 hr period, which are listed below: From this point forward each of these events will be referred to by their month and year, that is, the first event will be referenced as November 2001. A 24 hr period was selected to account for the various duration of each geomagnetic disturbance and to ensure the whole event is contained within the time window. These six events represent some of the largest GIC magnitudes measured at transformers in New Zealand. Some properties of these events are listed in Table 2. Note also that five of the six events (all but September 2017) appeared in the empirical New Zealand extreme storm GIC study investigating peak GIC magnitudes during the largest 25 geomagnetic disturbances from 2001 to 2015 .
The November 2001 geomagnetic disturbance is the most significant | ′ | value observed by the Eyrewell magnetic observatory from 2001 to 2020. The GIC observed during this event has been extensively discussed and analyzed in earlier studies (Béland & Small, 2004;Mac Manus et al., 2017;Marshall et al., 2012). The observed maximum GIC magnitude of 33.1 A was recorded at the Islington #6 transformer (ISLT6). A similar sized geomagnetic disturbance has yet to affect New Zealand in the years since, although we are confident larger events have occurred in the past, and will also happen in Note that the location of the maximum GIC reported varies across the events, being ISL T6 for the first two, and HWB T4 for the last 4. The maximum GIC is also prior to being down sampled to 1-min, hence it is larger than that given in Table S9 in the Supporting Information S1 Table 2 Properties of the 6 Significant Geomagnetic Disturbances the future. In contrast to November 2001, the September 2017 geomagnetic disturbance recorded the lowest | ′ | value of the 6 events analyzed, but was associated with the largest observed GIC magnitude measurement in New Zealand of 48.9 A measured at the Halfway Bush #4 transformer (HWBT4). The September 2017 disturbance has been examined earlier by Clilverd et al. (2018) who looked at long-lasting GIC and harmonic distortion observed in Dunedin, (the city hosting HWBT4) as well as by Rodger et al. (2020) who focused on harmonic distortions across the entire countries electrical network. Both of these studies utilized harmonic measurements also carried out by Transpower NZ Ltd. Unfortunately there are no GIC observations from ISLT6 during the September 2017 event and HWBT4 measurements only began in October 2012. Rodger et al. (2017) looked at a number of events with measurements at both transformers and found that the GIC is consistently approximately a factor of three times higher at HWBT4 than ISLT6. While measurements at HWBT4 were not being undertaken in November 2001, the transformer was damaged and replaced as a result of this event. It is believed that this transformer experienced ∼100 A of GIC at this time .
On top of the 6 geomagnetic disturbances we have also looked at 19 smaller geomagnetic disturbance events, requiring that they must have a maximum GIC magnitude of <15 A at every transformer, but >1 A for at least one transformer. To do this, one day each year from 2001 to 2019 was randomly selected that fits this criteria. As most days are very geomagnetically quiet it is highly likely that a randomly selected day out of any given year would have essentially zero GIC across all the measured transformers. As we want to include small geomagnetic disturbance events (i.e., not totally quiet) a randomly selected day that does not fit the two criteria above are discarded and a new one randomly selected. For all of these small disturbance events the full 24 hr period of the day selected is also modeled. The properties of these smaller events are given in Table 3. We add the smaller disturbance time periods to increase our data set of GIC activity, but also to span a wider range of power spectra (i.e., from these smaller events through to the larger disturbances). These events will be useful when creating average power spectral scaling curves in Section 5.1.

Modeling of Past Events
For all 25 events in Tables 2 and 3 the GIC was modeled at 1-min resolution for a 24 hr period at every earthed transformer and flowing through every transmission line. As the modeling relies on the thin sheet technique with inbuilt assumptions, it is not meaningful to undertake higher time resolution GIC modeling (see the discussion in Divett et al., 2020). As noted above, the measured GIC has a varying resolution (as high as 4 s during geomagnetic storms). For consistency when comparing the GIC observations with the modeled GIC and to perform Fourier transforms, we have down sampled the measured GIC to a constant sample period of 1 min. This was carried out using a Gaussian filter with a 90-s window centered around 1-min intervals. It should be noted that this is essentially the same approach as how the high time resolution magnetometer data is converted to 1-min values (St-Louis, 2020).
Throughout the paper the September 2017 geomagnetic disturbance is discussed and specific examples from some transformers are given. The full observations and modeling output for all six significant geomagnetic events and transformers can be found in the Tables in the Supporting Information S1. Those tables include Pearson correlation coefficients between the measurements and modeled GIC (S1-S4), Performance parameter (P) between the measurements and modeled GIC (S5-S8), absolute maximum GIC values (S9-S13) and absolute mean GIC values (S14-S18) throughout the 24 hr time period for each event analyzed. Results taking into account all 6 significant geomagnetic events are given throughout the paper and references to these tables are made. In this paper, the Performance parameter is defined as given in Equation 1 of Marsal and Torta (2019) and shown below.
Here, RMSE is the root mean square error between the observations and model, while σ o is the standard deviation of the observations.

September 2017
We will initially discuss the September 2017 geomagnetic disturbance in detail to demonstrate some of the issues when comparing the modeled and observed GIC data.
The spectra of the electric fields and GICs were modeled for the full 24 hr duration. An inverse Fourier transform is then applied to produce a time series of the transformer level GIC through each transformer winding. The top panel of Figure 3 shows a comparison between the absolute mean measured and modeled GICs for the September 2017 geomagnetic disturbance, for all transformers where observations are present. In red, we show the mean (average) absolute GICs modeled across the 24 hr period, while in blue is the mean absolute measured GIC at transformers that were monitored throughout the event. If a transformer has no measurements then a black cross is indicated at zero GIC. The dotted black lines separate the substations so dots within these lines represent GIC through a transformer at the same substation which are ordered approximately from north to south.
As seen in the top panel of Figure 3, for most transformers (51 of 66 with measurements) we find the mean modeled GIC (red circles) are equal or larger than the measured GIC (blue circles) for the September 2017 storm. In contrast, the bottom panel compares the maximum modeled and measured GICs. Here, only 28 of the 66 transformers have modeled GICs that have larger or equal maximum values when compared with the measured GIC.

Extension to Multiple Storms
We find that across all the six large geomagnetic disturbances the relationship between the modeled GIC at each transformer is consistent.  However some transformers have essentially the same GIC measured for all 6 events. We see that for the substations of Timaru, Tekapo B, Ohau A, and Ohau B the mean modeled GIC values are all ∼0.6 A when averaged across all six large storms.
Consistent relationships are also found when looking at the maximum modeled GIC. The full comparison of the mean and maximum modeled GIC for each transformer and large geomagnetic disturbance is shown in Tables S15 (mean) and S10 (maximum) of the Supporting Information S1. This consistency is expected, as the model uses the same network configuration, resistance, and ground conductivity for all six events. The only difference comes from the changing magnetic field variation, provided by the Eyrewell observatory, applied across the entire modeling domain.
It should be noted that there is much less consistency from storm to storm, in the GICs measured at each transformer location compared to modeled GICs. This can be seen in Tables S14 and S9 in Supporting Information S1, which provide the mean and maximum measured GIC (respectively) for each large geomagnetic disturbance and transformer. Unlike the model output, the observations made during the November 2001 geomagnetic disturbance do not show a consistent mean GIC double that of the September 2017 event. Across all transformers the mean measured GIC is 0.3 A for these two events however for individual transformers the relationship varies with some being lower for November 2001 while others are larger. There are however some transformers that show consistent mean measured GIC across all 6 geomagnetic disturbances. Transformers at the substations of Timaru, Tekapo B, Ohau A, and Ohau B have mean measured GIC values of ∼0.2 A. Clearly our initial modeling approach is not capturing the level of variability in the observed GIC. This is probably due to a range of factors, potentially including the use of a nonspatially varying magnetic field, changes in the transmission network or the varying range of frequencies driving GIC. Together, these factors could lead to these differences in mean GIC from disturbance to disturbance. In reality, the magnetic field disturbance will vary across the country. However, we have applied the Eyrewell magnetic observatory data across our entire modeling domain due to a lack of information from a spatially distributed magnetometer network.
We feel it is important to point out that some locations do not show the largely consistent pattern with the maximum and mean GIC being larger in the modeling than in the measurements. For example, some substations in the center of the South Island (i.e., Benmore, Aviemore, Waitaki) show near-zero mean GIC values for both the measured and modeled GIC (as seen in Figure 3 and Tables S14 and S15 in Supporting Information S1). We also see that transformers at Halfway Bush and South Dunedin have the largest maximum and mean GIC, both observed and modeled, when contrasted with the rest of the substations which are monitored (as seen in Figure 3 and also Tables S9, S10, S14, and S15 in Supporting Information S1). This suggests the model is approximately correct, consistent with the conclusions of Divett et al. (2020).
We now consider the average of all monitored transformers taken as a whole across all events. This information can be found in the last lines of Tables S9 to S18 in the Supporting Information S1. The mean modeled GIC, averaged across all monitored transformers and all storms is 0.4 A (Table S15 in Supporting Information S1). The maximum modeled GIC, averaged in the same way, is 6.4 A (Table S10 in Supporting Information S1). In contrast the mean observed GIC, averaged across all monitored transformers and storms is 0.3 A (Table S14 in Supporting Information S1). The maximum observed GIC, averaged in the same way is 5.0 A (Table S9 in Supporting Information S1). In both cases, the ratio between the observed value and modeled value is ∼1.3. We find that typically the model predicts the maximum GIC with a similar accuracy to the mean GIC modeled. In addition, this confirms that typically the model GIC is slightly higher than what is observed.
While this result applies for the majority of monitored transformers, it is strongly influenced by those transformers with rather low currents (a few amps, as is common even amongst these fairly large disturbances). If we look at transformers with large GIC values, for example, HWBT4, HWBT6, and ISLT6 the ratio of maximum modeled GIC to measured GIC is ∼0.4, but for the mean GIC this ratio is ∼0.7. This highlights some of the difficulties the model has in reproducing the larger GIC peaks. For these highly responsive transformers the model GIC is actually lower than the measured values. That identifies a difficulty in using the existing modeling approach for hazard identification.
We find that generally the thin sheet modeled GICs are typically adequate for the 6 significant geomagnetic storms, which is to say that for most transformers the modeled maximum GIC is typically similar to that observed. The detailed evidence of this can be found in Tables S9 and S10 in the Supporting Information S1. However for transformers with a maximum GIC larger than 10 A the model typically does not reproduce those large values. Throughout the 6 significant geomagnetic disturbances there are 37 occurrences of a transformer with a maximum observed GIC larger than 10 A. Of these the model GIC is lower for 25 (67%) of them, on average by ∼8 A (30%).
Unfortunately, as shown, for some transformers the model does a rather poor job of predicting the magnitude of the expected GIC. It is unclear if this is due to fundamental limitations in the thin-sheet code, differences between the conductance model employed and reality, inductance in the electrical network, small scale structure or the magnetic field forcing, or perhaps other as yet identified issues. In the next section, we present a technique to try and compensate for the identified weakness for locations with large measured GIC magnitudes.

Developing Corrected Power Spectra
As noted earlier, our study is part of a long term effort working toward a validated model by which we can predict GIC magnitudes across New Zealand during an extreme geomagnetic disturbance. This is necessary for the New Zealand power grid operator to plan for an extreme event and develop mitigation approaches. As reported in the previous section, the "basic" thin sheet modeling approach described by Divett et al. (2020) does a reasonable job of describing the time variation of GIC, but tends to produce GIC magnitudes which are typically too low at the most responsive transformers. However, as we will demonstrate below for a given geomagnetic storm we can use the modeled power spectrum of the GIC through a given transformer and scale it up to a similar magnitude to the power spectrum determined from the measured GIC. This scaling process allows us to empirically account for some inabilities of the thin sheet to accurately model particular frequencies of GIC, and should also correct other systematic modeling inaccuracies (as outlined in the previous section).
As an example we will consider the HWBT4 power spectrum during the September 2017 event. The Halfway Bush substation is located only ∼5 km away from the South Dunedin substation. However, in contrast to SDNT2, the thin sheet model fails to model the GIC at HWBT4 correctly (as earlier discussed in Divett et al., 2020, and also shown in Figure 3). We can attempt to correct this by scaling the power spectrum as described below, and as shown in the top panel of Figure 5. The power spectrum is corrected using a third-order polynomial fit. We fitted polynomial curves to the modeled and measured power spectra independently. We then find the difference between the polynomial fit for the modeled GIC and the polynomial fit for the measured GIC and shift the model power spectrum by this difference. In the upper panel of Figure 5, the black line represents the scaled model GIC power spectra that has been corrected to produce a similar magnitude as the measured GIC power spectra, which we term the "Corrected power spectra".
Using this scaled power spectra and performing an inverse Fourier transform produces the lower panel of Figure 5 which presents GIC at HWBT4 in the time domain. The corrected GIC model output is shown in black which we refer to as the "scaled GIC". This scaling approach is significantly better at modeling the long period of elevated GIC that occurred from 12:00 to 13:00 UT. Unfortunately, the initial few hours are significantly worse, indicating that this scaling method is still missing some important details. Overall the Pearson correlation coefficient has increased from r = 0.776 to r = 0.816. Although the magnitude in the initial few hours is overestimated with the scaled GIC, the overall shape of the time series has been improved, leading to this better Pearson correlation coefficient. This is further supported with the improved performance parameter from P = 0.216 to P = 0.443. The maximum measured GIC is 47.5 A, significantly larger than the initial maximum modeled GIC of 9.1 A. However using the corrected GIC power spectrum to create the scaled GIC gives an improved modeled maximum GIC of 40.8 A, clearly significantly closer to the measured GIC.
As for the rest of the transformers with measured GIC, the top panel of Figure 6 displays the change in the Pearson correlation coefficient when comparing the measured GIC with the scaled model GIC as opposed to the non-scaled modeled GIC. This information is also given in Tables S1 and S2 in the Supporting Information S1 For example, for HWBT4 the increase in Pearson correlation coefficient is 0.04 (r = 0.776 to r = 0.816). This panel shows that after scaling only one transformer (OHBT8) shows a significant decrease in Pearson correlation coefficient, which we take to be a decrease of more than 0.05, and only five transformers decrease, while 59 have increased Pearson coefficients. Similar results are found for the other 5 significant geomagnetic disturbances we have analyzed (as shown in Table S1 in Supporting Information S1 for the model GIC and Table S2 in Supporting Information S1 for the scaled GIC in the supplementary material). All the 5 transformers which show decreases in Pearson correlation coefficients have maximum GIC (both modeled and measured outputs) that are less than 2 A. For such small values, the extreme storm risk is clearly very low, and hence we are not concerned by the decrease in correlation for those transformers.
Likewise, the performance parameter for the September 2017 event also improves. Of the 64 transformers with measurements, 45 (mean P = 0.347) show positive performance parameters compared to 34 (mean P = 0.248) prior to scaling. Tables S5 and S6 in the Supporting Information S1 show similar results for the other 5 significant geomagnetic disturbances.
With regards to possible transformer damage during a geomagnetic storm, it is large GIC that is of concern, in particular the length of time for which the currents are large. While the Pearson correlation coefficient is a good indicator of the ability to model the overall GIC time series variation, it does not take into account the magnitude of these variations. In the bottom panel of Figure 6, we show the maximum absolute GIC for the transformers with GIC measurements available during the September 2017 event. With the red circles, we compare the measured GIC with the modeled GIC, both at 1-min resolution, while the black circles compare the scaled GIC against the measured quantities. Clearly the scaled GIC provides significantly better quality matches to the maximum absolute measured GIC as it is a closer fit to the perfect fit given by the dashed black line. This Figure demonstrates that in most cases the application of the scaling technique accurately predicts the peak GIC during the geomagnetic disturbance. Similar improvements are seen for the other 5 significant geomagnetic disturbances that we have modeled in this study (as shown in Tables S11 in the Supporting Information S1).
We now test the application of the corrected power spectra to smaller geomagnetic disturbances, to confirm the approach is valid over a range of disturbance sizes. Here we examine 22 April 2017 for HWBT4. Again the modeling is improved by applying the scaled power spectra, as shown in Figure 7. In this Figure, it is clear by eye that the scaled modeled GIC in the bottom panel has much better correlation with the measured GIC than is present in the upper panel. Here the Pearson correlation coefficient has increased from r = 0.781 (top panel) to r = 0.853 (bottom panel). In this Figure the maximum absolute GIC is also improved. The modeled scaled GIC value increases to 10.9 A compared to 3.7 A determined from the original modeling approach. As such it is significantly closer to the measured 1-min time resolution peak GIC of 12.3 A. Although we are more interested in modeling time periods of very large GIC and accurately predicting the peak GIC during very large geomagnetic disturbances, this figure shows that the model and correction technique used are still effective at capturing both the time variation and magnitude of smaller GIC values (GIC < 10 A) which occur throughout the time period considered. The improvement in the Performance parameter from P = 0.325 to P = 0.543 also supports this. The dashed black line shows the 1:1 fit between the modeled (y-axis) and measured (x-axis).

Scaling Transformers for Future Storms
One of our long term goals is to extend the GIC modeling across the entire New Zealand electrical transmission network, to allow identification of potentially "at risk" transformers during very large and extreme geomagnetic disturbances. The scaling technique described earlier relies on comparison and scaling of modeled results with information from GIC observations, which are limited to a subset of transformers that have observations for the event in question (see Table 1 in Rodger et al., 2020). The technique described in the previous section produces corrected power spectra separately for each disturbance. Hence we seek an approach to generalize that technique for all disturbances, which we will then extend again to un-monitored transformers. That is outlined in the current Sections 5.2 and 5.3 below.
In Figure 8, we show the corrected power spectra for all geomagnetic disturbances (including both Tables 2 and 3 events) where observations are available from HWBT4. The events are colored to identify geomagnetic disturbance periods from both the larger and smaller events. Recall that measurements are not always available for all 25 disturbances, hence Figure 8 shows only 3 corrected power spectra for the large storms, and 5 for the smaller disturbances. From the eight corrected power spectra, we then produce an average of all these power spectra, shown in Figure 8 as the black line. We will term this the "local multi-storm corrected power spectra" (LMSC power spectra) that can be utilized to represent the response at that specific transformer for any geomagnetic disturbance. We have followed the same approach to produce a LMSC power spectra for all 73 transformers for which there are GIC observations.
It is worth noting that over the 19-year data set the New Zealand electrical power network has changed. New transmission lines have been introduced and transformers have been both commissioned and decommissioned.
For the purpose of this study and to enable the creation of the various LMSC power spectra we have left the modeled network unchanged in the state it was at the start of 2018. This is useful in that the modern network state is most representative of the network which would be impacted by a future extreme storm. However, the approximation is also necessary as we do not have details for all of the changes which have occurred progressively over the two decades considered.  (Table 3). In blue is the measured GIC which is compared against the modeled GIC (red in top panel) and the output of the scaled model GIC (black in bottom panel).
In discussing Figure 5 (showing GIC at HWBT4 in September 2017), we noted that for most times the model GIC output values were improved after application of a corrected power spectra. We can now apply the LMSC power spectra, shown for HWBT4 in Figure 8, and determine the new GIC output produced with that LMSC power spectra. The results of these calculations are shown in Figure 9.
In Figure 9 we see that when using LMSC power spectra (magenta line), the new model GIC output is only slightly smaller than was seen in Figure 5 (the black line from Figure 9). Using the corrected power spectra that is unique for HWBT4 during the September 2017 event produces a maximum GIC of 40.8 A, rather close to the measured GIC maximum of 47.5 A. When we apply the more versatile LMSC power spectra correction the maximum modeled GIC is 32.1 A, a clear decrease, but significantly higher than for the basic modeling approach of 9.1 A. We also find similar Pearson correlation coefficients of 0.816 for the corrected power spectra and 0.822 for the LMSC power spectra correction. The smaller maximum GIC when using the LMSC power spectra is an acceptable compromise as the LMSC approach can be used to correct multiple different geomagnetic disturbances.

Scaling Transformers With No Measurements
As mentioned in Table 1, there are over 250 earthed transformers in the New Zealand electrical transmission network, all which could have potentially large GIC during geomagnetic disturbances. Around 75% of these transformers have no GIC measurements to give a direct comparison with the model GIC. It seems likely that some of these unmonitored transformers could be at risk due to an extreme geomagnetic disturbance.
In an attempt to create a way to accurately model the GIC at sites without existing GIC measurements, we have modeled a large number of geomagnetic disturbance events to determine a nationwide multi-storm corrected power spectra (NMSC power spectra) by averaging all the LMSC power spectra. This NMSC power spectra allows us to model all transformers across the New Zealand transmission network using the assumption that the un-monitored transformers will be reasonably well represented by the transformers for which there are measurements. Figure 8. Power spectrum of the modeled 1-min geomagnetically induced currents (GIC) fitted curves at HWBT4 during the events that have measured GIC. Blue lines are fits from the large significant Geomagnetic storms while the red lines are from the smaller disturbances. The mean power spectrum for all these curves is given by the black line. This is termed the local multistorm corrected power spectra (LMSC) power spectra.
The top panel of Figure 10 shows the LMSC power spectra for each of the transformers with measured GIC.
There are a total of 73 transformers for which we have GIC measurements and hence 73 LMSC power spectra are shown by the black lines in Figure 10 (top panel). From this we have determined the mean polynomial fit, producing the NMSC power spectra (blue line). The shaded red region represents the 1-sigma NMSC fit based on the various 1-sigma fits for the LMSCs. Figure 10 shows a large variation in the LMSC power spectra across all the transformers. For the purpose of modeling those transformers without any GIC measurements, the NMSC power spectra should help account for many limitations in modeling particular frequencies. While the variation in the LMSC are very large, the impact on calculated maximum and mean modeled GIC is less extreme. The maximum and mean GIC modeled using the NMSC power spectra are relatively close to those produced using the LMSC power spectra, and they are relatively close to observations. By looking at the maximum and mean values for the LMSC (in Tables S12, S17 in Supporting Information S1) and NMSC (in Tables S13, S18 in Supporting Information S1) corrections we can see that on average across all transformers the LMSC and NMSC values are similar. The maximum GIC for the LMSC scaled GIC averaged across all six significant geomagnetic disturbances is 4.1 A compared to 4.6 A for the NMSC. Likewise for the mean GIC, both the LMSC and NMSC produce 0.3 A across the same 6 disturbance events.
The bottom panel in Figure 10 compares the various maximum model GIC against the measured GIC. As expected, using the corrected power spectra that is unique to each transformer and geomagnetic disturbance (black) produces the closest maximum GIC to the measurements (average difference of 0.7 A), while the original noncorrected (red) model GIC is the worst (average difference of 3.4 A). The LMSC (magenta) and NMSC (green) power spectra scaled GIC have average deviations of 1.8 and 2.1 A. For values that better represent all transformers (e.g., below 15 A) the LMSC and NMSC models produce very similar results (as mentioned earlier in the text). The lower panel shows three transformers with significantly larger measured GIC. These are HWBT4 (47.5 A), HWBT6 (39.2 A), and SDNT2 (33.5 A). For the two at Halfway Bush (HWBT4, HWBT6) the NMSC Figure 9. Geomagnetically induced currents (GIC) at HWBT4 during the September 2017 event. In black is the scaled GIC using the corrected power spectra while in magenta is the scaled GIC using the local multistorm corrected power spectra (LMSC) power spectra. model produces significantly less (and further from the measured) maximum GIC than the LMSC model. For these locations we would continue to use the LMSC model.
In the lower panel of Figure 10, we show the upper and lower 1-sigma regions for the LMSC and NMSC model GIC. This involves calculating the corresponding LMSC for the upper and lower 1-sigma errors and using those fits for the corresponding 1-sigma NMSC fits. This results in a range of ∼±8% for the LMSC and ±12% for the NMSC.
The Redclyffe substation located on the east coast of the North Island near the city of Napier is one such location with no measured GIC. Figure 11 shows the difference between the original modeled GIC at transformer #1 (RDFT1) and the modeled GIC once the NMSC power spectra has been applied. The increase in GIC in this example is significant, particularly for the largest GIC values. The peak GIC increased more than 50% from 26 to 41 A. We have selected this transformer as it has one of the largest modeled GIC for a transformer without measurements. Rodger et al. (2020) showed a correlation between even harmonics and large GIC at sites with measured GIC in the South Island. That work also showed that Redclyffe had some of the largest even harmonics in the North Island suggesting that the modeling approach for nonmonitored transformers does improve the quality of the outputs. Mukhtar et al. (2020) discussed the Redclyffe substation while looking at using MT data to model GIC in the North Island. In that study they mentioned that Redclyffe has similarities with the Islington substation in terms of their physical location and electrical setup.
It is worth highlighting that when considering the North Island, particular the northern part, using magnetic field from Eyrewell is likely to be relatively over representing the true magnetic field. Intuitively, the electric fields and GIC are likely to be smaller than the South Island. However the NMSC is also itself an average which is lower than the spectrum that might be expected from a larger storm leading to an underestimate in the corrected GIC. It is unknown if these two effects would account for each other leading to a rather accurate GIC modeled with the NMSC. For future studies looking at extreme events a range of NMSC fits would need to be considered to produce a range of expected GIC. Figure 11. Modeled geomagnetically induced currents (GIC) for RDFT1 during the September 2017 geomagnetic event. In red is the original modeled GIC while the scale GIC output from the nationwide multi-storm corrected power spectra (NMSC) power spectra is given in green.
As shown in Figure 11, the NMSC power spectra can be applied to transformers without measurements. To further validate the use of the NMSC we compared the output from the NMSC approach for transformers with measurements to show that it still provides an improvement over the original model. To undertake this test, we excluded the ISLT3 LMSC when creating the NMSC and then determined the modeled ISLT3 GIC using the test NMSC. The result of this calculation is given in Figure 12. The Pearson correlation coefficients are improved (0.830 against 0.886) as is the performance parameter (0.316 against 0.567) and the Figure shows that the NMSC approach provides significant improvements over the original unscaled model GIC particularly at the largest peaks and is a valid approach to scaling transformers without GIC measurements.

Summary and Conclusions
Using a thin sheet model of the New Zealand electrical power network we model the GIC during 25 disturbed periods of 24 hr duration, ranging from large geomagnetic storms to weakly active periods. Initial analysis is undertaken of GIC during the September 2017 geomagnetic storm, and a technique of adjusting the model GIC output spectra in the frequency domain is applied. The power spectra of the thin sheet model output is compared against the measured GIC power at each of 73 sites where in-situ measurements are made, and a local scaling factor determined that corrects the model output. This technique improves the correlations between maximum model and observed GIC by between 10% and 40% depending on the transformer. Average local scaling factors were then found for all 25 disturbed periods analyzed giving a characteristic frequency domain scaling factor for each transformer.
This technique demonstrates that the GIC produced from the thin sheet model for New Zealand can be scaled in the frequency domain to produce a more accurate comparison with the measured GIC. We have shown how corresponding local scale fits of the power spectra can be created for each unique transformer. Once local scaling factors are determined, they can be averaged to describe a scaling factor for the network as a whole. For transformers with large GIC we find that the locally scaled model output provides a more accurate agreement with the maximum measured GIC. For smaller measured GIC values (below 15 A) the local and national scaled GIC results are very similar. The use of national scaling factors allow more accurate GIC modeling for transformers without in-situ measurements. We provide an example of a North Island transformer, Redclyffe, previously found to generate harmonic distortion. Here, the corrected thin sheet model peak GIC output increases by more than 50%, and for the September 2017 geomagnetic storm, maximum GIC of 41 A are calculated.

Data Availability Statement
The New Zealand electrical transmission network's DC characteristics and DC measurements were provided to us by Transpower New Zealand with caveats and restrictions. This includes requirements of permission before all publications and presentations and no ability to provide the observations themselves. In addition, we are unable to provide the New Zealand network characteristics due to commercial sensitivity. Requests for access to these characteristics and the DC measurements need to be made to Transpower New Zealand.