A reassessment of the anomalous geomagnetic diurnal variation two months prior to the 2011 Tohoku earthquake, Japan (Mw 9.0)

This study reassesses the claim of an anomalous geomagnetic diurnal variation two months prior to and northwest of the eventual epicenter of the 2011 Tohoku earthquake, Japan (Mw 9.0). Instead of discussing the ratios of geomagnetic diurnal variation ranges at observation and reference stations, which previous studies have analyzed, the ranges themselves are examined in the present work. For the reassessment, the technique to identify the range of typical variations is improved, leading to more appropriate values of the geomagnetic diurnal variation range than those obtained by previous studies. The variations that were deemed anomalous by previous authors can be explained by coupling of the ranges of three components of the geomagnetic diurnal variation; this coupling effect was not discussed in previous studies. The cause of the high ratio was a combination of low-amplitude quiet-day diurnal solar variations in the geomagnetic field in winter, day-to-day variations, and a weak magnetic disturbance during the period. It is therefore not strictly correct that the observed variations were anomalous and, as such, they cannot be unambiguously related to seismo-magnetic phenomena due to crustal activity prior to the megathrust earthquake. In addition, due to issues with the statistical treatment of samples in previous studies, a more appropriate statistical discussion of the matter is presented.


Introduction
Daily variations in the geomagnetic field when solar-terrestrial disturbances are absent are called solar quiet daily variations (Sq), whose characteristics were described in detail by Chapman and Bartels [1940]. Their amplitudes depend on the latitude and local-time longitude relative to the meridian through the sun, and increases with the sunspot number. Daily changes in Sq are much greater during daylight hours than during darkness; Sq also undergoes considerable annual variations, which particularly affect its amplitude, and some day-to-day variability in Sq has been found. The direction of change in Sq in the horizontal plane at mid-latitudes in the northern hemisphere rotates clockwise from north to east and upward during the daytime. In extra-tropical latitudes, the range of daily variations in each element is greatest in summer and lowest in winter.
Sq is regarded as being caused by both the ionospheric currents, which are external to the Earth, and Earth induction, which is internal [Campbell 2003]. Spherical harmonic analysis of geomagnetic data can separate Sq into components with external and internal origins. For Sq in the mid-latitudes, Campbell [2003] also showed that while the horizontal Sq with external and internal origins are parallel, their vertical components are anti-parallel. It is also well known that in the frequency domain, the vertical component of the geomagnetic field variation is linearly related to the two horizontal components of the variations, depending on the electrical resistivity structure surrounding the observation point and the frequency [see, for example, Vozoff 1991].
The geomagnetic diurnal variation range, hereafter called GDVR, can be defined as the difference between the maximum and minimum values of the northward (+X), eastward (+Y), and downward (+Z) components of the geomagnetic field for each day, taken from the records at each geomagnetic station. It is expected that, on quiet days, the GDVR is dominated by Sq. Based on statistical analyses, Xu et al. [2013], Han et al. [2015Han et al. [ , 2016, and Han et al. [2017] claim that the GDVR decreased prior to the 2011 off the Pacific coast of Tohoku Earthquake, Japan (Mw 9.0), which occurred at 05:46:18, March 11, 2011 UT with its epicenter 38°06.2'N, 142°51.6'E according to unified earthquake catalog by Japan Meteorological Agency ( JMA). The earthquake is hereafter referred to as the 2011 Tohoku earthquake. It is claimed that two months preceding the occurrence of the megathrust earthquake, a decrease was found only in the vertical component of GDVR at geomagnetic stations northwest of the epicenter, at distance approximately 180 km from the point of rupture initiation. Since they do not use GDVR itself, but instead compute the ratio of GDVR at a reference station divided by the corresponding value at a target station (not its reciprocal), the magnitude of the change of the GDVR itself is unclear. Han et al. [2016] ascribes the likely cause to the electrical resistivity change of the crust around the stations close to the epicenter, due to pre-seismic strain changes. By examining the spatial and temporal relationships between the GDVR decrease and geodetic data collected close to the epicenter, Han et al. [2016] concluded that the coupling of multiple pre-earthquake phenomena may help to understand the preparation process of mega earthquakes in subduction zones.
However, despite its insensitivity to long-term DClike tendencies in geomagnetic field variations, GDVR is sensitive to various factors that affect geomagnetic variations on a variety of timescales. Geomagnetic phenomena, such as magnetic storms, sudden impulses, bays, and solar flare effects, can perturb GDVR values. In addition, Xu et al. [2013], Han et al. [2015Han et al. [ , 2016, and Han et al. [2017] primarily focus on the vertical component of the GDVR. The variation of the vertical component couples with horizontal component variations due to geomagnetic phenomena which are dominant in the horizontal components. The coupling coefficients in the frequency domain depends on both the period and the electrical resistivity structure surrounding the geomagnetic station. A careful investigation of relative vertical component GDVR ratios is therefore required, especially with respect to locality.
This study aims to investigate basic features of GDVR that were not discussed by previous studies. As such, the spatial and temporal distribution of GDVR in Japan in early January 2011 are carefully reassessed. This work shows that temporal changes in GDVR in the epicentral area can be explained by geomagnetic phenomena without invoking seismo-magnetic phenomena by anomalous crustal activity.

Original data
Time series of geomagnetic data sampled in one-minute intervals at the 17 stations listed in Table  1 are analyzed in this study. The national network of geomagnetic stations in Japan is operated by the Geospatial Information Authority of Japan (GIAJ) and the JMA. Their locations, along with the reference location of the epicenter of the 2011 Tohoku earthquake, are shown in Figure 1. The original time series data are in D (declination), H (horizontal component), and Z in UT. The X (positive north) and Y (positive east) components can then be calculated as X = H cos D and Y = H sin D, respectively. Since only the definitive data are analyzed in the present study, the period of data analyzed depends on the status of data distribution for each station; these are also listed in Table 1. The total number of days from 1 January 1997 to 31 December 2015 during which geomagnetic data at ESA, MIZ, KNZ, MMB, KAK, and KNY are analyzed is 6939. The resolution of H and Z is 0.1nT. The Q component (where Q = X, Y, or Z) at the geomagnetic station whose code is ABC is hereafter referred to as ABCQ.

Filtering using wavelets
To eliminate short-period disturbances and noises in the X, Y, and Z components of the original 1-minute geomagnetic data, a discrete wavelet transform is used to filter the data as a low-pass filter. A length-10 Daubechies wavelet is convolved with the original time series, after Jach et al. [2006] based on the max-  Table 1 and the epicenter of the 2011 Tohoku earthquake (cross). imum overlap discrete wavelet transform. Wavelet coefficients at level j ≤ 5, which are most sensitive to energy at periods less than or equal to 64 min, are set to 0.0. The filtering procedure up to this point is the same as that of Xu et al. [2013], Han et al. [2015Han et al. [ , 2016, and Han et al. [2017]. In the present study, to evaluate the GDVR more precisely, the same thresholding procedure is applied to eliminate wavelet coefficients at transform levels j ≥ 11, which localize energy at periods longer than 2048 min. This part of the procedure thus behaves as a high-pass filter. Consequently, this technique applies a bandpass filter with corner periods of 64 and 2048 min to the original time series data. Figure 2 shows an example of the time series of filtered ESAX, ESAY, and ESAZ on quiet days.

Treatment of missing data
Missing data in the original time series are treated as follows. In the original time series, data gaps that consist of either isolated one-minute missing samples or sequences less than or equal to 6 min are padded by linear interpolation. Otherwise, an anomalous value 99999.0nT is substituted for each missing datum. Modified time series are obtained by these interpolation and substitution rules. A filtered modified time series is then obtained by applying the wavelet filter to the modified time series. Similarly, a corresponding dummy time series is generated with substituted data set to 1.0, and 0.0 otherwise. After applying an identical wavelet filter to the dummy time series, if the absolute value of the filtered dummy time series is larger than 10 -5 , the datum in the filtered modified time series at the corresponding time is regarded as influenced by the data missing value(s), and treated as missing. The threshold of 10 -5 corresponds to the assumption that the influence of missing data is equivalent to a tolerance of up to 1nT.   Figure 3 shows typical time series of geomagnetic data at stations ESA and KAK in a period from 00:00 2 January 2007 JST ( Japan Standard Time, UT+9) to 09:00 8 January 8 2007. The dates 2 and 3 January 2007 (UT) are two of the international five disturbed days, while 7 January 2007 (UT) is one of the international five quietest days. According to magnetic storm and bay disturbance catalogs by Kakioka Magnetic Observatory (KMO), JMA, no magnetic storms were identified during this period, while bays were recorded on 2, 3, 4, and 5 January. KMO identifies and reports the begin times and peaks of bays of ranks A and B, and the magnitude and polarity of each component; for the bays of rank C, only the begin time and polarity are identified and reported. A class B bay on 3 January 2007, identified by KMO, started at 22:31 JST. The peak field strengths at KAKX, KAKY, and KAKZ were estimated as +44, +15, and +22nT, respectively, at 23:51, 23:47, and 23:52 JST, respectively. On this day, the maximum value of KAKZ in the filtered time series was identified at 23:47 JST, close to the peak time (23:52 JST) at KAKZ. The difference of 5 min is not inappropriate, considering the characteristics of a bandpass filter with a corner period of 64 and 2048 minutes. The maximum value at ESAZ was attained at 03:19 JST and is apparently unrelated to the bay. Two cross-correlation coefficients between ESAX and KAKX, and between ESAY and KAKY, during the period shown in Figure 3, are 0.984 and 0.990, respectively. ESAX and KAKX show high positive correlation to one another, so do ESAY and KAKY. The obtained GDVR of ABCQ is hereafter expressed as R(ABCQ). It is therefore expected that R(KAKX) ∝ R(ESAX) and R(KAKY) ∝ R(E-SAY). In contrast, ESAZ and KAKZ are similar only if the geomagnetic field variations are quiet. On days on which bays are identified, maxima at KAKZ coincide with the timing of peaks of positive bays, while the minima occur just around noon, as on quiet days. ESAZ appears sensitive to abrupt change at ESAY, and therefore KAKY. A negative correlation appears to exist between ESAY and ESAZ, and therefore KAKY and ESAZ. Cross-correlation coefficient between ESAY and ESAZ is -0.166, while the coefficient between KAKY and KAKZ is -0.094. As a result, R(KAKZ) is expected to be sensitive to abrupt increases in KAKX and inflated, with a similar large value for R(KAKX). R(ESAZ) is expected to be sensitive to abrupt decreases in KAKY or normal Sq. In general, it should be remembered that the GDVR of the Z component is influenced by not only many geomagnetic factors, but also the resistivity structure surrounding the geomagnetic station. It is therefore obvious that the potential causes of temporal changes in an index defined as, for example, R(KAKZ)/R(ESAZ) are also many. R(KAKX) and R(KAKY) can be used to discuss the GDVR of the Z component as reference data. K-indices based on data at KAK (hereafter expressed as K(KAK)), the period of magnetic storms and the list of bays, all of which are reported by KMO are also referred to. The original method to derive K(KAK) is reported by Yumura [1951]. The revision of the scale of K indices, and the correlation of K(KAK) with Kp indices are assessed by Uesugi et al. [2005].

GDVR in local time
GDVR for each component at each station is determined for each day in JST in the present study, unlike Xu et al. [2013], Han et al. [2015Han et al. [ , 2016, and Han et al. [2017], which discuss GDVR in UT. Based on the description of Sq by Chapman and Bartels [1940], it is natural to define daily maxima and minima of geomagnetic variations corresponding to identical Sq current vortices in the ionosphere and the Earth, which move westward with the sun due to Earth's rotation. Because of the definition of JST, GDVR on solar quiet days in Japan is given by maximum and minimum values for an identical daytime in JST, while it is not guaranteed when using UT (see Figure 2). When discussing geomagnetic diurnal variations in Japan in UT, the determined GDVR can use maximum and minimum Sq values of current vortices that skip a night, and therefore are influenced by temporal changes in Sq current vortices on timescales of half a day. In order to discuss the components in GDVR that might possibly arises from crustal activity, or the locality of the GDVR, it is necessary to minimize the influence of the temporal changes in GDVR due to other possible sources, such as temporal changes in Sq current vortices due to changes in solar activity, which can be a cause of the day-to-day variability of Sq [Chapman and Bartels 1940]. For geomagnetic data recorded in Japan, adopting JST instead of UT is necessary for this specific reason. Figure 4 shows the times for each day when maximum or minimum values are recorded at ESAZ. R(E-SAZ) on quiet days is mainly characterized by the minimum. In UT, however, minima are often identified at 00:00 or 23:59 mainly in winter and summer. This is because, in Japan, 00:00 UT corresponds to the local morning, when Sq is already close to its peak; the exact Sq value, of course, varies day-to-day and seasonally. The identified maximum or minimum in a UT calendar day is therefore under-or over-estimated if the date changes in the middle of the daily geomagnetic variation, which is much greater during the hours of sunlight according to Chapman and Bartels [1940]. The numbers of days on which the minimum of ESAZ is identified at 00:00 and 23:59 are 206 and 1020, respectively, using UT; in JST, these values are 11 and 14, respectively. Thus, by adopting UT days, as in Xu et al. [2013], Han et al. [2015Han et al. [ , 2016, and Han et al. [2017], R(ESAZ) could be under-or over-estimated for as much as 18% of samples. Since it is natural to discuss the daily variations of the geomagnetic field based on local time, after Chapman and Bartels [1940], JST is used in the present study for further discussions of GDVR. The GDVR on a given day is regarded as missing in the present study when a maximum or minimum is identified at 00:00 or 23:59. In case where missing data in the filtered time series are found on a given day, the GDVR on that day is also treated as missing.   netic storms identified by KMO are eliminated from these plots. As a reference time series, the (JST) daily sum of K(KAK), expressed as Σ'K(KAK), is shown. The following four observations can be made from Figure 5.

Assumption of coupling among GDVRs
(1) The X components at all stations change similarly with Σ'K(KAK), especially with respect to abrupt changes that recover within several days and a long-term change possibly related to the period of 11-yr sunspot cycle.
(2) An annual change in the Y-component data is prominent.
(3) Variations in the horizontal components are similar from station to station.
(4) For ESA and MIZ, the Z-component is positively correlated with Y; for KNZ and KAK, Z is positively correlated with both X and Y (see Table 2).
Based on the features of the GDVRs shown in Figure 5, the following linear relationship is assumed for R(ABCQ), using KAK as the reference station: (1) The model coefficients A and B and the offset constant C are constants for each component of each station. A, B, and C can be estimated with an ordinary least-squares technique to minimize the sum of squares of the residual ΔR(ABCQ) ≡ R(ABCQ) -r(AB-CQ), where r(ABCQ) is the estimated R(ABCQ).

Averaging of daily ranges
For comparison with the results by Han et al. [2015Han et al. [ , 2016, and Han et al. [2017], a 15-day running average is calculated. This 15-day window was chosen by Han et al. [2015Han et al. [ , 2016 and Han et al. [2017] according to Liu et al. [2006] and Kon et al. [2011], claiming that it has been frequently used to obtain reliable background variations in previous ionospheric studies. The date of the center of each time window is given as the timestamp of the averaged value in the present study, while in Han et al. [2016] a backwards-running window was adopted. In the present study, in cases where eight or more observations in the original 15-day GDVR window are missing, the averaged value for the given date is treated as missing. In case where seven or fewer data from a given 15-day window are missing, the computation uses a weighted average, with a weight of 0 for missing data and 1 otherwise. No thresholding with respect to Σ'K(KAK) is applied to obtain the running average in the present study; Han et al. [2016] applied thresholding by excluding data on days when there were one or more Kp indices of > 5 determined and reported by Space Weather Prediction Center (SWPC), National Oceanic and Atmospheric Administration, USA, during the day, while by Han et al. [2015] > 6. Notice that during the period shown in Figure 3, the Kp indices are at most 5. The method to derive Kp by SWPC is based on Takahashi et al. [2001], assuming that Kp values derived by this method can be an appropriate proxy for Kp derived from worldwide geomagnetic data.
The averaged value of R(ABCQ) is hereafter expressed as R (ABCQ). Similar to Equation (1), are the coefficients and C is the offset constant for R (ABCQ), respectively. The residual ΔR (ABCQ) ≡ R (ABCQ) -r (ABCQ), where r (ABCQ) is the estimated R (ABCQ), is treated as the index of the anomalous value for R (ABCQ). The sum of squares of ΔR (ABCQ) is minimized to obtain Ā ,B _ , and C . To determine whether or not unique anomalous geomagnetic diurnal variations occurred in the period from 1 to 15 January 2011, Ā , B _ , and C were determined for a data set in which the samples during this period are excluded.
Notice that R (ABCQ) is mainly discussed in the present study with R (KAKX) and R (KAKY), based on the nature of GDVRs which appears in the geomagnetic data shown above. Xu et al. [2013], Han et al. [2015Han et al. [ , 2016, and Han et al. [2017] discussed averaged R(KAKQ)/R(ABCQ). This discrepancy is discussed after discussing the features of R (ABCQ) with respect to its relationships with R (KAKX) and R (KAKY). Figure 6 shows ΔR (ESAZ), ΔR (MIZZ), ΔR (KNZZ), and ΔR (KAKZ) from 1 January 1997 through 31 December 2015 without thresholding GDVRs with respect to Σ'K(KAK). Data recorded during the magnetic storm determined by KMO with the data at KAK are eliminated from procedures to determine constants and residuals. KMO monthly reports magnetic storm data including the date and time in UT of each storm. The period of a mag- netic storm is defi ned in the present study as the period in JST between the beginning and end date and time. A GDVR datum on a day on which magnetic storms begin, continue, or end is regarded to be infl uenced by magnetic storms and eliminated in the analyses. The number of analyzed samples varies slightly from station to station, because of missing data and those due to cases where 15day running averaged GDVR cannot be determined. The coeffi cient of determination R 2 obtained in determining Ā , B _ and C is between 0.69 and 0.77. Red symbols show values in the objective period between 1 and 15 January 2011. During this period, the ranges of ΔR (ESAZ) and ΔR (MIZZ) are -2.6 to -1.4nT and -2.2 to -0.8nT, respectively. The absolute values of the ranges are < 1s for each station, where s is the standard deviation. The coefficients and offset constants were determined without the GDVR during the objective period.

GDVR residuals: observed vs. expected values
The values of ΔR (ESAZ) and ΔR (MIZZ) in the objective period are as small as one standard deviation at each station and the values of ΔR (KNZZ) and ΔR (KAKZ) in the same period. This result suggests that R (ESAZ) and R (MIZZ) are well-explained, with R (KAKX) and R (KAKY) equal to R (KNZZ) and ΔR (KAKZ). The periods when either ΔR (ESAZ), ΔR (MIZZ), ΔR (KNZZ), or ΔR (KAKZ) exceed 3s are listed in Table 3. The anomalous periods when the residuals exceed 3s are grouped into three.
(a) The periods which appear immediately preceding or after a magnetic storm: the time window is 7 days, a half of 15 days. The length equals to the time window of averaging.
(b) The periods which appear simultaneously at all stations.
(c) The periods which appear at KNZ and KAK separately or simultaneously.
The anomalous periods at ESA and MIZ belong to either the fi rst or the second groups.

Features of model coeffi cients and offset constants
The analyses applied to the data from four geomagnetic stations can be applied similarly to data from other sites. Figure 7 shows the distribution of Ā , B _ and C obtained for the GDVR of the Z component of each station. Ā , B _ and C obtained with data   sets which include the data during the period from 1 to 15 January 2011 do not significantly differ from those obtained with the data sets which exclude the data during the period, considering their estimation errors. Seventeen geomagnetic stations are classified into three groups: Ā > B , Ā = B , and Ā < B . The relation Ā > B , at a station suggests that, at the station in question, the GDVR of the Z component is more sensitive to R(KAKX), while Ā < B is more sensitive to R(KAKY). Figure 7 also shows that, for the stations with Ā > B , C , values tend to be larger than at stations where Ā < B . As a concluding remark of these results, it should be remembered that stations ESA and MIZ with MMB are in a subgroup of stations for which the intensity of coupling between the Z-com-ponent GDVR and R(KAKX) and R(KAKY), defi ned as (Ā 2 + B 2 ) 1/2 is smallest; the subgroup belongs to a group with Ā < B . Figures 8 and 9 show the residuals at each station in the groups Ā < B and Ā > B , respectively. Again, there is no anomalous increase or decrease during the objective period from 1 to 15 January 2011. Generally, the standard deviations of the residuals are smaller for the group with Ā < B than for the group with Ā > B , possibly due to the insensitivity of the Z-component of GDVR to R(KAKX), which is sensitive to geomagnetic disturbances.

Cause of anomalously high R(KAKZ)/R(ESAZ) and R(KAKZ)/R(MIZZ)
I conclude that there is no anomalous decrease in R(ESAZ) and R(MIZZ) in the period between 1 and 15 January 2011, based on considerations of their cou-  pling, mainly with R(KAKY). Xu et al. [2013], Han et al. [2015Han et al. [ , 2016, and Han et al. [2017] neglected the coupling of R(ESAZ), R(MIZZ), and R(KAKZ) with R(KAKX) and R(KAKY), which could be one of the reason why they regarded R(KAKZ)/R(ESAZ) or R(KAKZ)/R(MIZZ) as anomalous during the same period. Xu et al. [2013] and Han et al. [2017] discussed the small amplitude of ESAZ and MIZZ during the objective period as causes of anomalously high averaged R(KAKZ)/R(ESAZ) and R(KAKZ)/R(MIZZ). Due to the fact that Sq varies day-to-day (Chapman and Bartels, 1940), and the finding that R(ESAZ) and R(MIZZ) couple mainly with R(KAKY), it is better to examine the horizontal components of Sq and its variations during this period more carefully. Figure 10 shows the geomagnetic variation at KAK, BMT (Beijing Ming Tombs, China), CNB (Canberra, Australia), SPT (San Pablo-Toledo, Spain), and HER (Hermanus, South Africa) during a period that includes the objective period. The time series data are filtered with the bandpass technique described above. The geodetic latitudes of these stations are 36. 232°N, 40.300°N, 35.314°S, 39.547°N, and 34.425°S, respectively. Qualitatively, the decrease in amplitudes of the horizontal component of Sq (X or Y) in early January 2011 (through 6 January 2011, UT) was global. After a geomagnetic disturbance on 6 January 2011, the amplitudes began to recover. The decrease in R(ESAZ) and R(MIZZ) in this period can be related to corresponding global decreases in amplitudes of the horizontal component of Sq via the coefficient B .
Although a characteristic spatial pattern can be found for the distributions of geomagnetic stations with Ā > B or Ā < B , it is necessary to discuss this with respect to resistivity structures on spatial scales as large as the induction scale length with a period of 1 day and its overtones, which principally comprise the GDVR. However, a discussion of the resistivity structures surrounding each geomagnetic station does not require GDVR. As was already described, since GDVR is an index influenced by many factors that affect geomagnetic variation, it is inadequate to explain GDVR in terms of resistivity structures alone. If discussion of resistivity changes due to strain changes prior to a mega-earthquake is intended, then a direct discussion of resistivity structures estimated with electromagnetic soundings (e.g., by magnetotelluric methods) would be far more appropriate.

Thresholding using geomagnetic disturbances
Another cause of the reported increase in averaged R(KAKZ)/R(ESAZ) and R(KAKZ)/R(MIZZ) in early January 2011 was the misleading use of thresholding by Han et al. [2016] when computing running averages of R(KAKZ)/R(ESAZ). For example, SWPC reports that Kp = 5 for 21:00-24:00, 6 January 2011 (UT). The other Kp indices on that day were 0 or 1. On 7 January 2011 (UT), the Kp value was 4 for 00:00-03:00 and between 1 and 3 for the rest of the day. There were no other days in January 2011 in which Kp > 3. Based on the thresholding by Han et al. [2015] or Han et al. [2016], in which geomagnetic data are excluded from the days when one or more Kp indices are greater than 6 or 5, respectively, none of the data in January 2011 are excluded from the running averaging. It is possible to change the threshold of Kp and confirm the effects on the result, the procedure of which is important to show the robustness of the original result, and therefore to confirm whether or not an anomalous GDVR is present.
Instead of such a confirmation, the present study simply confirms the contribution of the GDVR on 7 January 2011 JST by treating the corresponding datum as missing, because the period from 21:00 on 6 January to 03:00 on 7 January in UT belongs to 7 January 2011 JST. As references, results with a datum on 6 or 8 on January 2011 ( JST) treated as missing are compared. The comparison is a type of procedures to extract an outlier from a group of samples, which was also possible in Xu et al. [2013] and Han et al. [2015]. The present comparison is a test directly related to the result of Xu et al. [2013], which does not adopt thresholding in averaging, with differences in the lengths of running average time windows taken into consideration. The time series computed after running averaging of R(KAKZ)/R(ESAZ) from 1 December 2010 to 28 February 2011 ( JST) are shown in Figure 11. It is clearly shown that the datum on 7 January 2011 JST is an outlier, which enhances the averaged R(KAKZ)/R(ESAZ) from 31 December 2010 to 14 January 2011 by roughly a factor of two. It is also shown that the averaged R(KAKZ)/R(E-SAZ) between 31 December 2010 and 14 January 2011 JST have large standard errors in both the original time series and the time series in which the datum on 6 or 8 January 2011 are treated as missing. For the time series in which the datum on 7 January is treated as missing, the standard errors between 31 December 2010 and 14 January 2011 JST are smaller than those of other time series described above by roughly 50%. This analysis simply suggests that the datum of 7 January 2011 JST is an outlier. For the findings of Han et al. [2016], I now carefully investigate whether the anomaly in early January 2011 is statistically significant. Han et al. [2016] describes the space weather during the period from 1 to 21 of January 2011 as quiet except for a slight storm on 7 January ( JST), and claims that the anomaly is not likely to have resulted from a slight storm. On the contrary, the present test suggests that the anomaly was produced, or at least enhanced, by a so-called "slight" storm. The present investigation suggests that this anomaly is not due to decrease in the GDVR at ESAZ, but rather, an increase in the GDVR at KAKZ accompanying a weak disturbance with high Ā for R(KAKZ). The day-to-day and seasonal variability of Sq emphasized the abnormality of averaged R(KAKZ)/R(ESAZ) during the objective period by decreasing R(ESAZ). The present investigation also suggests that the index defined by the ratios of GDVRs of the Z components at different stations is easily influenced by a disturbance with Kp = 5, and therefore unreliable as an indicator of localized geomagnetic field variations, particularly those caused by crustal activities. Han et al. [2016] claims that there were a number of minor storms from 1997 to 2012, that the anomaly appears only in January 2011, and that the anomaly likely did not result from such a minor storm. However, the question of whether or not the corresponding minor storm caused the anomalous R(KAKZ)/R(ESAZ) has already been discussed. A magnetic disturbance manifests in many ways; Kp is one possible means of parameterizing it. In cases where seismo-magnetic phenomena are identified in geomagnetic data, which are influenced by magnetic disturbances, a careful discussion of the disturbances, paying attention to more parameters than Kp indices alone, is indispensable. The present study makes another simple interpretation: that the weak disturbance on 7 January 2011 was one of the rare geomagnetic disturbances that could cause an anomalous averaged R(KAKZ)/R(ESAZ).

Normality of the disturbance of 7 January 2011
It is natural to discuss the intensity of the geomagnetic disturbance in the context of amplitude variations of horizontal component data. In the present study, this is quantified by the GDVRs of the horizontal components. On 7 January 2011 JST, R(KAKX) and R(KAKY) were 54.2nT and 26.0nT, respectively; R(KAKZ) and R(ESAZ) were 29.5nT and 5.4nT, respectively. Based on the above discussion of Ā , B , and C for R (ESAZ) and R (KAKZ), R(KAKX) > 54.2nT and R(KAKY) < 26.0nT can cause higher R(KAKZ)/R(ESAZ) ratios. The days on which these two conditions are satisfied are hereafter referred to as anomalous days; here, the data are searched for such anomalies. The number of anomalous days, whether quiet or subject to disturbances, is eight including 7 January 2011 ( JST). The date, R(KAKX), R(KAKY), R(KAKZ), R(ESAZ), Kp indices, and Σ'K(KAK) on these anomalous days are listed in Table 4. It is noteworthy that all the anomalous days are in winter, when the amplitude of Sq is smaller: December, January, or February. Among the anomalous days, 26 January 2006 (UT) is the only day including a Kp index larger than 5, which means it was excluded from the analysis by Han et al. [2016] due to their thresholding criterion. The time se- ries of geomagnetic field variation on 3 January 2007 JST, another day amongst the anomalous days, is shown in Figure 3. Since a bay is a semi-global phenomenon, R(KAKZ) on this day can be overestimated from the perspective of discussing the locality of GDVR, while R(ESAZ) cannot. Thus R(KAKZ)/R(ESAZ) can be overestimated. The number of anomalous days for which calculations of R(KAKZ)/R(ESAZ) cannot be discarded by considering geomagnetic phenomena is therefore reduced from eight to six. However, large R(KAKZ)/R(E-SAZ) can have more causes than magnetic disturbances. It is possible to check the abnormal characteristics of the weak disturbance during the objective period, paying attention to the GDVRs of all the components. In contrast, the present study also finds that the minor disturbance of 7 January 2011 JST was typical. The minor disturbance caused a difference between the observed and expected GDVR values of the Z component. However, this difference was smaller than 3s limit by considering the coupling among the GDVRs of all three components.

Examination of statistical discussion
The first study of vertical component GDVR values to claim an anomalous increase two months prior to the 2011 Tohoku earthquake was Xu et al. [2013]. The present study shows that the anomaly was due to a disturbance on 7 January 2011 JST, not due to randomness of the data. Xu et al. [2013] also showed a similar feature of non-randomness, but derived a different conclusion that the anomalous GDVR can be related to the occurrence of the 2011 Tohoku earthquake. However, their statistical discussion did not demonstrate an anomalous R(KAKZ)/R(ESAZ) in January 2011. Their statistical discussion is inappropriate as follows. Xu et al. [2013] began their discussion with the statistical hypothesis that the anomalous increase was due to random sampling and running averaging of the ratio of R(KAKX)/R(ESAX), R(KAKY)/R(ESAY), and R(KAKZ)/R(ESAZ) during the period from 1 January 2010 to 31 December 2012. The significance level was set to 0.05. A stochastic test was performed and the probability of appearance of the anomaly was measured. An anomalous datum was defined in their work as a sample that exceeded the threshold µ ± 3s, where µ and s are the mean and standard deviation, respectively. As an anomaly was found with probabilities for averaged R(KAKX)/R(ESAX), R(KAKY)/R(ESAY), and R(KAKZ)/R(ESAZ) of 0.026, 0.001, and 0.013, respectively, these probabilities were below the significance level of 0.05 and there-fore the hypothesis of randomness was rejected. Although the assumed probability distribution was not explicitly declared, a threshold of 3s coincides with a significance level of 0.005 in a case where a Gaussian is assumed, while a level of 0.05 corresponds to a threshold of 2s. If the numerical results of their stochastic test are correct, then the hypothesis cannot be rejected at 3s: 0.026 for R(KAKX)/R(ESAX) and 0.013 for R(KAKZ)/R(ESAZ); only for R(KAKY)/R(E-SAY) can the null hypothesis be rejected. However, Xu et al. [2013] claimed that the non-randomness of R(KAKY)/R(ESAY) was due to seasonal variations. It is finally concluded that the clear anomaly before the mega earthquake in R(KAKZ)/R(ESAZ) is not a random anomaly but a reliable one with a statistical significance.
The statistical test of Xu et al. [2013] is implicitly based on an assumption that the ratio of GDVR obeys a normal distribution; however, Anderson-Darling test ( [Stephens 1986] can check the normality of the samples. To test this hypothesis in the present study, sample values of R(ESAZ) and R(KAKZ) from 1 January 2010 to 31 December 2012 UT are low-pass wavelet filtered to follow the GDVR preprocessing of Xu et al. [2013]. Since the samples on days with missing data at ESA are treated as missing in the present study, including a one-day window around the missing datum and the dates 1 January 2010 and 31 December 2012, the resultant number of samples is 1074 for both R(ESAZ) and R(KAKZ). Applying the test to the obtained R(ESAZ) and R(KAKZ) values, the Anderson-Darling statistics for goodness of fit of a normal distribution are 6.237 and 10.961, respectively. Since the hypothesis of normality is rejected if the statistic exceeds 1.159 at a significance level of 0.005, the hypotheses that R(ESAZ) and R(KAKZ) obey the normal distribution are rejected. The test statistics for the hypothesis that data follow a lognormal distribution are, on the other hand, 1.080 and 0.854 for R(ESAZ) and R(KAKZ), respectively. Assuming a lognormal distribution is natural because the samples are defined as positive, and a certain amount of skewness is found in the distribution of the samples, as shown in Xu et al. [2013]. Since the threshold to reject the hypothesis of a lognormal distribution is also 1.159 at the 0.005 significance level, the hypotheses that R(ESAZ) and R(KAKZ) obey lognormal distribution cannot be rejected. The ratio R(ESAZ)/R(KAKZ) thus appears to obey a lognormal distribution, not a normal distribution. It is possible to discuss the distribution of bivariate samples (R(ESAZ), R(KAKZ)) based on a bivariate lognormal distribution [Mostafa and Mahmoud 1964], since its marginal distributions are lognormal. Discussing the probability of obtaining a given ratio is then possible after, for example, Pham-Gia et al. [2006]. The present study only suggests the possibility, and follows the procedure above to determine the correct statistical distribution of the ratios.
The distributions of R(ESAZ), R(KAKZ), R(KAKZ)/R(ESAZ), and averaged R(KAKZ)/R(ESAZ) computed in 10-day windows as in Xu et al. [2013], can also be discussed in a similar manner to the above. It is concluded that Xu et al. [2013] does not demonstrate a statistically anomalous value of averaged R(KAKZ)/R(ESAZ) two months prior to the 2011 Tohoku earthquake, because the parameters are incorrect and the normality itself is not proven. Alternatively, the present study clearly shows the inaccuracy of discussing the ratio from its nature.

Conclusions
The claims by Xu et al. [2013], Han et al. [2015Han et al. [ , 2016, and Han et al. [2017] that unique GDVR values were observed on the vertical components of two geomagnetic stations two months prior to the 2011 Tohoku earthquake Japan (Mw 9.0) are reassessed. The present study shows that the observed variations can be explained by considering their coupling with the GDVRs of horizontal components at the reference station KAK. The apparent unique variations in the ratios of vertical component GDVR val-ues are caused by coupling due to the season, when the amplitude of Sq, which characterizes GDVR, is small, with contributions from a weak geomagnetic disturbance. At ESA and MIZ, R(ESAZ) and R(MIZZ) mainly couple with the GDVRs of Y-component data, which are less sensitive to geomagnetic disturbances. Since the amplitude of Sq is small in winter and was small globally in early January 2011, so were R(ESAZ) and R(MIZZ) during this period. At KAK, R(KAKZ) couples, more than R(ESAZ) and R(MIZZ) do, with the GDVR of the X component, which is more sensitive to geomagnetic disturbances. As a result, R(KAKZ)/R(ESAZ) and R(KAKZ)/ R(MIZZ) become large due to a weak geomagnetic disturbance that occurred on 7 January 2011 ( JST), in winter when the amplitude of Sq was small, and in a period when the amplitude of Sq decreased on a global scale. Since Xu et al. [2013], Han et al. [2015Han et al. [ , 2016, and Han et al. [2017] discussed such ratios without considering the coupling, it can be concluded that their interpretations are not robust.
By considering the coupling among the GDVRs of all three components at each station, statistically significant anomalies in the GDVRs of Z-component data at ESA and MIZ are not found two months prior to the mega earthquake. It is therefore not necessary to discuss seismo-magnetic phenomena appearing as GDVRs due to a possible preparation process of the mega earthquake in the manner in which Han et al. [2016] and Han et al. [2017] discuss. The conclusion of this study are derived from the high-pass filtering   (0)) and on the following day (Kp(1)); and Σ'K(KAK) on each day. "-" means that a GDVR value could not be obtained. Kp(0) and Kp(1) are reorganized to correspond to times in JST (UT+9).
of geomagnetic data for precise discussion of GDVR, computation of GDVR using windows based on local time, the assumption of coupling of geomagnetic data, and the statistical discussion of the samples.