The M7.0 Japan earthquake occurred on April 15, 2016, in Kumamoto region (Kyushu Island) in Japan. According to Japan Meteorological Agency (JMA) the magnitude is even larger (MJMA 7.3). Its epicenter was located at 32.76°N, 130.76°E, and its hypocentral depth was 12 km. The earthquake was preceded by a large foreshock occurred on 14 April near the location of the incoming mainshock.



Figure 1: Earthquake epicenter and Dobrovolsky Area.


This event occurred as the result of strike-slip faulting at shallow depth. Focal mechanisms for the earthquake indicate slip occurred on either a left-lateral fault striking to the northwest, or on a right-lateral fault striking northeast (USGS).The shallow depth and faulting mechanism of this earthquake indicate it occurred on a crustal fault within the upper Eurasia plate.

At the location of this event, the Philippine Sea plate converges with Eurasia towards the northwest at a velocity of 58 mm/yr. Deep earthquakes in the region are related to the subduction of the Philippine Sea plate. Shallow seismicity is less frequent and of lower magnitude then the deep one. Thirteen M5+ earthquakes have occurred at shallow depths (> 50 km) within 100 km of the April 2016 seismic event over the past century.




  1. Accelerated Moment Release (AMR) & Revised-AMR

The Catalogue earthquake analysis for this event started from January 2010 with a maximum distance of 1100 km from the epicentre. In the next step of the study a temporal cut was applied starting from May 2012 in order to exclude the large Tohoku earthquake (11 March 2011, M9). Figure 2 shows AMR results with the new catalogue starting from May 2012.


Figure 2: AMR curve corresponding to Catalogue starting from May 2012.

R-AMR method was applied in search of the optimum parameters summarized in Figure 3. It is evident the capability of this technique to better disclose acceleration covered by the far events.


Figure 3: Although R-AMR applied on Catalogue is already able to highlight the acceleration in data (a), several trials for seeking the best parameters lead to the best set (thick green border).





Two different approaches were developed to search EQ precursor using geomagnetic Swarm data: MASS method and Wavelet method. Both analyses are based on Level 1B MAGxLR Swarm product.


  1. MASS algorithm (Magnetic Swarm anomaly detection bySpline analysis)

The algorithm MASS (Magnetic Swarm anomaly detection by Spline analysis) was applied to M7.0 Japan 2016 earthquake with different thresholds (kt=2.0, 2.5 and 3.0), while the moving window was fixed at 3.0°. The algorithm analyses all tracks in DbA one month before and one after the EQ. The tracks are marked as “anomaly” only if the centre of the moving window is in DbA and if geomagnetic conditions are quiet. Figure 4 shows an example of anomalous tracks detected by MASS method for Japan EQ.


Figure4:Example of anomalous track in Y magnetic component (Satellite A- March 21, 2016).

The cumulative number of anomaly tracks detected by MASS one month before and after the Japan EQ (threshold  kt= 2.5) is shown in Figure 5. For the Y magnetic field component the number of anomalies detected after the EQ is bigger than that detected before (nEQ = 2.47) and this results seems to exclude the interpretation of the data as EQ-related phenomenon. It is also possible to describe the cumulative number of anomalies with S-shape function, but the critical point is after the EQ and probably not due to the EQ.


Figure5: Cumulative number of anomaly detected by MASS one month before and one month after the M7.0 Japan EQ. Threshold is kt = 2.5  and the anomalies are selected only with geomagnetic quiet time (|Dst| ≤ 20 nT and ap≤ 10 nT).



  1. Wavelet Analysis


The Wavelet spectral analysis has evidenced the existence of different anomalous families, each characterized by some particular featuresthat altogether do not clarify whether they are linked to LAIC or not.

Figure 6 shows a type of anomaly appearing as a short oscillation, i.e. spatially confined: even in absence of seismicity this kind of anomaly is frequent, but unfortunately not only inside the DbA.



Figure 6: Wavelet analysis results few days before the Japan M7.0mainshock (March 25, 2016).


Instead, Figure 7shows a different type of anomaly, which extends to a wide part of the tracks, thus involving a large sector of the ionosphere. Again, their presence even in regions far from the DbA does not help in clarifying their origin.


Figure7: Both these anomalies involve a large part of the tracks, although spatially confined as a packet; in addition, they appear both inside and outside the DbA.



  1. SWARM – IONOSPHERE from Satellite

Satellite-based data for the ionospheric characterization of the EQ-related events are mainly those referred to the LP (Langmuir Probe) instrument aboard the SWARM satellites. The electron density Ne is the relevant parameter used for the ionospheric characterization in the frame of SAFE project. Two different methods were developed to analyze Swarm ionospheric data: NeSTAD and NeLOG.


  1. Method I: NeSTAD

The NeSTAD analysis has been applied to the Swarm constellation data (LP and IBI data )available in the period from 16 MArchj to 7 May 2016 have been used to derive the track anomaly parameters.The NeSTAD has been initialized with the mild outliers mode (k=1.5) and with an “excess area” parameter equal to 0.1.

Then, to tag the interesting track anomalies for this particular event, the following criteria have been applied:

  • R>Rthr=0.7 and standard deviation of the filtered track below σthr=0.1 or standard deviation of the filtered track above σthr=0.1 independently of R value.
  • Only night time tracks have been selected (18-06 LT), because during night time and at mid latitudes the ionosphere is expected to be less turbulent.
  • Quiet ionospheric conditions (absolute value of Dst in the considered day not exceeding 20 nT).

An example of tagged anomaly is provided in figure 8.


Figure 8: Identified anomaly with the NeSTAD algorithm and tagged as interesting for the Japan M7.0 Japan EQ event. Tagged anomaly refers to Swarm Alpha satellite on 10 April 2016.



Figure 9: - Cumulative number of anomalies identified through the NeSTAD algorithm for Swarm satellite Bravo.


Figure 9 shows the cumulative number of the tagged anomalies through the NeSTAD algorithm for Bravo satellite. The black dashed line indicates the day in which the Japan M7.0 EQ event occurred. The red boxes indicate the days in which disturbed geomagnetic conditions have been recorded and that are not included in the tagging process. The two blue dashed lines indicate the time interval during which the given Swarm satellite covered the DbA of the event at night time. The purple crosses identify the tagged anomaly that are not considered as good, because likely due to EIA/plasma bubbles.


  1. Method II: NeLOG

The automatic search NeLOG of ionosphere electron density anomaly from Swarm data is applied to Japan M7.0 case study. A sample is classified as anomaly if the residual value exceeds a threshold kt times the RMS of the residual after the spline fit. A track is selected as an “anomaly” if it has more than 10 anomaly samples in DbA and if geomagnetic indices are |Dst|<20 nT and |Aap|<10 nT. Tracks are selected within a mean local time between 22 and 6.



Figure10: Anomaly detected by NeLOG (satellite Alpha, 21 March 2016) before the M7.0 Japan EQ.

Figure 10 shows an example of anomaly track 25 days before in the date of the mainshock, while Figure 11 reports the cumulative number of anomaly windows detected by NeLOG. For this event, R2 is not compatible with a linear fit, so the anomalies have not a random distribution over time. The factors nEq and C after normalization are equal 1.00 so they are not significant. The apparent positive response of the non-normalised factors is due in part to satellite coverage and in part to geomagnetic disturbances.


Figure11:Cumulative number of anomaly samples detected by NeLOGone month before and one month after M7.0  Japan EQ. Threshold is kt = 2.5, the anomalies are selected only with geomagnetic quiet time (|Dst| ≤ 20 nT and ap ≤ 10 nT) and in night time (22 ≤ LT ≤ 6).


  1. IONOSPHERE Ground-based


By means of the analysis on the ionosonde data (Kokubuni Ionospheric observatory, Latitude: 45.2° N, Longitude: 141.8° E) three anomalies have been identified, occurred on 5 April, 10 April and 11 April. As an example, the 5 April anomaly will be shown (Figure 12) and discussed.

During 5 April 2015, the geomagnetic conditions were quiet for the entire UT day. The maximum value of ap was 7nT and AE was low in the first part of the day, when the ionospheric anomaly occurred. Thus, in correspondence with the occurrence of the 5 April anomaly quiet conditions recorded, so it seems that such anomaly is not related to geomagnetic or auroral activity and can be a candidate for LAIC-related anomaly.

Similar analysis allows to conclude that both the anomalies occurred on 10 April and 11 April 2016 could be a candidate for LAIC-related event.


Figure 12: The ionospheric anomaly of the 5 April 2016; Δh’Es, δfbEs, and δfoF2 variations and 3-hour ap.




The “single-station GNSS analysis” is applied to vTEC values from data acquired by different IGS GNSS stations located in (or just outside) the area of Dobrovolsky from 9 April to 15 April 2016 (see Figure 13).


Figure13: IGS GNSS stations used to perform the single station GNSS analysis. Red point represents the EQ epicenter and blue circle shows the DbA.

Considering each station separately, the EEMD (Ensemble Empirical Mode Decomposition)is applied to obtain the IMF composing the vTEC signals. Then, the statistical test µ/σ is applied to separate the IMF contributions to the signal trend and oscillation.The CWT is then used to extract the time/frequency representation of the oscillation component. Figure 14 shows the wavelet spectrum of the oscillating signal from 9 April to 16 April for each of the stations considered.


Figure 14: Continuous Wavelet Spectrum of the oscillation component of vTEC over shao (panel a), aira (panel b), smst (panel c).

As for the EQ occurring in Japan on 16 February, for this case study the single station technique reveals a relationship between the anomalies detected and the distance from the epicenter.


The M7.0 Japan EQ was analyzed using both Swarm geomagnetic and ionospheric data searching for earthquake-related anomalies in the frame of LAIC theory.

Results from NeLOG method shows that the R2value obtained is not compatible with a linear fit, so the anomalies have not a random distribution over time. The factors nEq and C after normalization are equal 1.00 so they are not significant. In addition, NeLOG algorithm found a concentration of anomalies about one month before the EQ, but they probably are not due to LAIC because many of them have been detected also by Bravo satellite.

 The cumulative number of anomalies derived by tagging procedure on the NeSTAD output on Swarm satellites allowed identifying 3 anomalies for Alpha, 4 for Bravo and 2 for Charlie. The anomalies tagged for Alpha and Charlie occurred few days before the event, while those identified for Bravo are meaningfully earlier. Moreover, according to the Dst criterion used by the NeSTAD algorithm, many days before and during the event have been found to be disturbed, limiting the possibility to identify possibly LAIC-related anomalies.

Finally, the result of the GNSS analysis reveals a relationship between the anomalies detected and the distance from the epicenter. Two anomalies have been detected on:

  • 10 April from 03.00 UT to 14.00 UT, corresponding to a frequency of about 0.04 mHz (period≈7 hours) and visible on the spectrogram relative to SHAO and AIRA
  • 3 April from 01:00 UT to 05:00 UT, corresponding to a frequency of about 0.08 mHz (period≈3.5 hours).