Reproducibility of the Geomagnetically Induced Currents at Middle Latitudes During Space Weather Disturbances

Watari et al. (Space Weather, 2009, 7) found that the geomagnetically induced current (GIC) in Hokkaido, Japan (35.7° geomagnetic latitude (GML)), is well correlated with the y-component magnetic field ( B y ) (correlation coefficients >0.8) and poorly correlated with B x , z and d B x , y , z / d t . The linear correlation with B y would help predict the GIC, if we have capabilities of reproducing the magnetosphere–ionosphere currents during space weather disturbances. To validate the linear correlation with B y for any periods (T) of disturbances, we made correlation analyses for the geomagnetic sudden commencements and pulsations (T = 1–10 min), quasi-periodic DP2 fluctuations (30 min), substorm positive bays (60 min), geomagnetic storms (1–20 h), and quiet-time diurnal variations (8 h). The linear correlation is found to be valid for short periods (cc > 0.8 for T < 1 h) but not for long periods (cc < 0.3 for T > 6 h). To reproduce the GIC with any periods, we constructed one-layer model with uniform conductor and calculated the electric field (IEF) induced by B y using the convolution of d B y / d t and the step response of the conductor. The IEF is found to be correlated with the GIC for long periods (cc > 0.9), while the GIC- B y correlation remains better for short periods. To improve the model, we constructed a two-layer model with highly conductive upper and less conductive lower layers. The IEF is shown to reproduce the GIC with cc > 0.9 for periods ranging from 1 min to 24 h. The model is applied to the GIC measured at lower latitudes in Japan (25.3° GML) with strong B y dependence. The mechanism of the strong B y dependence of the GIC remains an issue, but a possible mechanism for the daytime GIC is due to the zeroth-order transverse magnetic (TM0) mode in the Earth-ionosphere waveguide, by which the ionospheric currents are transmitted from the polar to equatorial ionosphere.


KEY POINTS INTRODUCTION
Geomagnetic disturbances have been known to induce electric fields on the surface of the Earth, which create a potential difference between transformers in the power transmission line system. The potential difference drives electric currents [geomagnetically induced currents (GICs)] in the power lines through the Earthing lines of the transformers (Pirjola, 1983). The GIC is a quasi-steady current, compared to the frequency (50 Hz or 60 Hz) of the power system, and has the larger magnitude at higher latitudes, particularly at the auroral latitudes where the auroral electrojets cause large-amplitude disturbances in the northward magnetic field (B x ) on the ground. The magnetic disturbances often go over 2000 nT and occasionally cause blackouts of the power system, as actually occurred in Canada and USA in March 1989 (Bolduc, 2002).
The GIC is derived from formula relating the GIC and the surface electric field; GIC aE x (t) + bE y (t), where E x,y is the surface electric field measured or calculated from the surface magnetic fields, and a and b are the parameters, which depend on the topology and the electrical characteristics of the system (Pulkkinen et al., 2007;Viljanen et al., 2012;Wei et al., 2013). The surface impedance has been widely used to estimate the electric fields from the magnetic fields through the relationship; E x,y Z·B y , x /μ, where B y , x , μ, and Z are the horizontal magnetic field, magnetic permeability, and surface impedance, respectively. The surface impedance is derived from the ground conductivity and layer thickness through the complex-image method (Boteler and Pirjola, 1998) and from the measured GIC and surface magnetic fields (Pulkkinen et al., 2007).
On the other hand, it was reported that the GIC is much more closely related to time derivatives (dB/dt) than B (deflection from the pre-event value) (Viljanen, 1997;Trichtchenko and Boteler, 2006) and the GIC has been evaluated by quite a few papers from dB/dt using the Faraday's law (Viljanen, 1997;Pirjola, 2000;Carter et al., 2016;Kozyreva et al., 2018). It should be noted, on the other hand, that dB/dt is related to the spatial derivative of the electric field, E, in the Faraday's law. No clear correspondence between the GIC and dB x,y /dt reported by Pirjola (1983) may be an example of an inappropriate application of dB/dt to the GIC. The dB/dt method has been used to assess the GIC even under the equatorial electrojet where the GIC has never been measured (Ngwira et al., 2013;Carter et al., 2016). It should be noted that the induction theory tells us that dB/dt should be convolved with the function of 1/ t √ (t: time) (Cagniard, 1953;Viljanen and Pirjola, 1989) that is a response of the conductor to step function-like magnetic field changes (Cheng, 1959). Thus, the sole usage of dB/dt does not meet the induction theory except that the model is composed of two layers with less conductive upper layer over the highly conductive lower layer (Pirjola, 2010). Watari et al. (2009) demonstrated that the middle latitude GIC in Hokkaido (35.7°GML), Japan, is not correlated with dB x,y,z /dt nor with the x-component magnetic field, B x (deflection from the pre-event value), but very well correlated with the y-component magnetic field, B y , with the correlation coefficients (cc) > 0.8. This result meets the two-layer model with the highly conductive upper layer (Pirjola, 2010). Concerning the strong B y dependence, Watari et al. (2009) suggested that the GIC is a return current of the ionospheric currents carried by the TM 0 mode waves in the Earth-ionosphere waveguide, which was applied to explain the instantaneous transmission of the polar electric field and currents to the equator (Kikuchi et al., 1978). Brändlein et al. (2012) also suggested that the GIC is closely associated with the ionospheric currents by showing diurnal and seasonal variations of the GIC observed at low latitude in northern Chile.
The Hokkaido GIC has been reproduced from the ground magnetic fields using the surface impedance (Pulkkinen et al., 2010). Furthermore, Love and Swidinski (2015) reproduced the geoelectric field (GEF) measured at Kakioka, Japan (27.8°GML), from the magnetic fields at Kakioka using the convolution of dB/ dt and the response of the semi-infinite one-dimensional flat Earth. Love and Swidinski (2014) solved the diffusion equation using the Laplace transformation and applied the function of t √ named linear ramp function for the convolution with dB/dt. The reproduced IEF was plotted in good shape with the observed GEF, but the correlation is not evaluated quantitatively. As overviewed above, there are various methods to reproduce the GIC/GEF from the surface magnetic field such as from B y , dB/ dt, surface impedance and from the convolution of dB/dt and response functions. The variety of the method may be due to many factors affecting the GIC such as directions of power lines and coastlines, 3-D structures of Earth's conductivities (Goto, 2015;Nakamura et al., 2018;Ivannikova et al., 2018), and the propagation mode that transports magnetic disturbances from the ionosphere and magnetosphere into the Earth. In this study, we revisit the GIC in Hokkaido to construct a model that reproduces the GIC from the observed magnetic field, B y . The model is not to clarify the structure of the Earth's conductivity that has been made by other methods like the magneto-telluric (MT) method but is rather a tool designed so as to reproduce the GIC from the surface magnetic field as accurately as possible. As a next step to the accomplishment of the good correlations between the GIC and B y (Watari et al., 2009), we examine if the GIC-B y correlations are valid for any space weather disturbances with different period/time scales ranging from 1 min to 24 h. As shown in the following sections, we found that the GIC-B y correlation depends on the period of disturbances such that cc > 0.8 for short periods (<1 h) and cc < 0.3 for long periods (>6 h). To construct a model capable of reproducing the long period GIC, we calculated the IEF, E x,y using the convolution of dB y,x /dt and the step response of the semi-infinite one-layer conductivity model. The GIC-E x correlation is shown to be much better (cc (GIC-E x ) > 0.9) than cc (GIC-B y ) for long periods, while cc (GIC-B y ) is still better than cc (GIC-E x ) for short periods. To construct a model covering both short and long periods, we built the two-layer model composed of highly conductive upper layer over less conductive semi-infinite lower layer. The two-layer model is shown to reproduce the GIC with cc > 0.9 for periods ranging from 1 min to 24 h. In Derivation of the IEF from the Observed Magnetic field, we formulate equations that derive the IEF from the observed magnetic field in one-and two-layer models. Then, we calculate correlation coefficients of the GIC with B x,y and E y,x induced in the one-and two-layer models in Correlations among Observed GIC, B x,y and E y,x . In order to evaluate the capability of the model for various types of space weather events, we analyzed impulsive geomagnetic sudden commencement (SC), short-period (1 min) geomagnetic Pi2 pulsations, longer-period (30 min-8 hours) DP2 fluctuations, and solar quiet diurnal variations (Sq), isolated substorm magnetic bays, and long-lasting storm disturbances (1-24 h). To examine the generality of the model, we applied the model to the GIC measured at the Shin-Yamaguchi (SYG) substation of the Chugoku Electric Company located at lower latitudes in Japan (25.3°GML). We found that the model well reproduced the GIC at SYG with high correlation coefficients (cc 0.87-0.95) for DP2 and SC events with strong B y dependence. In Discussion, we discuss that the daytime GIC can be connected with the ionospheric currents by the TM 0 mode waves in the Earth-ionosphere waveguide, which carry the ionospheric currents from the polar ionosphere to the equator (Kikuchi et al., 1978;Kikuchi, 2014). We further stress that large-amplitude GICs tend to occur around the midnight during substorms, which raises an issue on the propagation mode from the magnetospheric currents to the ground on the nightside.

Convolution Theorem
The magnetic field, B y , propagates downward in the Earth (a conducting medium) as described by the diffusion equation derived from the Faraday's law, Ampere's law, and Ohm's law. To solve the equations, we use the Laplace transformation that transforms differential equations into algebraic equations and convolution into multiplication. In the transformed equations, the time derivative is multiplication of s (s is the complex number used in the Laplace transformation), and integration is multiplication of 1/s. Following the theory of response of the linear system (Cheng, 1959), the induced electric field (IEF), E x (t), is a response of the linear system (conductor) to the external excitation (applied B y (t)). Letting the Laplace transforms of E x (t) and B y (t) be e x (s) and b y (s), respectively, we express e x (s) as a product of the excitation transform, b y (s) and the transfer function,f(s) (Laplace transform of the impulse response function, F(t)) as shown below. e x (s) f(s) · b y (s). (1) Here, we note that the impulse response is a response of the system to the external excitation in a form of the delta function, δ(t), of which Laplace transform is 1. The inverse Laplace transformation of Eq. 1 gives E x (t) in a form of the convolution integral of the impulse response, F(t) and the external excitation, B y (t) as given by where p refers to the convolution of two functions. The convolution (2) implies that E x (t) is a sum of impulse responses of the system excited by B y recorded from τ 0 to t. We may write the Eq. 1 in a different form including 1/s (time integral) and s (time derivative) as Using Eq. 3, we can write the convolution (2) in a different form as where G(t) denotes the step response of the conductor, a response to the excitation function in a form of the unit step function, U(t)( 1 for t > 0 and 0 otherwise) and B y (0) is the initial value of B y . In the following, B y (0) is assumed to be zero since t 0 is set to the quiet time before the arrival of the disturbances. The convolution (4) is identical to the Eq. 12 of Cagniard (1953). It should be stressed that the IEF is obtained from dB/dt convolved with the step response of the conductor.
In the following, we use Eq. 4 to derive the IEF, while Eq. 2 works the same way. Since the GIC and magnetometer data at Memambetsu (MMB), Hokkaido are sampled every one second, the time t is discrete, and the time derivative and integral in Eq. 4 are replaced with a difference and summation, respectively, as follows.
where we start the summation from i 1 since B y (t ≤ 0) 0 is assumed.

IEF in One-Layer Model
The diffusion equation in the conductor is derived from the Faraday's law, Ampere's law, and Ohm's law as listed below.
where σ and J are electric conductivity and current in the conductor, respectively. The Eq. 6 leads to the following diffusion equations for B y and E x propagating toward the z-direction in the coordinates; x, y, and z directed toward the north, east, and down, respectively.
The Laplace transforms of Eq. 7 are given by where the initial value of B y is assumed to be zero as mentioned above. Transformed solutions are obtained in the following form: We give the unit step function for B y (t) at z 0 to obtain the step response function, which is used to derive the IEF from the convolution with dB y /dt. The coefficients a 1,2 are determined from the following boundary conditions: B y (t) U(t) at z 0 and B y (t) 0 at z ∞. We note that B y in the convolution integral (4) is the magnetic field observed at the surface of the Earth, which therefore includes effects of induced currents as well as external currents flowing in the ionosphere and magnetosphere. Whatever the source currents are, the observed B y is the boundary value for the diffusion equation in the ground.
We thus have the transformed solutions as Substituting z 0 for Eq. 10 and applying the inverse transformation, L −1 {1/ s √ } 1/ πt √ , we obtain the step response function, G(t) for E x as Substituting (11) for (5), we obtain the IEF as E y induced by B x is calculated by replacing B y with -B x in the Eq. 12. We calculated E x and E y with the conductivity, σ 10 −4 mho/ m, and as will be shown, the GIC is well correlated with E x (cc > 0.9) for long-period disturbances, whereas the correlation is still better with B y than with E x for short periods. This result would raise a problem that requires us to use two models to reproduce the GIC, depending on the period of disturbances. To address this problem, we construct a twolayer model as shown in the next subsection.

IEF in the Two-Layer Model
Using the earth-currents and magnetometer data, Owada (1972) showed that the subterranean electric conductivity at Memambetsu has a structure of three layers with depths of 8-20 km, 20-90 km, and 90-170 km and with conductivities higher in the top layer than in the lower layers. The MT method has revealed inhomogeneous distribution of the Earth's conductivity in Hokkaido not only in the vertical but also in the horizontal directions (Satoh et al., 2000;Uyeshima et al., 2001;Uyeshima, 2007). However, since our purpose is to construct a model that is capable of reproducing the observed GIC, we pay our attention to the vertical profile of Owada (1972)'s results and construct a two-layer model with thickness, d 20 km and σ 1 10 −4 mho/m in the upper layer (layer 1) over the semi-infinite less conductive (σ 2 10 −8 mho/ m) layer (layer 2). The parameter dependence of the model will be discussed in the discussion section. We assume the magnetic field be a fixed value at z 0 in the same way as in the one-layer model, and the magnetic permeability is common in both layers. The Laplace-transformed solutions in the layer 1 and layer 2 are given as follows: We give the unit step function for B y at z 0 and employ the boundary conditions at z d as b y and e x being continuous across the boundary. Then, we have the following relations among the coefficients: We obtain Under the condition, k 12 e −2d μσ1s √ < 1 (t ≥ 1), we can use the following series expansion that represents reflections at the boundary between the two layers.
where j refers to the number of reflections and n is chosen so that the summation approaches a steady value (n 50 in the calculation below). Then, we have the coefficients as follows: Substituting Eq. 18 for Eq. 13, we obtain the transformed solutions as follows: The Laplace transform of the step response function, g(s), is obtained by substituting z 0 for e x1 as Using the inverse Laplace transform, we obtain the step response of the two-layer model as Using Eq. 21 in Eq. 5, we obtain the IEF from the convolution,

Pi2 and SC (1-10 min)
To confirm the high correlation between the GIC and B y (Watari et al., 2009), we picked out a Pi2 event with the period of 1 min recorded at Memambetsu (MMB) in the morning sector (0600 MLT) and a geomagnetic sudden commencement (SC) with the preliminary impulse (PI, 1 min) followed by the main impulse (MI, 5-10 min) in the morning sector (0630 MLT). Figure 1 shows B x and B y at MMB (top left) and the GIC (bottom left) observed during the Pi2 event. The GIC is well correlated with B y (cc 0.90), while almost nothing with B x (cc -0.21). The Pi2 often occurs on the nightside during substorms, while also observed on the dayside as the event in Figure 1 is the case (e.g., Han et al., 2004). The daytime Pi2 has been attributed to ionospheric currents flowing from the polar ionosphere to the equator carried by the TM 0 mode waves in the Earth-ionosphere waveguide (Sutcliffe and Lühr, 2010;Imajo et al., 2015). The TM 0 mode waves propagating southward (-x direction) have the magnetic field B y perpendicular to the propagation plane and the electric fields,E x,z in the propagation plane, which transport the ionospheric and ground surface currents with north-south direction (Kikuchi and Araki, 1979). The good correlation between the GIC and B y may indicate that the GIC is the ground surface current transported by the TM 0 mode waves as suggested by Watari et al. (2009) andBrändlein et al. (2012). Figure 1 also shows the IEF in one-layer model, E xI and E yI induced by B y and B x , respectively (top right), and the IEF in the two-layer model, E xII and E yII (bottom right). The correlation coefficient of the GIC with E xI is cc(GIC − E xI ) 0.77, less than cc(GIC − B y ) 0.90, but the correlation with E xII is Frontiers in Astronomy and Space Sciences | www.frontiersin.org October 2021 | Volume 8 | Article 759431 5 cc(GIC − E xII ) 0.96, indicating that the GIC can be reproduced almost perfectly by the two-layer model. E xII is plotted in the frame of the GIC with the dotted curve, whereE xII is scaled to the GIC so that one can see the correlation with the GIC visually. On the other hand, the correlations with E yI and E yII are almost nothing (cc -0.18 and -0.07) in the same way as the correlation with B x , indicating that the GIC has no relations with B x . Figure 2 shows the SC event with the positive PI followed by the negative MI in B x and negative PI followed by the positive MI in B y . The GIC is well correlated with B y (cc 0.91) in the same manner as the Pi2 event. The PI and MI in B x are caused by ionospheric Hall currents driven by the dusk-to-dawn and dawnto-dusk electric fields, respectively, while those in B y are due to north-south Pedersen currents flowing from the polar to the equatorial ionosphere (Kikuchi et al., 2001). The MI of SC in B x is primarily composed of a stepwise increase caused by the magnetopause currents, superimposed by negative deflections due to the ionospheric Hall current in the morning sector (2130 UT, 0630 MLT) (Kikuchi et al., 2001). Since the Pedersen currents of the PI and MI are transmitted by the TM 0 mode waves (Kikuchi, 2014), the GIC is consistent with being the ground surface currents carried by the TM 0 mode waves. Figure 2 (right panels) shows the IEF in one-and two-layer models. The correlation of the GIC with E xI , cc(GIC − E xI ) 0.88, is less than the correlation with B y , cc(GIC − B y ) 0.91, but the correlation with E xII is extremely good as cc(GIC − E xII ) 0.99. The GIC can be reproduced almost perfectly by the two-layer model as shown with the solid and dotted curves in the frame of GIC. In the following sections, we examine the correlations for longer period/time scale disturbances, ranging up to 24 h.  (Sutcliffe and Lyons, 2002). Therefore, the fluctuations are quasimagnetostatic field of the substorm current wedge that can be calculated using the Biot-Savart formula, although the 30-min period is shorter than the typical recurrence period (1-2 h) of substorms (Akasofu, 1964;Freeman and Morley, 2004;Borovsky and Yakymenko, 2017).

Quasi-Periodic DP2 Fluctuations (30 min)
The correlation of the GIC with B y is cc(GIC − B y ) 0.54, and the correlation with B x is almost nothing as cc(GIC − B x ) −0.07. The cc(GIC − B y ) is much less than the previous event, probably because the fluctuations are superimposed by the background gradual increase that may not have affected the GIC. However, the correlations with E xI  , 1973). The positive and negative bays in B y are due to the location of MMB station being close to the upward and downward FACs in the pre-and postmidnight, respectively. The GIC (bottom left) resembles the B y in both events, and their correlation, cc(GIC − B y ) 0.67, is better than cc(GIC − B x ) 0.23 in the same manner as for the shortperiod disturbances, while the correlation is not so good as for the SC and Pi2. In contrast, the correlation with E xI is much better as cc(GIC − E xI ) 0.94, and furthermore, the two-layer model almost perfectly reproduces the GIC as cc(GIC − E xII ) 0.97. The change in sign of the GIC across the midnight again indicates that the GIC has no association with B x . The B y -dependence of the nighttime bay events raises a question on the TM 0 mode wave scenario since the magnetic bays on the nightside are caused not by ionospheric currents but primarily by field-aligned currents. It remains an issue what kind of propagation mode explains the B y dependence of the GIC on the nightside.  Figure 6 shows an example of the solar quiet geomagnetic variations (Sq). The period of Sq is 24 h, while significant changes occur over 8 h in the daytime (00-08 UT, 09-17 MLT). B x and B y are caused by the ionospheric currents driven by the thermospheric tidal motions (Kelley, 1989). It is remarkable that the correlation of the GIC with B y , cc(GIC − B y ) 0.27 is much lower than those for the shorter period disturbances. Furthermore, the correlation with B y is even less than the correlation with B x , cc(GIC − B x ) 0.66. The better correlation with B x does not necessarily mean that the GIC was caused by B x , since the temporal variations of the GIC resemble those of B y , if the time of B y is shifted ahead. On the other hand, the correlations with the IEF are extremely good as

DISCUSSION
The GIC in Hokkaido, Japan, can be reproduced from B y with high correlation coefficients as shown by Watari et al. (2009). We have further shown that the reproducibility strongly depends on the period of disturbances. As summarized in Table 1, the correlation with B y is high for short periods, e.g., SC (cc 0.91), but not for long periods, e.g., geomagnetic storm (cc 0.65) and Sq (cc 0.27). In particular, the correlation with B y of the Sq is even lower than the correlation with B x , cc(GIC − B x ) 0.66. Using the one-layer model composed of the semi-infinite uniform conductor with the flat surface of the ground, we have calculated the IEF, E x induced by B y . The IEF is found to be highly correlated with the GIC as cc 0.92 and 0.97 for the geomagnetic storm and Sq, respectively. This result implies that the long-period disturbances penetrated deep into the Earth, and the Earth can be considered to be uniform conductor. However, despite of the success with the one-layer model for long periods, the linear correlation with B y (cc 0.90) is still better than that with E xI (cc 0.77) for short period Pi2. This raises an issue on the period dependence of the reproducibility of the GIC.
To address this issue, we constructed the two-layer model composed of higher conductivity in the upper layer, following the TABLE 1 | Correlation coefficients between the GIC and surface magnetic fields, B x and B y and the electric field, E y and E x induced by B x and B y , respectively, calculated in one (I)-and two (II)-layer models for space weather (SW) events with periods ranging from 1min to 24 h. previous works on the geoelectric conductivity at Memambetsu (Owada, 1972). Fujii et al. (2015), using the MT method, clarified that the apparent resistivity of the Earth increases with an increasing period of geomagnetic disturbances at Memambetsu. This result is qualitatively consistent with the two-layer model with lower conductivity in the lower layer. The two-layer model with higher conductivity in the upper layer well explains the linear relationship with B (Pirjola, 2000). As summarized in Table 1, the induced electric field in the two-layer model, E xII , well reproduces the GIC with cc > 0.95 for both short-and long-period disturbances. To confirm the reproducibility with more events, we plotted the GIC and E xII for other nine events in Figure 7, where impulsive, periodic, isolated, and long-lasting disturbances on both the day and night sides are shown. The GIC (solid curve) well coincides with the E xII (dotted curve) that is scaled to the GIC. In particular, the peak of the GIC is well reproduced, which would help predict the GIC responsible for serious damages in the power transmission line.

SW events
We here check parameter dependence of the correlation coefficients (cc) in the two-layer model. Provided that the conductivities are fixed, major parameters responsible for cc are the number of reflections (n) and the depth of the upper layer (d 1 ) in the Eq. 22. Using the DP2 event in Figure 3, we calculated cc with d 1 20 km and different n. The cc increases as n increases such that cc 0.94, 0.96, 0.96, 0.97, 0.97 for n 20, 30, 40, 50, 100, respectively. We then calculated cc with fixed n 50 and different d 1 . The cc increases as d 1 increases such that cc 0.94, 0.96, 0.97, 0.97 for d 1 10, 15, 20, 30 km, respectively. Thus, we fixed n 50 and d 1 20 km in the two-layer model used for the calculation of the correlation coefficients.
Here, we make a brief comment on the singularity of 1/ t √ at t 0 included in the step response function. Love and Swidinsky (2014), Love and Swidinsky (2015) introduced the ramp function, t √ , derived from the inverse transform of 1/s s √ , to avoid the inconvenience in manipulating the singularity of 1/ t √ . By using the time difference of t √ , Love and Swidinsky (2015) reproduced the geoelectric field from the surface magnetic field in their twolayer model. The observed and calculated electric fields show fairly good coincidence, which may indicate success in using the ramp function. In our calculations, we replaced t with t + 0.0001 to avoid 1/ t √ ∞ at t 0. This approximation worked well to achieve the excellent correlations between the IEF and GIC, while it is just technical so that 0.0001 can be replaced with another small value.
We next examine if we can estimate the GIC that could have occurred during the past major storms. For this examination, we fix scale factors, k 1 (GIC/E x1 ) 8.0 [A/(mV/m)] and k 2 ( GIC/ E x2 ) 0.17 [A/(mV/m)], derived from the isolated substorm event in Figure 8. The observed GIC is well reproduced by both the one-layer and two-layer models with cc (GIC-E x1 ) 0.97 and cc (CIC-E x2 ) 0.99 as shown in the bottom left and right panels of Figure 8, respectively. For the sake of visual comparison, the observed GIC is plotted with dotted curves in each of the panels. Then, we used the scale factors, k 1 and k 2 , to reproduce the GIC observed during the SC event ( Figure 2). As shown in Figure 9, the GIC is well reproduced by the two-layer model with the same amplitude and high cc ( 0.99), whereas the GIC is not well reproduced by the one-layer model with lower cc ( 0.88) and overestimation of the rapid changes at the onset of the SC. The good correlation between the GIC and E x2 for the bay and SC events would allow us to use the scale factor k 2 to estimate the GIC for the past major storms. Figure 10 shows two examples of the estimated GIC during the storms on November 6, 2001 (panel (a)) and October 30, 2003 (panel (b)). It is interesting to note that the GIC estimated for the October 2003 storm has the largest magnitude at 20 UT because of the large magnitude of B y , when the storm ring current had not fully developed yet. This result suggests strong local time dependence of the GIC at MMB, which raises an important issue from the space weather forecasting point of view.
The parameters used in the two-layer model may not represent the ones estimated by the MT method (Fujii et al., 2015), but the excellent correlations in Table 1 allow us to use the model to reproduce the GIC from the observed magnetic field disturbances. Therefore, the model should be referred to as an empirical model that works for MMB. Although the model is not a commonly applicable model, we check the model with GICs measured at the Shin-Yamaguchi (SYG) substation of the Chugoku Electric Power Company in Yamaguchi prefecture in the western-southern part of Japan (34. 16°N,131.09°E GR;25.25°N,201.67°E GM). The power transmission line extends in the east-west direction along the coastline. Figure 11 shows a DP2 fluctuation event (T 80 min) observed at Kakioka (KAK,36.23°N,140.19°E GR;27.95°N,209.77°E GM) and the GIC at SYG, where the high frequency components are removed by applying the moving average over the window of 10 min. The KAK observatory is separated from SYG by 2.7°in GML, but the GIC is well correlated with the IEF such that cc (GIC-E xI ) 0.85 and cc (GIC-E xII ) 0.87. The model parameters are the same as used in the calculations for MMB except for the depth of the upper layer of the two-layer model being 15 km. Figure 12 shows an SC event with time scales of 1-10 min, where the window for the moving average is 30s. The correlation coefficients are better than those of the DP2 event such that cc (GIC-E xI ) 0.93 and cc (GIC-E xII ) 0.95. The correlation coefficients of the GIC with B x /E yII are not so good; cc 0.54/0.14 and 0.27/0.51 for the DP2 and SC events, respectively. It is remarkable that the GIC at SYG is strongly dependent on B y /E x , similarly to the GIC at MMB. Furthermore, there is no big difference between the one-and twolayer models, suggesting us to use the simple one-layer model to estimate the GIC at SYG during the past major storms. Using the models constructed in the present study, we would be able to predict the GIC during space weather disturbances with the aid of the global simulations. Ebihara et al. (2014) successfully reproduced ground magnetic disturbances due to the equatorial electrojet driven by the penetration electric fields during substorms. Furthermore, Tanaka et al. (2020) reproduced magnetic disturbances due to field-aligned currents as well as the ionospheric currents during the SC and substorm.
The power transmission line in Hokkaido is directed southwestward (Watari et al., 2009), which would predict that the GIC is affected equally by both B x and B y . However, as shown above, the GIC depends only on B y or E x . Furthermore, the GIC depends on the B y /E x at SYG, where the power line and coastline are in the east-west direction. Two possible mechanisms may explain the strong B y dependence. One is that the GIC is a ground surface current induced by the TM 0 mode waves in the Earthionosphere waveguide (Watari et al., 2009;Brändlein et al., 2012). The TM 0 mode wave transmits the B y and E x,z perpendicular and parallel to the (x-z) propagation plane, respectively, carrying the ionospheric currents and ground surface currents from high latitude to the equator (Kikuchi et al., 1978;Kikuchi, 2014). The TM 0 mode propagates at the speed of light and explains the simultaneous occurrence of the PI of SC (Araki, 1977) and DP2 fluctuations (Kikuchi et al., 1996) at high latitude and equator. In Figure 13, we show the equatorial electrojet (EEJ) defined as difference in B x between Yap, Micronesia (YAP, 0.5°GML), and Okinawa, Japan (OKI, 17.0°GML) (Kikuchi et al., 1996) together with B x (dots) and B y (solid) at MMB, KAK, and OKI. It is remarkable that the EEJ is well correlated with B y at MMB, KAK, and OKI, of which amplitude decreases as the latitude decreases. The latitudinal features may indicate that the Pedersen currents responsible for B y at middle latitudes flow into the equatorial ionosphere. Since the TM 0 mode waves induce ionospheric currents and ground surface currents (Kikuchi, 2014), it would be reasonable to attribute B y and the GIC to the TM 0 mode waves. It should be noted, however, that the TM 0 mode wave scenario may not be valid on the nightside, since ground magnetic fields are caused by magnetospheric currents in addition to the ionospheric currents during substorms (Ritter et al., 2008). Among these currents, the field-aligned currents transport the electromagnetic energy from the magnetosphere to the polar ionosphere. Therefore, a question arises, what kind of propagation mode transports the electromagnetic energy from the foot of the field-aligned currents or directly from the magnetosphere to the power transmission line at middle latitudes on the nightside? This will be a challenging issue of the magnetosphereionosphere coupling at middle latitudes.
FIGURE 12 | (A) B x and B y observed at KAK in the morning (2220 UT, 0720 MLT) and GIC at SYG during the SC event (T 1-10 min) with the same parameters as in Figure 11, except that the GIC data is smoothed over 30s. (B) IEFs in the same format as in Figure 1. The other possible mechanism is the effects of the geometry such as the direction of power lines and coastlines and of the 3-D structures of Earth's conductivities (Fujii et al., 2015;Goto, 2015;Nakamura et al., 2018). Ivannikova et al. (2018) found that much of Great Britain was affected by coastal effects owing to the strong conductivity gradient between the land and the ocean. The coastline effects on the GIC are also significant in Hokkaido as deduced from the model calculations (Nakamura et al., 2018). Furthermore, Fujii et al. (2015) clarified that the MT response at Memambetsu shows that B y affects the induction in x-direction more strongly than B x does in y-direction. The MT-deduced anisotropy is explained by means of the spatial inhomogeneity of the Earth's conductivity. Thus, both effects of the coastline and inhomogeneous distribution of the Earth's conductivity should have affected the anisotropic response of the GIC to the surface magnetic field. We would need to take into account the inhomogeneous conductivity distribution even in a thin layer model (e.g., McKirdy and Weaver, 1984). However, the horizontally uniform models employed in the present study well explain the GIC-B y /E x correlations. The consistency between the MT and GIC results and the inconsistency between the nonuniform and uniform models remain a question to be addressed in future studies. CONCLUSION 1) We have shown that the GIC at Memambetsu in Hokkaido (35.7°GML) is linearly correlated with the y-component geomagnetic field, B y , for the short-period disturbances such as the geomagnetic sudden commencements (cc 0.91) and Pi2 pulsations (cc 0.90), while the correlation was found to become worse as the period of disturbances increases, such that cc 0.67 for the substorm and cc 0.27 for the solar quiet diurnal variations. 2) The induced electric field in the one-layer model with the semi-infinite conductor (σ 10 −4 mho/ m) well reproduces the GIC with cc 0.94 for the substorm and 0.97 for the solar quiet variations. But, the correlation with B y (cc 0.91) is still better than the correlation with the induced electric field (cc 0.88) for short period, SC. 3) We constructed the two-layer model with higher conductivity in the upper layer (σ 1 10 −4 , σ 2 10 −8 mho/ m), which is found to be capable of reproducing the GIC with high correlations for both short periods, cc 0.99 for SC, and for long periods, cc 0.95 for Sq. 4) The GIC at Shin-Yamaguchi, Japan (25.3°GML) is well correlated with B y /E xII similarly to the GIC at MMB, such that cc 0.87 and 0.95 for the DP2 and SC events, respectively. 5) The strong B y dependence of the GIC could be associated with the TM 0 mode in the Earth-ionosphere waveguide, which carries B y and ionosphere-ground surface currents from high latitude to the equator. This mechanism should be valid on the dayside, but it remains an issue to explain the B y dependence on the nightside, where the magnetospheric current effects dominate over the ionospheric current effects.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
TK built the one-and two-layer models and wrote the whole sections of the manuscript. YE contributed to creating the GIC database and to the calculation of the IEF in the models. KKH performed space weather event studies with magnetometer data and contributed to establishing the data acquisition system from SYG substation. KK installed the GIC meter and calibrated the data at SYG. S-IW contributed to recording the GIC at MMB. All authors contributed to manuscript revision, read, and approved the submitted version.