Prediction of geomagnetically induced currents (GICs) flowing in Japanese power grid for Carrington-class magnetic storms

Large-amplitude geomagnetically induced currents (GICs) are the natural consequences of the solar–terrestrial connection triggered by solar eruptions. The threat of severe damage of power grids due to the GICs is a major concern, in particular, at high latitudes, but is not well understood as for low-latitude power grids. The purpose of this study is to evaluate the lower limit of the GICs that could flow in the Japanese power grid against a Carrington-class severe magnetic storm. On the basis of the geomagnetic disturbances (GMDs) observed at Colaba, India, during the Carrington event in 1859, we calculated the geoelectric disturbances (GEDs) by a convolution theory, and calculated GICs flowing through transformers at 3 substations in the Japanese extra-high-voltage (500-kV) power grid by a linear combination of the GEDs. The estimated GEDs could reach ~ 2.5 V/km at Kakioka, and the GICs could reach, at least, 89 ± 30 A near the storm maximum. These values are several times larger than those estimated for the 13–14 March 1989 storm (in which power blackout occurred in Canada), and the 29–31 October 2003 storm (in which power blackout occurred in Sweden). The GICs estimated here are the lower limits, and there is a probability of stronger GICs at other substations. The method introduced here will be immediately applicable for benchmark evaluation of low-latitude GICs against the Carrington-class magnetic storms if one assumes electrical parameters, such as resistance of transmission lines, with sufficient accuracy.


Introduction
Geomagnetically induced current (GIC) is one of the important aspects of space weather because it can give rise to a severe, irreversible damage on power grids, which are fundamental infrastructure for our daily lives. It is well known that the power grids at high latitudes are exposed to the risk. Wik et al. (2009) summarized problems of ground-based systems arising from GICs in Sweden, showing that GIC of more than 300 A for a transformer was measured on 6 April 2000 and a power blackout occurred in Malmö on 30 October 2003. A few minutes before the blackout, the calculated GIC reached approximately 330 A (Wik et al. 2009). In the Finnish power system, GIC of 201 A was measured at the Rauma substation on 24 March 1991. This was the largest GIC in Finland as of 2007 ). However, problems caused by GICs had not been reported on the Finnish power grid ). On the same day, GICs of more than 110 A were recorded at the Radisson substation in Canada, which led to the loss of a high-voltage direct current link (Bolduc 2002).
GIC flows in electric circuits between the power transmission lines and the ground, so that GIC is closely associated with the geomagnetically induced electric field (GIE). GIE can be directly observed by dividing the voltage between two electrodes in the ground by the distance between them (Fujii et al. 2015). There is a clear tendency that variations of the GIE are larger at high latitudes and smaller at low latitudes. The border is suggested to be located at geomagnetic latitudes of approximately 50° (Pulkkinen et al. 2012;Ngwira et al., 2013). Although the occurrence frequency is low, relatively large-amplitude GICs are often observed even at the lowlatitudes; 22.5 A at the Islington substation (− 46.6° geomagnetic latitude) in New Zealand on 6 November 2001 (Marshall et al. 2012), 55.8 A at the Ling'ao nuclear power plant (13.1° geomagnetic latitude) in China on 9 November 2004 (Liu et al. 2009), ~ 40 A in Japan on 6 November 2001 (Kappenman 2004) and 129 A in Japan on 31 (METI 2015. (Note that the Japanese facility at which 129 A of GIC was recorded is not open to public.) Such large-amplitude GICs have been recorded during the intervals of large magnetic storms. A question is how large GICs will flow in the low-latitude power grids, and what will happen if severe magnetic storms occur.
The Carrington event that took place on 1-2 September 1859 is one of the most severe magnetic storms recorded by magnetograms (Tsurutani et al. 2003;Cliver and Dietrich 2013;Kumar et al. 2015;Hayakawa et al. 2019). A white-light flare was observed on 1 September 1859 (Carrington, 1859), which might correspond to X45 (± 5)-class flare (Cliver and Dietrich 2013), resulting in an extremely large magnetic storm on 2 September 1859 (Hale 1931). The recent study evaluated the occurrence probability of the Carrington event for the next decade of 0.46-1.88% with a 95% confidence (Morina et al. 2019). If the Carrington event happens again, large-amplitude, rapidly changing geomagnetic disturbances will occur, which will affect the power grids seriously. Pulkkinen et al. (2008) estimated the maximum amplitude of the GIE to be about 4 V/km for the Carrington-class event on the basis of statistics for the occurrence of the GIE at high latitudes. Introducing a technique for scaling the spectral characteristics of the geomagnetic fields, Winter et al. (2017) calculated the maximum amplitude to be 4-20 V/km in UK. The purpose of this paper is to estimate the maximum amplitude of GIEs in Japan and the lower limit of the GICs flowing in the Japanese power grid for the Carrington-class event on the basis of timeseries data of the geomagnetic disturbances (GMDs) recorded at Colaba, India, in 1859.

Data
We used 1-min geomagnetic field (B) and geoelectric field (E) data acquired at Kakioka (36° 13′ 56″ N, 140° 11′ 1″ E). The use of the 1-min data may underestimate results, in particular, for rapid variations in comparison with that of 1-s data. Because we focus primarily on the geomagnetic variations associated with storm main and recovery phases, the underestimation would be negligible. The x-and y-directions refer to the geographical north and east directions, respectively. The geoelectric voltage had been observed at Kakioka since 1932 until 2021. Since the 1990s, electrodes (1 m × 0.5 m) made with copper had been placed with a face-to-face distance of 180 m in the north-south direction and 190 m in the east-west direction in Kakioka (Fujii et al. 2015). The geoelectric field is obtained by dividing the geoelectric voltage by the distance between the electrodes. The detailed information about the measurements and characteristics of the geoelectric field are described by Fujii et al. (2015). We used GIC data observed at 3 substations, Shin-Fukushima (SFS), Shin-Tsukuba (STB), and Shin-Fuji (SFJ). These substations are located in a suburban area of Tokyo (27°-29° geomagnetic latitudes) as shown in Fig. 1, and are connected with the 500-kV transmission lines. A clamp-type current probe was attached to a lead line between neutral of a transformer and the ground at each substation (Watari et al. 2021). The GIC data were sampled at 10 Hz, and was automatically transferred to Kyoto University via a cellular network.  Figure 2 summarizes B and E observed at Kakioka, as well as GICs observed at the SFS, STB and SFJ substations during the magnetic storm that occurred on 20-21 April 2020. The original B and E values were subtracted by those observed at 0220 UT on 20 April 2020, just before the beginning of the magnetic storm. Hereinafter, we refer B and E to geomagnetic disturbance (GMD) and geoelectric disturbance (GED), respectively. Shortlived spikes, which are probably due to non-geophysical reasons, are included in the original E data. We removed the spikes by extracting the data exceeding 10 3 mV/km, and filling the gap by linearly interpolating the values in between. Spikes are also included in the GIC data, which were removed by the moving average method with a window size of 30 s. B x showed a rapid increase at 0230 UT on 20 April 2020. This rapid increase was recognized to be a storm sudden commencement (SSC) by the Kakioka Magnetic Observatory (https:// www. kakio ka-jma. go. jp/ obsda ta/ datav iewer/ en). After the SSC, B x showed a negative excursion with its peak taking place at ~ 1300 UT on 20 April 2020. According to the real-time Dst index provided by World Data Center for Geomagnetism, Kyoto (http:// wdc. kugi. kyoto-u. ac. jp/ dstdir/), a magnetic storm occurred with the minimum realtime Dst of − 59 nT. In response to the SSC, E x and E y decreased to − 6 and − 40 mV/km, respectively. After the SSC, E x and E y showed a positive excursion associated with the development of the ring current. In the course of the positive excursion, there were, at least, two largeamplitude impulsive intensifications of E x and E y , peaking around 0852 and 1248 UT. The maximum amplitudes of E x (E y ) were 13 (75) and 12 (66) mV/km, respectively. The dominance of |E y | (east-west component) over |E x | (north-south component) is a common feature at Kakioka during the magnetic storms (Watari et al. 2021;Fujii et al. 2015). We also note that the minimum peak of B x during this storm did not coincide with the moment at which the GED was zero. This implies that the GED is not simply proportional to the time derivative of the geomagnetic field (∂B/∂t) because of the finite underground conductivity (Cagniard 1953). The amplitudes of the GICs reached the maximum values of 3.4, 0.8 and 1.6 A at SFS, STB and SFJ, respectively, around 1248 UT. The polarity of GIC at SFS is opposite to that at STB, which can be explained by current continuity (Beggan 2015). SFS is located near the eastern edge of the power grid, whereas STB is located west of SFS as shown in Fig. 1. When the eastward GED (positive E y ) is dominant, some GICs flow into the grid from the ground at STB, and flow out of it at SFS. The amplitude of GIC at SFS is larger than at STB. The difference may be explained by the number of simultaneously operating transformers at these substations, and the "edge effect" (Boteler and Pirjola 2017). Transmission lines are further connected to the southwest of STB. This means that STB is located between the edges, and that the current flowing into the grid from the ground is, in part, canceled by the current out of grid. The GIC at SFJ seems not to be correlated with those at SFS nor STB. This tendency is pointed out by Nakamura et al. (2018), who explained it in terms of uneven, localized ground conductivities, together with topology of the network.

Results
We calculated E from B on the basis of the convolution theory (Cagniard 1953; Love and Swidinsky 2014) as (1) where μ 0 and σ are the magnetic constant and the representative ground conductivity, respectively. The matrix C is given by The matrix G is called a time-constant galvanic tensor (Groom and Bahr 1992), which describes local distortion of the geoelectric field due to inhomogeneity of ground conductivity. Based on the observations of magnetic and geoelectric fields at KAK during the Halloween magnetic storm on 29-31 October 2003, Love and Swidinsky (2014) suggested σ and G as We multiplied σ by 1.5 to obtain the best fit with the geoelectric field observed in April 2020. The physical meaning of the ad hoc factor of 1.5 is unknown. The convolution was performed for 10 days. The results are indicated by the red lines in Fig. 2c, d. The relationship between the calculated E and the observed E is shown in Fig. 3. The correlation coefficients between them are 0.72 and 0.93 for the x-and y-components, respectively. The calculated E is not necessarily equal to the observed one since the observed E is known to include contributions from small-scale anomalies in the conductivity of underground, long-term trend, and random noise component (Fujii et al. 2015). The difference between the calculated and observed E is beyond the scope of this paper.
We assumed that GIC is represented by a linear combination of E x and E y as where a, b, and c are coefficients Viljanen et al. 2004). a refers to an offset, probably including noise ). b and c are parameters (3) σ = 5.13 × 10 −4 S/m, depending on the topology of the network and resistances of it. We fitted E x , E y , and GIC observed on 20-21 April 2020 to Eq. (4), and obtained the coefficients a, b, and c for the substations SFS, STB and SFJ. They are summarized in Table 1.
Using the coefficients, we calculated GICs from the observed E, which are indicated by the red lines in Fig. 2e-g. The relationships between the calculated and observed values are shown in Fig. 4. The correlation coefficients between them are 0.91, 0.91 and 0.81 for SFS, STB and SFJ, respectively. It can be said that the GICs are fairly well reproduced by the method using the convolution theory and the linear combination.
Finally, we applied our method to the magnetic records of the Carrington magnetic storm acquired at Colaba (18° 54′ 36″ N, 72° 48′ 36″ E). Figure 5a shows the horizontal component of the GMD (ΔH) observed at Colaba on 2 September 1859. The temporal resolution varied from 5 to 60 min, depending on time. For example, ΔH was recorded at 15-min interval from 0600 to 0700 UT, and 5 min from 0700 to 0910 UT. A sharp, negative excursion was recorded at 0600-0800 UT on 2 September 1859. ΔH reached its minimum value of approximately − 1600 nT at 0700 UT. The features of the ΔH variations were described by Tsurutani et al. (2003) and Siscoe et al. (2006) in detail. Siscoe et al. (2006) discussed two possibilities for the cause of the sharp, negative excursion of ΔH. One is the magnetospheric current, and the other is the ionospheric current. They suggested that the magnetospheric current is likely to cause it. The rapid decrease in ΔH (storm main phase) has been attributed to the development of the ring current (Tsurutani et al. 2003;Keika et al. 2015;Li et al. 2006), which is one of the major magnetospheric currents during the storms. As for the rapid increase in ΔH (storm recovery phase), Li et al. (2006) explained it in terms of the enhancement of the magnetopause current (flowing eastward), whereas Keika et al. (2015) explained it in terms of the recovery of the ring current resulting from a rapid replacement of the dense plasma by the tenuous one in the inner magnetosphere (Ebihara and Ejiri 2003). The contribution from the tail current is probably small in comparison that from the ring current (Ohtani et al. 2001). The field-aligned currents can   et al. Earth, Planets and Space (2021) 73:163 result in GMDs at mid-latitudes, in particular, in the geomagnetic east-west component (D-component) (Ohtani et al. 2018;Nakano and Iyemori 2005). The contribution from the field-aligned current is taken into account in the calculation as mentioned below. If the sharp, negative excursion recorded at Colaba (at ~ 10° magnetic latitude) at 0600-0800 UT was caused by the ring current and/or the magnetopause current, similar GMDs would occur in Japan (at ~ 15-38° magnetic latitude). As was discussed by Tsurutani et al. (2005), the contribution from the ionospheric current is probably small because Colaba was far away from the regions where the equatorial electrojet and the auroral electrojet flows. In addition, the variations of the GMDs at Colaba are fully different from those at high latitudes (Hayakawa et al. 2019). This means that the GMDs at Colaba are hardly explained by the equatorward expansion of the ionospheric current. Because of the reasons mentioned above, we assumed that the magnetospheric current (the ring current and/or the magnetopause current) could have contributed to the extremely rapid negative excursion of ΔH with the large amplitude.
For a benchmark test, we made the following two assumptions. Firstly, the variations of ΔH observed at Colaba on 2 September 1859 were primarily caused by the magnetospheric currents (the ring current and/or the magnetopause current). The magnetospheric currents flow relatively far from the ground, so that it is reasonably considered that the magnetospheric currents could contribute to ΔH in Japan as well. We corrected the latitudinal difference as mentioned below. The longitudinal difference between Colaba and Japan is attributed to the difference in the magnetic local time (MLT). Japan could undergo similar GMD variations if the Carrington event began ~ 4.5 h earlier. Secondly, the amplitude of the D-component of the GMD (ΔD) is assumed to be one half of that of ΔH. For the 13-14 March 1989 magnetic storm (in which the minimum Dst value was − 589 nT), the maximum amplitude of the declination was less than 0.5° at Kakioka, that is, |ΔD| was within 260 nT. Ohtani et al. (2018) summarized the 10 largest hourly ΔD values measured at Memambetsu for 1987 to 2016. According to their results, this assumption is valid for the magnetic storms with the Dst values being less than − 250 nT. The large-amplitude ΔD is probably caused by the field-aligned current (Ohtani et al. 2018;Nakano and Iyemori 2005).
Based on the above considerations, we estimated the x-and y-components of the GMDs at Kakioka as follows: where δ is the declination angle at Kakioka, ΔH is the GMDs recorded at Colaba, and a is a coefficient to correct for the magnetic latitude effect. δ is set to be − 7° (positive eastward) for 2021. ΔD is assumed to be 0 and ± aΔH/2, according to the above consideration. Following Sugiura (1991), we corrected for the magnetic latitude effect, in which the ring current is thought to be the largest at the magnetic equator on the ground. a is given by where λ c and λ k are the magnetic latitudes at Colaba and Kakioka, respectively. Assuming the dipole magnetic field with the geomagnetic pole locations determined by International Geomagnetic Reference Field (IGRF) 1900 and IGRF 2020 (Alken et al. 2021), we obtained λ c and λ k to be 10° and 28°, respectively. That is, a is 0.9.
With the convolution theory, we estimated the GEDs at Kakioka, which are shown in Fig. 5b, c. In Fig. 5b-f, the solid lines indicate the results without applying ΔD, whereas the vertical bars indicate the results obtained when we applied ΔD (= ± aΔH/2) in Eq. (5). The maximum amplitudes of E x and E y are 338 ± 27 and 1990 ± 454 mV/km, respectively. The magnitude of the GEDs (= (E x 2 + E y 2 ) 1/2 ) could reach 2018 ± 445 mV/ km. These amplitudes are about 50-56 times higher than observed in the 20-21 April 2020 storm. Applying the obtained E to Eq. (4), we calculated GICs, which are shown in Fig. 5d-f. The maximum GICs are 89 ± 30, 20 ± 7 and 19 ± 27 A at SFS, STB and SFJ, respectively. The polarity of the storm-time ΔD is known to depend on MLTs Iyemori 2003, 2005). In the worst case, it can be said that the GED could reach ~ 2.5 V/km and that the GIC could reach ~ 120 A. We should note that the calculated GEDs and GICs are probably underestimated since the sampling interval of the GMD at Colaba was 5, or 15 min near the storm maximum.
We repeated the same analysis for the two storms on 13-14 March 1989 and29-31 October 2003. For these storms, the geomagnetic and geoelectric field data at Kakioka are available. The first storm is the largest in terms of the minimum Dst value as of 2021, and caused a power blackout in Québec, Canada (Boteler 2019; Allen et al. 1989). The results are shown in Fig. 6. The (5) B x = a�H cos δ − �D sin δ, B y = a�H sin δ + �D cos δ, maximum value of the observed |E y | was 380 mV/km at 1126 UT on 13 March 1989. At this moment, the maximum values of GICs are estimated to be 14.4, 3.3 and 16.7 A at SFS, STB and SFJ, respectively. Figure 7 shows the results for the second storm, which is known to cause a blackout in southern Sweden (Pulkkinen et al. 2005). The maximum value of the observed |E y | was 521 mV/km, taking place at 0534 UT on 31 October 2003. This moment exactly corresponded to that at which 129 A of GIC was observed at another electric facility in Japanese power grid (METI 2015). The maximum value of |E y | on 31 October 2003 was approximately 4.2 times smaller than that estimated for the Carrington magnetic storm. The maximum values of the calculated GICs were 24.4, 5.2 and 13.8 A at SFS, STB and SFJ, respectively, at which |E y | reached the maximum. The large-amplitudes of |E y | and GICs occurred when B x rapidly increased. A plausible explanation for this rapid increase is an arrival of interplanetary shock, resulting in the intensification of the magnetopause current. However, the solar wind data are unavailable for this interval, so that we cannot confirm it. We note that the maximum value of |E y | took place during the storm recovery phase. During this phase, B x gradually increases (probably associated with the decay of the ring current), resulting in a gradual decrease in E y . When the interplanetary shock arrived at the Earth during the storm recovery, E y further decreased probably due to the combination of the decay of the ring current and the development of the magnetopause current.

Discussion
If the Carrington magnetic storm occurs again, GICs will reach, at least, 89 ± 30 A in the Japanese power grid. A question arises if the power grids face the risk against the Carrington-class magnetic storm. The use of digital protective relays removes harmonics and DC current, and overcomes vulnerability against such failures (METI 2015; Power System Relaying and Control Committee K 2019). Pulkkinen et al. (2005) investigated the cause of the power blackout taking place in Sweden on 30 October 2003, and concluded that the failures related to GIC were primarily caused by unwanted operation of the protective relays. Tokyo Electric Power Company and Japanese transformer makers conducted an experiment on DC excitation of transfers due to GICs, and investigated the relationship between GICs and local heating in a transformer (Takasu et al. 1994). They imposed DC current with amplitude corresponding to GICs of approximately 200 A per 3 phases on an experimental transformer. The temperature increased to 107 °C and 110 °C in the core plate and the core support of magnetic steel, respectively. They concluded that the effects of the local heating on transfers are small in terms of lifetime reduction since the occurrence frequency of such GIC level is low, and the duration of it could be short in Japan. Indeed, the temperature obtained by Takasu et al. (1994) is lower than the maximum temperature limits suggested in IEEE C57.91-2011 (https:// stand ards. ieee. org/ stand ard/ C57_ 91-2011. html) for normal life expectancy loading condition. According to the documents published by the North American Electric Reliability Corporation (NERC 2016), if the maximum effective GIC in the transformer is higher than 75 A per phase (225 A per 3 phases), transformer thermal impact assessments are required to be performed on power transformers (NERC 2016). From these considerations, the predicted GICs shown in Fig. 5 are most likely below the threshold to be concerned. Of course, a careful investigation is necessary to confirm a possible influence on the power grids against the Carrington-class event.
We cannot conclude that the Japanese power grid is safe against the Carrington-class event because we focused on only 3 substations. More than 100 substations and power plants are in operation for the 500-kV class power grid in Japan. There is a high probability that GICs with amplitude exceeding 120 A flow in other substations and power plants. In fact, 129 A of GIC was observed for 1-2 min at another electric facility in the Japanese power grid on 31 October 2003 (METI 2015). Installing the clamp-type current probe to all the substations and power plants is a difficult task and unrealistic. One of the most reasonable ways is to solve the equation of the current continuity for given geoelectric field distribution (Lehtinen and Pirjola 1985). The geoelectric field distribution over Japan has been obtained by 3-dimensional modeling and simulation (Püthe et al. 2014;Fujita et al. 2018;Nakamura et al. 2018). Using the geoelectric field distribution, Nakamura et al. (2018) calculated GICs in the modeled 500-kV power grid that consists of 134 nodes and 160 links. They used arbitrary chosen values for electrical parameters including the transmission line resistance and ground resistance since such electrical parameters are not open to public. To evaluate the GICs at all the substations and power plants against the Carrington-class event with sufficient accuracy, we need to know the correct electrical parameters and geoelectric field distribution. This is beyond the scope of this paper, and will be addressed in future.
A question remains whether the Carrington-class event is regarded as a worst-case event. Severe magnetic storms with expected minimum Dst values less than   et al. 2006). Observations also show that the speed of CMEs increases with the peak flux of soft X-ray, which is a proxy of the solar flare magnitude (Takahashi et al. 2016). When the southward component of interplanetary magnetic field is embedded in the fast magnetic cloud and/or in the sheath region ahead of the cloud, a large magnetic storm is expected to occur (Tsurutani et al. 2003). Using data from the Kepler space telescope, Maehara et al. (2015) found 187 extremely large flares named 'super flares' on 23 solar-type stars. For example, the averaged occurrence rate of the X-100 class flares is about once in 500 to 600 years. ICMEs associated with a super flare may cause an extraordinarily large magnetic storm. Another scenario will be necessary to estimate GICs for super magnetic storms associated with X-100 class flares. It is reasonable to consider that the magnetic storms being much severe than the Carrington event could have occurred in the past, and will occur in the future.

Conclusions
We reached the following conclusions: 1. The method using the convolution theory and the linear combination of the geoelectric disturbances (GEDs) was confirmed to reasonably reproduce geomagnetically induced currents (GICs) at 3 substations in the Japanese ultra-high-voltage (500-kV) power grid. 2. We applied the time-series data of the magnetic field disturbances (GMDs) acquired in Colaba, India, on 2 September 1859 to the method presented in this paper, that is, a convolution theory to estimate GEDs and a linear combination to estimate GICs. GEDs could reach ~ 2.5 V/km at Kakioka, Japan. GICs in the Japanese power grid could reach 89 ± 30 A, at least, at one transformer of the Japanese power grid. The GICs were probably below the threshold level to be concerned, although a technical investigation is necessary to confirm a possible influence on the power grids. The estimated GICs estimated here are the lower limit. We need to evaluate the GICs at all the substations and power plants in the extra-highvoltage power grids against the Carrington-class event. 3. We applied the same analysis to the 13-14 March 1989 and the 29-31 October 2003 storms. The former storm is the largest since 1957 in terms of the minimum value of the Dst index. The maximum values of the observed |E y | are 0.38 and 0.52 V/km, respectively, which are 3.8 and 5.2 times smaller than that estimated for the Carrington magnetic storm.
The maximum values of calculated GICs at SFS could be 14.4 and 24.4 A, respectively, which are 4-6 times smaller than the GICs estimated for the Carrington magnetic storm.