Geomagnetically Induced Current Model Validation From New Zealand's South Island

Geomagnetically induced currents (GICs) during a space weather event have previously caused transformer damage in New Zealand. During the 2015 St. Patrick's Day Storm, Transpower NZ Ltd has reliable GIC measurements at 23 different transformers across New Zealand's South Island. These observed GICs show large variability, spatially and within a substation. We compare these GICs with those calculated from a modeled geolectric field using a network model of the transmission network with industry‐provided line, earthing, and transformer resistances. We calculate the modeled geoelectric field from the spectra of magnetic field variations interpolated from measurements during this storm and ground conductance using a thin‐sheet model. Modeled and observed GIC spectra are similar, and coherence exceeds the 95% confidence threshold, for most valid frequencies at 18 of the 23 transformers. Sensitivity analysis shows that modeled GICs are most sensitive to variation in magnetic field input, followed by the variation in land conductivity. The assumption that transmission lines follow straight lines or getting the network resistances exactly right is less significant. Comparing modeled and measured GIC time series highlights that this modeling approach is useful for reconstructing the timing, duration, and relative magnitude of GIC peaks during sudden commencement and substorms. However, the model significantly underestimates the magnitude of these peaks, even for a transformer with good spectral match. This is because of the limited range of frequencies for which the thin‐sheet model is valid and severely limits the usefulness of this modeling approach for accurate prediction of peak GICs.


Introduction
Transformers in high-voltage electrical transmission networks have been damaged by geomagnetically induced currents (GICs) during space weather events. During these events, GICs are induced in transmission lines by an electromotive force along each line as a result of geoelectric field variations associated with changes in the geomagnetic field. The geomagnetic field varies in response to the interaction of solar wind disturbances with the Earth's ionosphere and magnetospheres. These physical processes are described in relevant textbooks (e.g., Bothmer & Daglis, 2007). GICs have been reported in many low-latitude to midlatitude countries including the United Kingdom (Erinmez et al., 2002), South Africa (Gaunt & Coetzee, 2007), Brazil (Trivedi et al., 2007), China (Liu et al., 2009), Spain (Torta et al., 2012), Australia (Marshall et al., 2013), and New Zealand. In New Zealand, a transformer at Dunedin's Halfway Bush (HWB) substation was written off after it was damaged during the November 2001 storm (Béland & Small, 2004;Marshall et al., 2012;Mac Manus et al., 2017). The location of Dunedin is shown in Figure 1.
Modeling GICs flowing through transformers in the high-voltage transmission network is needed to increase understanding of the impact of GICs on transformers within the grid and allow improved hazard planning. Models to calculate GICs in the transmission network require two steps: (i) modeling the geoelectric field due to the combined effects of the magnetic field variation during a storm and the varying ground conductance and (ii) using a network model to calculate the GICs flowing through each transformer due to the geoelectric field. The thin-sheet modeling approach of Vasseur and Weidelt (1977) (VW77 hereafter) has been successfully used to calculate the geoelectric field in the United Kingdom (e.g., Beggan, 2015;Beggan et al., 2013;Mckay, 2003), New Zealand , and Austria (Bailey et al., 2017(Bailey et al., , 2018. However, these studies have only simulated the response of the Earth to a magnetic field of a single frequency. In the present study, we use the same VW77 thin-sheet modeling approach but use a range of frequencies, limited to those for which VW77's approach is valid. GICs have been calculated from the electric field using the matrix method of Lehtinen and Pirjola (1985) (LP85 hereafter) by Koen and Gaunt (2003), Torta et al. (2012), Beggan et al. (2013), Blake et al. (2016), Bailey et al. (2017), and Divett et al. (2017). These studies calculated substation-level GIC by assuming that each substation within the electrical network can be approximated as a single DC resistance to Earth. In reality, most substations contain several transformers connected to multiple voltage levels. Often, each transformer has a different DC resistance. Boteler and Pirjola (2014) showed how autotransformers and normal transformers connected in a simple test substation can be represented within LP85's nodal network model. Boteler and Pirjola (2017) further showed how to use this representation to calculate the transformer level current in an extension of the LP85 approach. Divett et al. (2018) used such a representation to develop a transformer level network model for the transmission network of the South Island of New Zealand and calculate the transformer level GICs for an idealized magnetic field variation. In the present study, we use this network model to calculate transformer level GICs during the St. Patrick's Day Storm of 17 March 2015, with the goal of model validation.
Once a network model has been developed, the final step toward validation requires comparison with measurements of GICs at transformers in the network. It is particularly important to compare against a broad selection of substations in the network because there can be substantial variation in the GIC flowing through different transformers at different locations in the network and even within a substation (Divett et al., 2018;Mac Manus et al., 2017;Rodger et al., 2017). In previous model-measurement comparisons, GIC measurements have been sparse, limiting the validation to a small number of transformers, in substations isolated from one another. Bailey et al. (2018) validated their model by comparing against observed GIC at three transformers in three different substations in Austria. They used a similar VW77 and LP85 approach to model GICs to that which we apply in the present study. However, they simulated the response to a single frequency of the rate of change of the magnetic field variation, assuming that the energy of the whole spectrum of frequencies responds in the same way as this single frequency. Butala et al. (2017) had access to measurements of GICs at 23 transformers in Wisconsin, USA but compared modeled GICs against only five of these transformers (located in four substations) due to poor data quality at the other locations. Their modeling approach used a transfer function approach to calculate electric fields and the commercial package PowerWorld to calculate transformer level GICs. Nakamura et al. (2018) compared measured GICs at two substations in the Tokyo area against GICs modeled using a Finite Difference Time Domain approach to calculate the electric fields and a substation-level network model to calculate GICs. However, they did not know how many transformers were at each substation, so they assumed 10 and divided the substation-level GICs by 10. This assumption limited their validation against measurements to two transformers as all transformers within a substation are effectively identical. In New Zealand, Transpower has measured GICs at up to 58 individual transformers since 2001 , providing a particularly large data set to compare against modeled GICs. This is the data set we have considered in the present study.
In this study, we start by describing the magnetic field variations and measured GICs for the St. Patrick's Day Storm of 17 March 2015 in section 2. We use the same technique as (Divett et al., , 2018 to model the electric fields in the New Zealand region and GICs in New Zealand's South Island, which we will describe in section 3. We also rely upon the same ground conductance model as those earlier studies but with a larger domain for the numerical grid. This increases the range of frequencies that the model is valid for, under VW77's modeling assumptions which are described in section 3.1.2. We calculate the spectrum of the geoelectric field for the full range of valid frequencies, by applying the Fourier coefficients of the spectrum of the magnetic field variation, for each valid frequency, to the VW77 model. This technique is a significant deviation from previous studies (e.g., Beggan et al., 2013;Blake, 2017;Bailey et al., 2017;Mckay, 2003), who used the same VW77 and LP85 models to calculate GICs. Those studies ran the simulations for a single frequency, by applying the magnitude (or in some cases the rate of change) of magnetic field variation, at each point in time, to the VW77 code. In the present study, we compare the modeled GICs with observations in the modeled frequency range and in the time domain in section 4. Finally, we test the sensitivity of our modeled GICs to variations in the inputs in section 5. We believe our study is the most extensive attempt to date of GIC model validation using a real storm.

Observations During the St. Patrick's Day Storm of 2015
We undertake model comparisons with an extensive set of GIC observations made at multiple locations in the lower half of New Zealand's South Island. The geoelectric field modeling relies upon geomagnetic field changes interpolated from magnetic field observations made by magnetometers and variometers in the region.

Magnetic Field Measurements
In order to model GICs during the St. Patrick's Day Storm, we utilized the magnetic field observations from four sites during this event. These measurements allowed the generation of a spatial distribution of the magnetic field across the South Island. These magnetic field observations were made at the four locations shown in Figure 2e. Eyrewell (EYR) and Macquarie Island (MCQ) are DI fluxgate magnetometers. Both instruments are part of INTERMAGNET, a global network of geomagnetic observatories, established to monitor the Earth's magnetic field around the world and provide rapid magnetic observatory data exchange between the international scientific community. The variometers at Middlemarch (MDM) and Te Wharau (TEW) are operated by Osaka Electro-Communication University, Japan and are part of the CRUX array. These instruments were described by Obana et al. (2015). The MDM data have previously been used in New Zealand GIC research by Clilverd et al. (2018), who studied the September 2017 storm.
During the geomagnetic storm which occurred on 17 March 2015, the global Kp index reached a maximum of 8-, while the maximum local EYR K index was 6. This storm resulted in significant geomagnetic variations at the EYR observatory, beginning with a sudden commencement at 04:46 UT, as shown in Figures 2a and 2b. At this time, a comparatively large horizontal rate of change (H′) of 68 nT/min was observed (using 1-min resolution data). This is the 13th largest H′ value reported in New Zealand since 1 January 2001 . The magnetic field variations at MDM and TEW show similar behavior to that at EYR over the day. However, the variations at TEW are smaller than EYR, and those at MDM are larger than EYR, as expected due to their relative latitudes. Due to MCQ's location in the auroral zone, the variations at MCQ are considerably larger than the other three locations with the largest variations in B x being in the opposite sense to the other locations. This difference in sense, for example, around 9:00 UT, would appear to be the effect of the auroral electrojet moving into a location between MCQ and MDM, so B x is enhanced on one side (positive) and depressed on the other (negative). A substorm current wedge typically shows such variation.
The four magnetometers used here all provide 1-min magnetic field data at a resolution of 0.1 nT in the X (positive to geographic north), Y (positive toward east), and Z (positive vertically downward). For this analysis, the 1-min X and Y components from the four magnetometers were turned into a spatially varying magnetic field as described in section 3.1.3. Figures 2c to 2j will also be discussed in this section.

South Island GIC Observations
Transpower New Zealand Limited has measured DC currents in multiple transformers across multiple substations in New Zealand's South Island and has archived the data since November 2001. The primary purpose for the DC observations was to monitor stray currents when the high-voltage DC link between the South and North Islands operates in Earth return mode. A detailed description of the New Zealand GIC measurements in the South Island can be found in Mac .
The largest GIC observed during this storm also occurred at 04:46 UT, that is, the time of the sudden commencement. At this time, there are 23 transformer windings in the South Island with operating, reliable GIC measurements. Observations at four of these transformer windings are shown in Figure 3a: T4L at HBW (HWBT4L), T2H at South Dunedin (SDNT2H), T6H at Islington (ISLT6H), and T6H at Manapouri (MANT6H). Figure 3b shows the variation in the peak absolute GICs across the lower South Island on a map and as a bar plot from South to North in Figure 3c at the time of sudden commencement. At this time, the observed GIC for all 23 transformer windings is the largest. These figures demonstrate both the spatial variation as well as differing GIC magnitudes within a given substation. The largest GIC measurement during this storm was at HWBT4L with a maximum absolute GIC of 48.85 A. The transformer with the smallest maximum is Ohau A transformer T7 (OHA T7), where only 0.86 A was seen at sudden commencement.
After the sudden commencement, large, fluctuating GIC was observed throughout the rest of 17 March 2015. These longer durations of lower peak amplitude GICs (i.e., 12:00-14:00 UT) are associated with the long-lasting main phase magnetic activity seen in the magnetic field observations over the same duration.

Modeling Method
We have used a two stage modeling approach to calculate geoelectric fields across the NZ region and GICs flowing through transformers in the South Island, similar to the approach used previously by (Divett et al., , 2018. The main difference is that in the present study, we calculate the spectra across a broad range of frequencies for a real storm. In the first stage, we used the thin-sheet model of VW77 and the thin-sheet conductance model developed in  to calculate the geoelectric field spectra around the New Zealand region. In the second stage, we used Divett et al.'s (2018) transformer-level network model to calculate GICs flowing through every earthed transformer winding in the South Island's high-voltage transmission network.

Thin-Sheet Model to Calculate Geoelectric Field
We apply the VW77 thin-sheet model to calculate the electric field in the New Zealand region. While this region is mostly ocean, we are interested in the calculated geoelectric field on land in the South Island because the field elsewhere in the domain is not necessary to calculate GICs. However, due to the dominance of the geomagnetic coast effect in New Zealand, we need to include the ocean in the model domain.
For our model domain, we represented the 28°× 28°region (roughly 3,360 km × 3,360 km from 158°-186°E to 27°-55°S) around New Zealand's North and South Islands and the surrounding ocean as a horizontal grid of 168 × 168 square cells. In this way, each cell is one sixth of a degree long (roughly 20 km) in both the north and east directions.

Conductance Model of the New Zealand Region
The domain of our thin-sheet map is 1.75× larger than that used by (Divett et al., , 2018, covering a wider area of ocean in the north and east directions compared to those previous studies. The map domain is wider in all directions so that New Zealand's North and South Islands remain in the center of the map. have used the ocean depth and an assumed seawater conductivity of 3 S m −1 to calculate the conductance in the same way as in Divett et al. (2017).

Range of Valid Periods Limited by the Thin-Sheet Assumptions
There are four assumptions in the VW77 thin-sheet model that effectively limit the range of valid periods, is the short-period cutoff and T long C is the long-period cutoff. These are explained in more detail in (Divett et al., , 2018. Only two of the four assumptions limit the range in our grid, due to the other two assumptions (d≪δ, and p ≤ δ/3) being less restrictive for our choice of numerical discretization and thin-sheet thickness. Here, d = 20 km is the thickness of the thin-sheet layer, δ is the skin thickness in the medium lying beneath the thin sheet, and p = 20 km is the length of each cell in the numerical grid.
The short-period cutoff, T short C , is determined by the highest conductance in the thin-sheet layer due to the assumption that d η where η is the skin depth in the thin sheet. Or in terms of period, T short c ≫πμ 0 d 2 σ ts . μ 0 is the permeability of free space, and σ ts is the conductivity in the thin sheet. This assumption leads to a different short-period cutoff depending on which region we are most concerned about. Given that we are mostly interested in the electric fields in the South Island, but we know that the coast effect dominates the electric field for most of the South Island , we use a relaxed interpretation of this assumption (where we interpret "≫" to mean "greater than three times") to the region of the South Island, such that T short C = 5 min. We will discuss the validity of this assumption, in the context of our results, in section 6.1.
The long-period cutoff for the range of valid periods (T long C ) is determined by the resistivity and the length of the numerical domain. This limit is due to the assumption that ( 2) where, N = 168 is the number of cells in the numerical grid. In terms of period, this can be written as T long c ≤ ðN−1Þ 2 p 2 πμ 0 σ u . Where σ u is the average conductivity of the underlying medium as a function of the period, over all, the layers in our layered resistivity structure down to the skin depth. σ u is dependent on the period, due to the distance to the skin depth (which depends on the period) changing with the period. We solved this equation by iterating the vertical dissipation of a horizontal plane wave source at ground level over the required number of depth layers, for a range of periods. This yields T long C ≤ 80.5 min. The short-and long-period cutoffs described here define the range of valid periods, and hence, valid frequencies, we calculate using the VW77 thin-sheet model.
There are two reasons why we used a larger domain compared to Divett et al.'s (2017) earlier thin-sheet model in the New Zealand region. The first reason is that the larger domain allowed us to increase T long C from 26.8 min, for the domain used by Divett et al. (2017), to 80.5 min for the present study as a direct consequence of Equation 2. The second reason is that the larger domain removes the influence of an effect associated with the edge of the modeled domain from impinging onto land, for the range of periods we have used. Near the edges of the domain, the calculated E y is unrealistically decreased (or increased), in an area extending up to 4°or 5°from east and west (or north and south) edges of the domain.

Ground-Level Magnetic Field Spectrum Interpolated Across Domain from Sparse Observations
The locations of measured magnetic field variations within the New Zealand region that we have used are shown in Figure 2e where it is clear that they are spatially sparse. Only the South Island is covered with sufficient observation locations to make firm inferences about the magnetic field variations. While the South Island is the region we are most interested in, we also need to ensure that the magnetic field is sufficiently realistic poleward of MDM and also across the rest of the domain, such that we have confidence in the calculated electric field near the coasts. Including MCQ in the interpolation is useful because it constrains the magnetic field south of MDM to a realistic level. Hence, the South Island is sufficiently covered by the variometers in the lower North Island and South Island which restrict the magnetic fields at the top and bottom of the South Island to observed values.
However, all of these locations are at a similar geomagnetic longitude. The significance of this to calculating the geoelectric field is that we do not have any knowledge of longitudinal differences in the magnetic field. Consequently, we have not been able to use the Spherical Elementary Current Systems (SECS) approach that has commonly been applied previously (e.g., Beggan et al., 2013) to interpolate the magnetic field. Instead, we used a linear interpolation of the magnitude and phase of the Fourier coefficients of the northward and eastward components of the magnetic field for each frequency (b x (f) and b y (f), respectively), in the range of valid frequencies. These interpolations of b x for the single period, T = 20 min, are shown in Figures 2e and 2f for the whole domain and as a function of geomagnetic latitude, respectively.
We are interested in the Fourier coefficients for each frequency because we have applied the thin-sheet modeling approach for each of the valid frequencies. There are 1,440 frequency bins for our 1 day of data at 1-min sample period. Of these, 540 are within the valid range. The power spectrum, |b x (f)| 2 of B x at TEW, EYR, MDM, and MCQ, in Figures 2c and 2d, respectively, shows that there is significantly more power in the magnetic field variations at MCQ compared to observations further north, matching the larger variations in the time series. There is also more power in the longer periods than shorter periods for all four locations, 10.1029/2020SW002494

Space Weather
providing another reason to increase the size of the domain so that the range of valid periods includes more of longer periods that contain the highest power. However, there is still significant power in the periods greater than T long C . There is also significant power in the periods below T short C that will drive much of the short term behavior in the GIC time series.
A snapshot of the magnetic field variation is shown at a rapid change in the magnetic field at 16:29 UTC, during the main phase of the storm in Figures 2g-2j. This snapshot shows the interpolated and recombined B x (t) and B y (t) across the NZ region, representative of the spatial variation of the magnetic field across the NZ region as a result of the linear interpolation in the frequency domain. The recombined B x and B y plotted against geomagnetic latitude demonstrate that these recombined curves match observed values (shown by crosses in Figures 2h and 2j, as expected. These B field variations are consistent with an ionospheric electrojet located between MCQ and the southern most edge of the South Island, as expected for an event such as the St. Patrick's Day Storm.

Transformer-Level Network Model
We used the transformer-level network model for New Zealand's South Island, shown in Figure 1. This was developed by Divett et al. (2018), and we apply it with one modification-in the present study, we have changed the path of the transmission line connecting ISL to Livingstone (an unearthed substation southeast of AVI (see enlarged box in Figure 1) so that the modeled transmission line does not traverse the ocean. This modification changed the GICs by less than 1% compared to GICs calculated with the earlier network model. We used this model for the analysis in section 4 and as the base case of the sensitivity analysis in section 5. We calculated the transformer-level GICs flowing through each winding of every transformer in the South Island's high-voltage transmission network using the approach described in Divett et al. (2018). This approach uses the calculated electric field to determine GICs using a modified version of the LP85 matrix method. We have used the same two modifications to the LP85 matrix method as Divett et al. (2018). We use a "virtual node" to represent each voltage bus within a substation so that each transformer winding is represented as a connection between nodes.
Because we are calculating currents due to a broad spectrum of frequencies, we perform the LP85 matrix calculations on Fourier coefficients in the frequency domain. Accounting for this difference, in the LP85 approach, we calculate the spectrum of the vector of currents flowing to Earth from each node, i s (f), using, i s =(1+YZ) −1 j(f). Y is the network admittance matrix, and Z is the earthing impedance matrix, built using electrical DC resistance values measured by Transpower. In that expression, j(f) is the Fourier coefficient of the current sources along each transmission line found by integrating the Fourier coefficient of the electric field, e ! , along the path of each transmission line with transmission line resistance, R line , using jðf Þ¼ In the second modification to the LP85 approach, we find the transformer-level GICs flowing through each transformer winding. We first calculated the node voltage, v n e ðf Þ¼i s ðf ÞR e , relative to local substation Earth for each nth node. R e is the Earth ground resistance measured by Transpower for each substation or 10 10 for each virtual node. We then determined the GICs flowing through each transformer winding using i trans ðf Þ¼ðV n e ðf Þ−V m e ðf ÞÞ=R m n . R m n is the resistance of the transformer winding connecting nodes n and m, supplied by Transpower for all transformer windings in the South Island. i trans refers to the GIC flowing through a low-or high-side winding of a specific transformer (e.g., i ISLT6H is the spectrum of GIC flowing through the high-side winding of the T6 transformer at ISL. In this sense, transformer winding means the high-or low-voltage winding of a transformer, as described in Divett et al. (2018). The high (or low) side winding refers to the windings on the high-voltage (or low-voltage) side of the transformer core of a two-winding or normal transformer, denoted "H" (or "L"), respectively. For autotransformers, the series winding between 220-and 110-kV buses is denoted "H", while the common winding between the 110-and 0-kV buses is denoted "L," also following Divett et al. (2018).

Filtering Measured and Modeled GICs
In order to compare the results of the modeled GICs with measured GICs, we bandpass filtered both spectra to remove contributions from frequencies outside the model's valid range. Before Fourier transformed the 10.1029/2020SW002494

Space Weather
observed GICs to the spectral domain, we downsampled to 1-min sample period. This downsampling was to ensure constant sample period and for consistency with the modeled GICs which arises from the 1-min sample period of the magnetic field observations from INTERMAGNET observatories.
We calculated the electric fields and GICs for several frequencies outside the range of validity, as well as the valid frequencies. We then applied a Gaussian-tapered bandpass filter to the spectra of both modeled and observed GIC. We applied this as a mask in the frequency domain. The Gaussian-tapered filter reduces nonphysical spikes and ringing at both ends of the time series, associated with an otherwise square-edged (or boxcar) bandpass filter. The long and short cutoff periods for this filter, T long C and T short C , are defined by the range of validity of the thin-sheet model, as described earlier. Outside the bandpass region, this filter exponentially decays toward zero.
After filtering, we applied an inverse Fourier transform to i trans (f) to produce the time series of the modeled transformer-level GICs through each transformer winding, I trans (t). We applied the same transform to the observed GIC spectra, after filtering. These modeled and filtered observed time series are discussed in section 4.

Coherence Analysis
We used a coherence analysis to try to quantify the similarity between the measured and calculated GIC spectra for each substation. Coherence, is a common statistical method in spectral analysis that is often used in climate research (e.g., von Storch & Zwiers, 1999) and weather analysis (e.g., Biltoft & Pardyjak, 2009). Coherence is calculated from the absolute value of the cross spectrum, |A obs,model (f)|, and the autospectral densities of the observed, |i obs (f)|, and modeled, |i model (f)|, GIC. Coherence is essentially a correlation coefficient that is frequency dependent. In calculating the coherence, the data in each spectra must first be binned. In our coherence analysis, we have used frequency bins with 10 frequencies per bin. We only calculate coherence for the range of frequencies that are valid for the thin-sheet model. The coherence is always between 0 and 1, showing the degree of linearity between the amplitudes and phases of two spectra within each frequency bin. A value closer to 1 indicates stronger linearity or closer similarity between the two signals amplitude and phase.
We calculated the 95% significant threshold for these coherence indicators, following Biltoft and Pardyjak (2009). F 0.95 is the upper tail critical value of the F(2,ϵ−2) distribution, and ϵ is the number of independent spectra samples used. Coherence greater than this threshold indicates that the coherence estimate is different from zero with 95% confidence for this frequency bin.

Electric Field
The magnitude and direction of calculated electric field, at sudden commencement is shown in Figure 4a and as a time series for Arthur's Pass, in Figure 4b. Arthur's Pass is a location representative of the strong electric field in the center of the South Island. Both figures have been converted to the time domain from the calculated Fourier coefficients after filtering by the same Gaussian-tapered bandpass filter that we applied to GICs to remove contributions from outside the range of valid periods. Due to the limitations in the range of valid periods, the electric fields presented here are very likely weaker than the real electric field would have been during this storm. However, we do not have any measurements for comparison.
The dominant spatial pattern of the electric field at the snapshot at sudden commencement seen in Figure 4a shows a much stronger electric field on land compared to offshore, as expected, due to the lower conductance on land relative to the ocean. The electric field on land is aligned roughly perpendicular with the main northeast-southwest axis of New Zealand, showing that the geomagnetic coast effect dominates the electric field for the whole island, as Divett et al. (2017) described for an electric field driven by a uniform magnetic field.

Space Weather
The highest magnitudes of the calculated electric field occur near the middle of the South Island, although the region of strongest electric field does not cover the full length of the Southern Alps. The region at the south of the South Island where the conductance model is least certain, lack of observational data on the conductivity structure, has a lower electric field than the rest of the Southern Alps.
The largest peak in the electric field shown in Figure 4b occurs at sudden commencement, as expected.
There are also substantial peaks during the main phase of the storm (12-14 UT) that varies in a similar way to the changes seen in the measured magnetic field and GICs. However, the main phase peaks in the calculated electric field are almost as high as the calculated electric field peak at sudden commencement which is not seen in the measured magnetic field or GICs. This is probably because of the broader spectral content in the rapid pulse at sudden commencement compared to the broader peaks during the main phase.

Comparing Modeled and Measured GIC Spectra
For the St. Patrick's Day Storm of 2015, the calculated power spectrum of the GICs through most transformer windings is of similar magnitude to the power spectrum of the measured GICs, for the majority of valid frequencies. There are 23 transformers that we have reliable measurements of GICs for during this storm. For 18 of those transformers, the modeled power spectra have similar magnitudes to the measured power spectra, and the coherence exceeds the 95% threshold, for most frequencies. Each of those 18 comparisons between modeled and observed spectra shows similar features to SDNT2H and MANT6H, presented in Figures 5a and 5c, respectively.
We show a comparison of the GIC power spectra calculated from our model in Figure 5 for four transformers. We show SDNT2H and MANT6H because they are examples that show good agreement between modeled and measured GICs. We show HWBT4L and ISLT6H because they are the transformers that have shown the largest observed GIC magnitudes in the past (Mac . These two transformers also show some of the poorest agreement between modeled and measured GICs of the 23 we compared, the remaining three "poor" comparisons being transformers at Benmore (BEN in Figure 1).
The three Benmore transformers historically have low GIC. It should be noted that low GIC is not an indication of poor correlation between modeled and observed GIC. Cromwell and Ohau A have equally small observed GIC, yet there is a good agreement for these substations. In contrast, HWBT4L and ISLT6H are two of the transformers that typically have high measured GICs flowing through them during geomagnetic storms. An example of this is during the 2001 storm when HWBT4 tripped out of operation, and the static VAR compensator at ISL tripped due to high GICs (as described by Béland & Small, 2004;Marshall et al., 2012;. The highest GICs measured in NZ's South Island have occurred at HWBT4 and ISLT6 (referred to as ISLM6 by Mac . The modeled GIC spectra at both of these transformers are consistently lower by at least an order of magnitude than the observed GIC spectra, across all frequencies, although the difference is a little larger at higher frequencies at ISLT6H.
There is more power in the spectra of GICs for long periods than for short periods, in both the modeled and observed spectra across all of the transformers we have observations for during this storm. This matches the trend seen in the magnetic field spectra for this storm, as shown in Figures 2c and 2d. The match between observed and modeled GICs at short periods appears to be a little better than for the longer periods, within the range of valid periods for most of the 23 transformer windings.
For each of the selected transformers, we show in Figure 5 the power spectra of GIC, the coherence between modeled and measured GICs. The coherence between the modeled and observed GICs is shown in Figures  5b, 5d, 5f, and 5h for transformers SDNT2H, MANT6H, HWBT4L, and ISLT6H, respectively. The dashed blue line in each of these figures shows the 95% significant threshold for these coherence indicators. The time series of modeled and measured GIC are shown in Figure 6.
With the exception of a few individual frequencies, where coherence also dips below the 95% threshold, the spectra of observed GICs are similar to the modeled spectra at the same transformer. Clearly, the coherence is well above the 95% significance threshold for most frequencies for the GICs flowing through SDNT2H and MANT6H, although the coherence is below the threshold for four or six frequency bins of the 22 bins for SDNT2H and MANT6H, respectively. These bins with coherence below the threshold correspond to the 10.1029/2020SW002494

Space Weather
few individual frequencies where the measurements are significantly different to the modeled GIC. SDNT2H had the second highest peak GIC during sudden commencement (as shown in Figure 3b). This is reassuring evidence that the model can provide a similar GIC spectrum to the measured GIC spectrum for a transformer where high GICs have been measured.
The coherence for HWBT4L and ISLT6H is below the threshold for nine frequency bins, reflecting the larger differences between modeled and observed spectra, across a broader range of frequencies than for the other two transformers. Further, coherence does not give a good indication of whether the magnitude of the two signals is similar across the whole spectrum, because a constant mean offs et alone does not prevent the two spectra from being coherent. This can be seen in the spectra and coherence for HWBT4L and ISLT6H in Figure 5. While the coherence for both of these transformers exceeds the threshold at many of the frequency bins, the power in each modeled spectra is significantly less than the observed power, for the

10.1029/2020SW002494
Space Weather entire valid range of frequencies. Hence, for a complete picture of how well the modeled GIC spectrum represents the observed spectrum, both the coherence values and the two power spectra need to be compared.
We also show the sum of coherence per sample, Σγ 2 /N, across the range of valid frequencies for each transformer in Figures 5b, 5d, 5f, and 5h. This sum can be compared with the same 95% significance threshold, with caution, because we have normalized the sum by the sample length, N. Of the transformers in Figure 5, the sum is highest for SDNT2H (Σγ 2 /N = 0.117) and MANT6H (0.112). This sum is considerably lower for HWBT4L (0.0973) and ISLT6H (0.0859). This sum of coherence is useful for giving a quick, quantitative estimate of the similarity of the modeled and observed spectra of the GIC flowing through a given transformer, across the whole range of valid frequencies, with a single metric. However, it is not by any means a complete comparison of the two spectra and also needs to be used in combination with the power spectra. By looking at the power spectra in Figure 5, the reader should notice by eye a better agreement between the modeled and observed GIC spectra at SDNT2H and MANT6H. This agrees with the larger sum of coherence seen at these transformers.
It is concerning that the modeled spectra is least like the observed spectra for the two transformers with two of the three highest peak measured GICs that occurred during this storm (see Figure 3c). At ISLT6, we suspect that either (1) the Transpower-supplied resistances we used in the model are incorrect for one or more of the neutral Earth resistors, (2) the neutral Earth resistor is faulty or installed incorrectly, or (3) there is a fault with the measurement installation or calibration. This is because the three other transformer windings for which we have observed GICs at ISL show very good agreement between modeled and measured GICs and are electrically similar to ISLT6. Those other transformers (ISLT3H, ISLT7H, and ISLT9H) are connected to the rest of the network in the same way as ISLT6; they are all connected in parallel between the 220-kV bus and the Earth bus; see

Space Weather
ISLT9H have nominal resistances and neutral Earth resistances within 11% of the others. Due to the similarity of these four transformer windings, and their identical connections to the rest of the electrical network, we would expect that differences between measured and modeled GICs would also be similar for each of the four transformer windings. This expectation is based on the understanding that the only difference between the level of GIC flowing through each of these four transformer windings should be due to Ohm's Law, which provides a linear relationship between GIC and the resistance of each transformer winding.
We are working together with Transpower to understand whether any of these three possibilities could be correct. Because each of these four transformers is connected to the same transmission lines in parallel, the induced voltage difference across each transformer winding is identical. This rules out any differences between the voltage across each of the four transformers due to the calculated electric fields, leaving only the application of Ohm's Law across the transformer windings to explain the difference in GIC between each of these transformers.
Understanding the reason why the modeled GIC spectra is lower than the measured GIC at HWBT4 is more speculative. While we only have measurements of GICs for one of the four transformer windings at HWB, we do have observations at the nearby SDN. At SDN, the modeled GICs match the measured GICs very closely (Figure 6a). These substations are both within Dunedin City and are spatially separated by only 4.6 km. Although in reality there are significant local variations in geoelectric structure between the sites-SDN is on low-lying, swampy, coastal sediment that is saturated by highly conductive seawater; HWB is on a steep, rocky hillside-in the thin-sheet model, they are both in the same conductance cell with a value more representative of HWB than of SDN. While some interpolation of electric fields has been carried out the close proximity of HWB and SDN means that the calculated electric fields are essentially the same at the two substations. The two substations are directly connected by 220-kV transmission lines with resistance of 0.26 Ω and via the unearthed TMH substation (also in Dunedin) by transmission lines that total 0.25 Ω. At first sight therefore, due to their close proximity and this strong coupling between the two substations, we expect comparison between measured and modeled GICs to be similar.
There are two differences between the two substations which we speculate might help to account for the difference in the model results and the mismatch between modeled and observed GICs at HWB. First, the difference in surface geology is reflected in the measured resistance between the Earth mat at each substation and local Earth. The values are 0.03 Ω at SDN and 0.23 Ω at HWB. With the same potential difference across the ends of transmission, lines leading into the substations this difference in grounding resistance would argue for GIC at SDN to be almost 10 times larger than at HWB, although the similarity in observed GIC magnitudes at the two substations suggests that this difference is in fact clearly balanced by other factors. Second, in addition to the common transmission lines connected to SDN and HWB from Gore (GOR) and Roxburgh (ROX), HWB is also connected independently to GOR and ROX by additional lines. Of these, the line between HWB and ROX takes a significantly different route than does the common line connecting both SDN and HWB to ROX and twice crosses a relatively significant conductance boundary in the thin-sheet model. The contribution of this line to the overall model GIC at HWB may also therefore partly account for the observed mismatch. The line from HWB to GOR may similarly have an effect although this line does not cross conductance boundaries and follows a relatively similar path to the common line between GOR and the two Dunedin substations.
To fully represent the variations in conductance on a finer scale will require further magnetotelluric (MT) studies in this region. We have recently taken MT measurements at one location near HWB . Utilizing the Fourier coefficients for the St. Patrick's Day Storm we can calculate electric fields from these MT measurements. This is shown in Figure 7, and they highlight the significant influence of local structure and topography on the electric field direction that we do not see in our modeled geoelectric fields. At HWB, we find that the E x thin-sheet electric fields are smaller than the MT measurements while the E y thin-sheet fields are a little larger. Rough calculations suggest that these differences may lead to a difference of 20-30°in the orientation of the total electric field given by the thin-sheet model compared to the MT measurements. Such differences between the thin-sheet calculated electric fields and those from the MT measurements are an example of how the two approaches differ in resolution.

GIC in the Time Domain
The comparison of the modeled GIC time series, reconstructed from the spectrum by inverse Fourier transforming, and the observed GIC time series highlights the strengths and weaknesses with the thin-sheet modeling approach for GICs. These comparisons are shown in Figures 6a-6d for SDNT2H, MANT6H, HWBT4L, and ISLT6H, respectively. We show the time series of modeled GIC after applying the bandpass filter. The modeled time series is compared with the downsampled, bandpass filtered observations, shown in these figures. The full scale of the raw observations is shown in Figure 3.
For many peaks in the GIC time series seen in this storm, the timing and duration of the observed GIC peaks are recreated well in the modeled GIC. We conclude we can have reasonable confidence that we can model the onset time and duration of peaks in the GIC time series, even with the limited range of valid frequencies.
There are a few notable exceptions to this (e.g., around 9:30 UTC for SDNT2H). However, the main limitation of the technique is that the magnitude of the modeled peaks is usually significantly lower than the magnitude for the raw observations (e.g., HWBT4L and ISLT6H). The magnitudes of the modeled peaks are closer to the filtered peaks because both of these time series exclude the same low and high frequencies.
Without the contributions of the high and low frequencies to the time series, the peaks in the observed time series are significantly lower than the peaks in the raw observations.
For instance at sudden commencement, the peak observed GIC at HWBT4L was 48 A in the "raw observations" but the peak modeled GIC of 6.8 A was only one seventh of that. Even the filtered measured GIC was 16 A, twice the modeled value. In contrast, the model overestimates the GIC at sudden commencement for SDNT2H. Here, the 18 A peak in modeled GIC is twice the peak in the filtered measured GIC. But even for this transformer where the modeled and observed spectra look similar, the model cannot hope to predict the raw measured peak of 45 A without contributions from power from outside the range of valid frequencies. Throughout the day of storm, modeled GIC similar or a little lower than the filtered observations for each later peak is seen, but neither is as large as the peaks in the raw observations. The difference between the scale factors at sudden commencement and the rest of the storm period is probably due to the different spectral content between the sharp change at sudden commencement and the broader peaks later.
This underprediction of the peak GIC at sudden commencement when compared to the raw observations, even in the transformer with the best spectral match, is extremely significant to our ability to determine the storm time hazard from GICs to transformers. These differences highlight just how significant the assumptions of the thin-sheet model are to GIC modeling, given that power in the frequencies outside the range of validity clearly contribute a significant proportion to the magnitude of peak GIC, both at sudden commencement and at later times.

Space Weather
Further differences between the observed and modeled GICs are likely due to the inherent approximations in our approach. The uncertainties in conductance and spatial variation of the magnetic field variations are significant, particularly in the south of the South Island where the conductance model is based on local geology. We have already discussed how the 20-km discretization of our grid misses some local variation in geoelectrical structure around HWB and SDN. This choice of discretization may further be significant to the modeled GICs in other parts of our domain.
To summarize, the modeled and observed spectra and time domain for HWB and ISL GIC show the least agreement, which we believe is due to two issues. The first of these are the limits in the valid frequency range for the modeling, while the second are difficulties with modeling the electric field as discussed in section 4.2.
In contrast, the modeled and observations between SDN and MAN are quite good inside the frequency range limits for which the modeling is valid.

Trends in Spatial Variation of GIC
For the St. Patrick's Day 2015 Storm that we have modeled, we see the same large spatial variation in magnitude of GIC that have been seen in previous studies of the South Island of New Zealand  and the United Kingdom (Beggan et al., 2013). For instance the maximum modeled GIC at MAN (3.4 A) is significantly less than at SDN (18.0 A), as shown in Figures 6a and 6b, respectively. This is due to the combined effects of regional variation in the electric field along the path of the transmission lines, the difference in the way transmission lines connect into transformers within a substation, and the different resistances of components of the network.

Sensitivity to Variations in Input
We ran our model with variations to the inputs to ascertain which input affects the calculated GICs the most. The four inputs that we have tried varying are the conductance values for land, the magnetic field, the transmission line paths, and electrical network resistances. We note that in most cases we have made large changes to provide a qualitative indication of the relative importance of each factor, rather than undertaken a detailed quantitative test. The later would be very challenging due to the uncertainties in multiple parameters. The sensitivity analysis presented below should help guide future research by suggesting where we should focus the most effort to improve modeling. We present the change in the spectrum of GIC at SDNT2H for each of the four modeled variations compared with the base model and observations in Figures 8a, 8c, 8e, and 8g, respectively. We selected SDNT2H for this analysis because the differences between the base case and each variation are representative of those for other transformers. We also show the coherence and relative difference, |i(f) base −i(f) var |/i(f) base , between the GIC spectrum for the base mode, i(f) base , and each of the four variations, i(f) var , in Figures 8b, 8d, 8f, and 8bh, respectively.

Uniform Conductance on Land
The first variation we looked at was changing the thin-sheet conductance map to assume that the whole of New Zealand's land has a constant conductance. In this test, we still used the varying depth of the ocean to calculate conductance in the ocean. By assuming that the land has a constant conductance, we are attempting to assess how big the impact of the spatial variation of the conductance of land is to modeled GIC.
We initially used a value of conductance of τ = 1 S, representative of the rocky mountain backbone of New Zealand. This constant value has previously been used by Nakamura et al. (2018) for Japan's main islands to model electric fields and GIC in absence of any more detailed information about the ground conductance. Using this value for New Zealand resulted in very small differences between the base model of GIC and the model with constant conductance. The differences were less than 5% for all frequencies and on the order of 1% for most frequencies.
Because this difference was so small, we also tested a more extreme case where we set the constant conductance to τ = 500 S on land. This represents saturated sediments or swamp. The spectrum for GIC at SDNT2H for this test, shown in Figure 8a, is also similar to the base modeled GIC. The relative differences, shown in Figure 8b are significantly bigger than when the land is assumed to be all rock, around 50-60% for most frequencies. Higher frequencies show a smaller relative difference than lower frequencies.
The coherence of the two GIC spectra, also shown in Figure 8b, is close to 1 for all frequencies, giving a sum of coherence of 0.981 for all transformers. This simply indicates that the variation in GIC with frequency is almost the same for both cases, as would be expected as we have not changed the spectral content of the magnetic field input, only the spectral response of the electric field due to the different thin-sheet map. However, it seems that the significance of the varying thin-sheet conductance on the land is significant, at least when rather extreme variations are considered. In the future, it would be interesting to test a smaller scale variation of the thin-sheet conductance such as moving the boundary between regions of different conductance or changing the value of a single conductance region.

Real Transmission Line Paths
Finally, we built a network model using the true paths of the transmission lines instead of assuming that each transmission line follows a straight path between substations. This has been possible because

Space Weather
Transpower has recently supplied us with the detailed path of each transmission line. Applying these paths, we find the power spectrum is reduced by around 5% for SDNT2H by removing the assumption that the transmission lines follow straight lines between substations, as shown by the spectrum and relative difference with our base case for this variation at SDNT2H in Figures 8c and 8d, respectively. We find the reduction at transformers with relatively straight lines anyway to be smaller.

Uniform Magnetic Field Input
We also ran the model with a spatially uniform magnetic field that is the same as that measured by INTERMAGNET at EYR for the day of the 2015 St. Patrick's Day Storm. This removes the spatial variation of the magnetic field that we determined from the three other magnetic field observation sites around the region, while retaining the changes in time. This variation resulted in a similar difference in the GIC spectrum for each transformer, as shown for SDNT2H in the power spectrum in Figure 8e and more clearly in the relative difference and coherence between the base model and this B EYR model in Figure 8f. The relative difference and coherence vary significantly more with frequency compared to the previous conductance variation test. This is because the power spectrum of the magnetic field is different between the uniform magnetic field and the base case. The relative difference is up to 100% of the GIC for several frequencies.
The coherence ranges from 0.6 to 0.9 as the frequency increases over the range of valid frequencies, as shown in Figure 8f. The sum of coherence over this range is 0.847 for SDNT2H. This sum of coherence is lower for transformers further away from EYR, as we would expect. This reflects the increasing impact of the difference in magnetic field between the two cases on transformers further from EYR, such as MANT6H (sum of coherence = 0.761) compared to those closer to EYR, such as ISLT6H (sum of coherence = 0.967).

Simplified Network
For this variation, we assume that we know the network structure but we do not possess detailed knowledge about the resistance of each element. Hence, in this simplified network model, every transformer, every earthing resistance, and every transmission line is represented by a 0.5-Ω resistor. There are no neutral Earth resistors in this network. A similar assumption was used by Beggan et al. (2013), although they represented each substation as a single resistance instead of representing each transformer. We chose to represent every transformer winding rather than an aggregated substation resistance so that we can compare transformer-level calculations of GIC with our base case.
The power spectrum of GICs flowing through SDNT2H in this simplified network is 40% lower than the base case, averaged across the spectrum, as shown in the spectrum and relative difference to the base case in Figures 8g and 8h, respectively. This is not as large as predicted as we expected network resistances to have a significant impact. This could be due to the fact that the resistance of SDNT2H was 0.23 and the earthed resistance of the substation was 0.03, so changing this to 0.5 Ω only doubles the total resistance. In contrast, the power spectrum is increased by a factor of 10 at ISLT6H, largely due to removing the neutral Earth resistor in the simplified network (2.1-0.5 Ω). Clearly, this is a very significant change.

Summary of the Sensitivity Analysis
Overall, the biggest impact on modeled GICs comes from assuming that the magnetic field has no spatial variation. The second biggest change comes from knowing the true land conductance values compared to assuming a uniform conductance on land, followed closely by knowing the true resistance values of transformers and transmission lines compared to assuming they all have a constant resistance. Changing the transmission line paths to follow their true paths instead of a straight lines have the least impact on the modeled GICs. This shows that to improve accuracy of a GIC model, the most significant information is the magnetic field. This is consistent with the findings of the sensitivity analysis conducted by Beggan (2015) for modeled GICs in the United Kingdom. They found that their modeled GICs were most sensitive to variations in the magnetic field, although they did not know the true network resistances at the time so were not able to compare that variation. Despite this result, it should be emphasized that we picked an extreme conductance of 500 S. If a less extreme value was selected, the relative importance of having access to the true land conductance would likely be reduced.

Discussion
The originally intended purpose of the VW77 thin-sheet model was using measured electric and magnetic fields at a given frequency to test for a suspected conductance anomaly. Thus, we, along with previous studies that use VW77's thin-sheet model for GIC modeling, are using the thin-sheet model outside its originally intended purpose. As such, there has been some confusion in the community about how the thin-sheet model should be applied to help understand GICs in transmission lines. Global geoelectric field models based on a similar integral equation approach to VW77's thin-sheet approach have explicitly been described as being applicable in the frequency domain (Sun & Egbert, 2012). That study further stated that their method was similar to VW77's but with a more modern technique, suggesting that VW77's thin-sheet model should be applied in the frequency domain as we have in the present study.
Other GIC modeling studies have used the thin-sheet model in a different way to that in which we have used it in this paper. Bailey et al. (2018) compared the modeled GIC in the Austrian network to observed GIC at a single transformer in the network using a similar thin-sheet and network modeling approach to the present paper. However, Bailey et al. (2018) applied the thin-sheet method in the time domain by driving it with a series of snapshots of the rate of change of the magnetic field at each minute of the model for a single frequency. They then compared the model output to a time average of the GIC. It is not clear to us how this approach could work in New Zealand. To start with, they are only modeling the response of the ground to a single frequency but assuming that all of the energy in the magnetic field signal is at this single frequency.
As the spectra of the magnetic field in Figure 2 show, there is energy in a wide spectrum of frequencies. If we were to apply the energy as if it occurred at a single frequency, the magnetic field would be a time-varying sine wave with magnitude and phase coming from the Fourier coefficient for that frequency. The time-varying rate of change that Bailey et al. (2018) applied to their model was clearly not a sine wave. So it is difficult to understand how any approach that uses the time-varying signal of B or dB/dt as the input to the thin-sheet model is going to get a result that can be sensibly compared to GIC. This is why we chose to run the simulations in the frequency domain in the present study.

Modeling Assumptions About the Period
As described in section 3.1.1, we chose a short-period cutoff of T short C = 5 min. We chose this period based on our interpretation of VW77's modeling assumptions and our assumption that it is valid to use the calculated geoelectric field across the South Island for periods shorter than a more strict short-period cutoff for the whole domain. The whole domain includes 8,000-m-deep ocean with a conductance 144× greater than the greatest conductance on land in the South Island. Hence, restricting the range of the model's validity based on the deepest ocean in the region, combined with the long-period cutoff of T long C = 80.5 min, would mean that VW77's model is not valid, for any periods at all, in regions of our model domain that have even a small area of deep water.
This strict interpretation of the limits of validity does not seem to be true around New Zealand, because for most transformers, the modeled GICs appear to compare reasonably well with observations. Divett et al. (2017) showed that the same thin-sheet modeling approach we apply resulted in 10-min induction vectors which were close to the induction vectors determined from measurements (Chamalaun & McKnight, 1993). Further, the GICs we calculated in section 4 demonstrate that the electric fields calculated using the thin-sheet approach result in GIC power spectra that are a reasonable match to those of measured GICs, at least for the range of periods between 5 and 80.5 min.
Further, the results of the present paper appear to be valid for a larger range of periods than should strictly apply due to the conductance of the shallow ocean, where the short-period cutoff would be 12 min. We made the choice to make T short c = 5 min even though  concluded that the direction of the coastline dominates the direction of the modeled electric field. This effect can dominate the 3-D modeled geoelectric field on similar islands (i.e., the United Kingdom and Ireland, Ivannikova et al., 2018; and Japan's main islands, Nakamura et al., 2018, respectively) and even the peninsula of Sweden (Rosenqvist & Hall, 2019). Ivannikova et al. (2018) also showed that the coast effect spreads across the whole of Ireland for periods longer than 50 s and that the coast effect is stronger for longer periods. Nakamura et al. (2018) also found that the coast effect is stronger near curved coastlines. In this context, we would expect that the thin-sheet modeling approach would only be valid for modeling GICs at periods greater than 40 min, where the thin-sheet modeling assumptions hold for the shallowest bin of coastal ocean around the South Island as well as on land. However, our results suggest that the thin-sheet modeling approach is probably valid for calculating GICs at most locations in the South Island, for periods shorter than 40 min. By extension, this is possibly also the case in other locations where the coast effect dominates, such as Japan, Sweden, the United Kingdom and Ireland, and even up to 100-200 km from the coast in large continents such as North America.
A more modern geoelectric model such as that developed by Püthe and Kuvshinov (2013) or Nakamura et al. (2018) would avoid the frequency limitations of the VW77 code. Those models appear more appropriate for GIC modeling than the VW77 model we have used in the current study. This is particularly apparent through our demonstration that the frequency limitations in the VW77 approach significantly reduce the predicted peaks in a GIC time series relative to those seen in observations. However, the more modern approaches are significantly more computationally intensive than the model used in the present study.
The transfer function method of Ingham et al. (2017) demonstrates an alternative method to calculate the GIC spectrum for each transformer directly from the magnetic field spectrum. This method requires a time series of GIC for each transformer we want to model, unlike the method in the present paper that can calculate the GICs at other transformers as long as enough detail is known about the whole network. The transfer function method is also unable to reflect any changes to the network that are made for operational reasons or to mitigate the effects of GICs during a storm. The method of Ingham et al. (2017) does, however, have the advantage that it can be applied without detailed knowledge of the electrical transmission network or ground conductance and is considerably less computationally intensive than the approach in the present paper.

Conclusion
We have developed and validated a model for calculating transformer-level GICs in the electrical transmission network of the South Island of New Zealand in two stages. We use VW77's thin-sheet modeling approach to calculate the spectrum of electric fields in the region from the spectrum of the ground-level magnetic field variations during the St. Patrick's Day Storm of 2015. We interpolate this magnetic field variation from observations at four locations that span the South Island. We do this in the frequency domain rather than the time domain because the response of the Earth is frequency dependent. The thin-sheet model is only valid for midrange frequencies, due to modeling assumptions built into the VW77 approach. Our ground conductance model is based on MT surveys, rock-type maps, and bathymetry for the New Zealand region.
We then calculate the spectrum of GICs flowing through each transformer winding in the network by imposing the calculated electric field spectra on our network model of the South Island's transmission network. In our network model, we represent each transformer winding, neutral Earth resistor, Earth to ground resistance, and transmission line in the network by a DC electrical resistance value, supplied by Transpower, assuming transmission lines follow straight lines between substations for the base case.
We have validated this modeling approach, in the frequency domain, for the St. Patrick's Day Storm by comparing calculated GICs with the GICs measured through each of 23 transformers in which Transpower had reliable GIC observations. We found that the thin-sheet model is good at calculating the GIC spectrums for most of the midrange periods, between the high and low cutoff periods. We believe this study is the most extensive attempt to date to attempt validation of GIC modeling during a real-world storm. We found that there is reasonable agreement between modeled and measured GICs for 18 out of those 23 transformers (78%), within the range of valid periods.
In the time domain, this translates to being able to hindcast the timing and duration of the peaks in a GIC time series fairly well. The magnitude of most peaks compares favorably with the bandpass-filtered measured GICs at the 18 transformers with similar spectra. This gives confidence that the model is representing the relevant physics correctly, within the valid range of frequencies. However, the modeled GIC time series significantly underestimates the magnitude of GIC peaks (during both sudden commencement and substorms), compared to the raw measured GICs. This underestimation demonstrates the importance of the power in the spectrum outside the range of valid frequencies to calculating GICs. This is a very significant limitation to the VW77 thin-sheet modeling approach we have used, because understanding the peak magnitude of GICs is very significant to the potential damage to transformers during significant space weather events. The VW77 thin-sheet modeling approach is still useful for comparing the relative impact on transformers in a network or changes to the magnitude of GICs due to mitigation efforts.
We also conducted a small sensitivity analysis to explore which potentially unknown inputs to the modeling have the biggest impact on modeled GICs. We found that the modeled GICs are most sensitive to a change in magnetic field input. Representing the land as a constant, swampy, or rocky conductance made the second biggest change to GICs. Using the true transmission line paths instead of a straight line approximation had the least change to GICs In future work, it is potentially possible to find a scale factor to adjust the modeled time series of GIC to match observations for each transformer by comparing measured and calculated GIC over several storms. This approach would need to take into account the different spectral content in the sudden commencement and main phase and would potentially be quite different for each storm. However, a more modern geoelectric model such as Püthe and Kuvshinov (2013) or Nakamura et al. (2018) would avoid this issue.