Earthquake Report: Hokkaido, Japan

Following the largest typhoon to strike Japan in a very long time, there was an earthquake on the island of Hokkaido, Japan today. There is lots on social media, including some spectacular views of disastrous and deadly landslides triggered by this earthquake (earthquakes are the number 1 source for triggering of landslides). These landslides may have been precipitated (sorry for the pun) by the saturation of hillslopes from the typhoon. Based upon the USGS PAGER estimate, this earthquake has the potential to cause significant economic damages, but hopefully a small number of casualties. As far as I know, this does not incorporate potential losses from earthquake triggered landslides [yet].

This earthquake is in an interesting location. to the east of Hokkaido, there is a subduction zone trench formed by the subduction of the Pacific plate beneath the Okhotsk plate (on the north) and the Eurasia plate (to the south). This trench is called the Kuril Trench offshore and north of Hokkaido and the Japan Trench offshore of Honshu.

The okhotsk plate is considered part of the North America plate on some maps. The location of the plate boundary of the Okhotsk plate are not well understood (e.g. using GPS plate motion velocities, it is difficult to find the northern boundary with the North America plate).

Many of the earthquakes in this region are related to the subduction zone. Most notably is the 2011 Tohoku-oki M 9.1 tsunamigenic earthquake. More background information about the 2011 earthquake can be found here and information about the tsunami can be found here.

The 2011 earthquake had lots of aftershocks and was quite complicated. One interesting thing that happened is that there was an extensional earthquake in the Pacific plate to the west of the Japan Trench. This M 7.7 earthquake happened along faults formed as the Pacific plate bends near where it meets the trench. Similar subduction zone / outer rise earthquake pairs are known, including some along the New Hebrides Trench in the western equatorial Pacific ocean, as well as further north along the Kuril subduction zone. I spend time discussing the 2006/2007 Kuril earthquake pair in this report.

There was also a subduction zone earthquake in 2003, the Tokachi-oki earthquake, that triggered submarine landslides. These landslides transformed into turbidity currents and these were directly observed with offshore instrumentation.

One of the interesting things about this region is that there is a collision zone (a convergent plate boundary where two continental plates are colliding) that exists along the southern part of the island of Hokkaido. The Hidaka collision zone is oriented (strikes) in a northwest orientation as a result of northeast-southwest compression. Some suggest that this collision zone is no longer very active, however, there are an abundance of active crustal faults that are spatially coincident with the collision zone.

Today’s M 6.6 earthquake is a thrust or reverse earthquake that responded to northeast-southwest compression, just like the Hidaka collision zone. However, the hypocentral (3-D) depth was about 33 km. This would place this earthquake deeper than what most of the active crustal faults might reach. The depth is also much shallower than where we think that the subduction zone megathrust fault is located at this location (the fault formed between the Pacific and the Okhotsk or Eurasia plates). Based upon the USGS Slab 1.0 model (Hayes et al., 2012), the slab (roughly the top of the Pacific plate) is between 80 and 100 km. So, the depth is too shallow for this hypothesis (Kuril Trench earthquake) and the orientation seems incorrect. Subduction zone earthquakes along the trench are oriented from northwest-southweast compression, a different orientation than today’s M 6.6.

So today’s M 6.6 earthquake appears to have been on a fault deeper than the crustal faults, possibly along a deep fault associated with the collision zone. Though I am not really certain. This region is complicated (e.g. Kita et al., 2010), but there are some interpretations of the crust at this depth range (Iwasaki et al., 2004) shown in an interpreted cross section below.

I present more about the basics behind ground shaking, triggered landslides, and possible earthquake triggering on Temblor here:

Below is my interpretive poster for this earthquake

I plot the seismicity from the past month, with color representing depth and diameter representing magnitude (see legend). I include earthquake epicenters from 1918-2018 with magnitudes M ≥ 6.5 in one version.
I plot the USGS fault plane solutions (moment tensors in blue and focal mechanisms in orange), possibly in addition to some relevant historic earthquakes.
I also include active crustal faults from the Coordinating Committee for Geoscience Programmes in East and Southeast Asia (CCOP). Note the abundance of north-northwest oriented yellow lines to the east of today’s earthquakes. While today’s earthquake was not on those crustal faults, the earthquakes and these faults are responding to similarly oriented tectonic stresses.

  • I placed a moment tensor / focal mechanism legend on the poster. There is more material from the USGS web sites about moment tensors and focal mechanisms (the beach ball symbols). Both moment tensors and focal mechanisms are solutions to seismologic data that reveal two possible interpretations for fault orientation and sense of motion. One must use other information, like the regional tectonics, to interpret which of the two possibilities is more likely.
  • I also include the shaking intensity contours on the map. These use the Modified Mercalli Intensity Scale (MMI; see the legend on the map). This is based upon a computer model estimate of ground motions, different from the “Did You Feel It?” estimate of ground motions that is actually based on real observations. The MMI is a qualitative measure of shaking intensity. More on the MMI scale can be found here and here. This is based upon a computer model estimate of ground motions, different from the “Did You Feel It?” estimate of ground motions that is actually based on real observations.
  • I include the slab contours plotted (Hayes et al., 2012), which are contours that represent the depth to the subduction zone fault. These are mostly based upon seismicity. The depths of the earthquakes have considerable error and do not all occur along the subduction zone faults, so these slab contours are simply the best estimate for the location of the fault.

    Magnetic Anomalies

  • In the map below, I include a transparent overlay of the magnetic anomaly data from EMAG2 (Meyer et al., 2017). As oceanic crust is formed, it inherits the magnetic field at the time. At different points through time, the magnetic polarity (north vs. south) flips, the north pole becomes the south pole. These changes in polarity can be seen when measuring the magnetic field above oceanic plates. This is one of the fundamental evidences for plate spreading at oceanic spreading ridges (like the Gorda rise).
  • Regions with magnetic fields aligned like today’s magnetic polarity are colored red in the EMAG2 data, while reversed polarity regions are colored blue. Regions of intermediate magnetic field are colored light purple.
  • Note the parallel magnetic anomalies to the east of Japan. These were formed about 150 million years ago at the spreading center where this portion of the Pacific plate was created. More can be found about the creation of the Pacific plate in Boschman and van Hinsbergen, 2(016).

    I include some inset figures. Some of the same figures are located in different places on the larger scale map below.

  • In the upper right corner is a low angle oblique view of the tectonic configuration in this region. Note how many subduction zones are that are interacting in different ways. This is from the AGU blog, “Trembling Earth.” I place a blue star in the general location of today’s earthquakes (same for other figures in this poster).
  • In the lower right corner is a plate tectonic map of this part of the world (Liu et al., 2013). The major plate boundary faults are shown, along with the volcanoes in the magmatic arcs. Also, seismicity is shown (the 2011 earthquake as a small blue star) and the slab contours for the Pacific and Philippine Sea plates. Color shows the age of the oceanic crust. These authors place the southern boundary of the Okhotsk plate further to the south (dashed black line), where the Izu Collision Zone intersects Japan (near the intersection of the magmatic arc associated with the Izu-Bonin Trench, with Japan).
  • In the lower left corner is a geologic map of Japan (van Horne et al., 2016). Note the orientation of the rocks in Hokkaido as they are oriented in a northwest-southeast direction in the area labeled Hidaka Collision. These rocks are oriented this way due to the northeast-southwest convergence. This map places the southern boundary of the Okhotsk plate near where the Hidaka Collision is. Compare this with the Liu map to the right.
  • In the upper left corner is a large scale portion of a figure from NUMO (Kurikami et al., 2009), a publication put together by the N to evaluate the suitability of sites for high level radioactive waste. They considered various geologic hazards in this report. This map shows some key tectonic features and geologic data. I include the legend to the right of the map. The magmatic arc is shown as a red line. The Hidaka Collision Zone is shown as a dashed blue line with arrows showing the direction of collision. The blue arrows show the direction of maximum stress, the stress field. These arrows are pointed in the direction of compression. The convergence direction along the collision zone is oriented well with today’s earthquakes, but the stress field data are not perfectly oriented.
  • Here is the map with a month’s seismicity plotted.


    •

  • Here is the map with a centuries seismicity plotted.


    •

Some Relevant Discussion and Figures

  • This map shows the current tectonic configuration of this region, along with some inherited features from the tectonic past (e.g. green lines). This is from NUMO’s report: “Evaluating Site Suitability for a HLW Repository (Scientific Background and Practical Application of NUMO’s Siting Factors), NUMO-TR-04-04.”


    •

  • Also from the NUMO report, this shows the Niigata-Kobe fold and thrust belt. In addition, this map shows a northwest striking convergent plate boundary along the southeastern boundary of Hokkaido. However, it cannot explain the interesting orientation of the M 6.2 deep (240 km) earthquake.


    •

  • Here is a great figure from Itoh et al. (2005) that shows how they interpret the Hidaka Collision Zone.


    • Maps showing tectonic context around the Japanese Islands (a) and geologic belts in Hokkaido (b; after Kato et al., 1990).

  • This map (also from Itoh et al., 2005) shows the active faults and folds mapped in the region, along with the geology.


    • Geologic map around the Umaoi anticline redrawn from Geological Survey of Japan (2002). Location of active fault and/or fold scarps (after Ikeda et al., 2002) are also shown. buQ and bdQ attached on fault traces are upthrown and downthrown sides of faults, respectively. Sampling points of surface paleomagnetic data is after Kodama et al. (1993).

  • Here is more evidence for the thrust faults associated with the Hidaka Collision Zone (Iwasaki et al., 2004). These authors used seismic refraction and seismic reflection experiments to interpret the deep crustal structures associated with the collision here. The profile shown in the next figure is denoted by the east-west oriented black arrows in the lower part of this figure.


    • Geological map of Central Hokkaido with our seismic refraction/wide-angle reflection profiles and shot points (stars). Seismic reflection lines of the Hokkaido Transect were laid out from shot L-2 to M-5 on the wide-angle line. Reflection lines carried out from 1994 to 1997 in the southernmost part of the HCZ and refraction/wide-angle reflection lines in 1984 and 1992 are also shown. SYB: Sorachi-Yezo Belt; KMB: Kamuikotan Metamorphic Belt; IB: Idon’nappu Belt; HMB: Hidaka Metamorphic Belt; HB: Hidaka Belt; YB: Yubetsu Belt; TB: Tokoro Belt; HMT: Hidaka Main Trust.

  • Here is the interpreted cross section from Iwasaki et al. (2004). Note (1) the thrust faults and (2) the depths for these different structures. There are still regions that are poorly understood. Recall the depth of the M 6.6 earthquake is about 33 km.


    • Geological interpretation of the seismic model. KMB: Kamuikotan Metamorphic Belt; IB: Idon’nappu Belt; HMB: Hidaka Metamorphic Belt; Yz: Yezo Super Group; Sr: Sorachi Group; HMT: Hidaka Main Thrust.

  • Here is the cool tectonic map from Liu et al. (2013). We all like cool maps! (right?)


    • Tectonic settings of the study region (black box). The solid sawtooth lines and the black dashed line denote the plate boundaries (Bird 2003). The red triangles denote the active volcanoes. The blue dashed lines and the pink lines denote the depth contours to the upper boundary of the subducting Pacific slab and that of the subducting Philippine Sea slab, respectively (Hasegawa et al. 2009; Zhao et al. 2012). The topography data are derived from the GEBCO_08 Grid, version 20100927, http://www.gebco.net. The ages of oceanic plates are from M¨uller et al. (2008).

  • This is a very cool figure (also from Liu et al., 2013) that shows a plot of earthquakes from 3 different perspectives. First is the map view. To the right of the map is a plot of earthquakes shown as viewed from the east of the map and this shows the hypocenters. The profile below the map shows a cross section of seismicity as viewed from the south looking north. The original figure includes more maps (A and B).


    • (c) Distribution of the 4803 earthquakes used in
    this study. The black crosses denote 3818 events (Group-1) that occurred under the seismic network. The green dots show 228 events (Group-2) that occurred outside the seismic network, selected from the events relocated by Gamage et al. (2009) using sP depth phases. The red dots denote 757 suboceanic earthquakes (Group-3) that are newly relocated in this work using P-wave, S-wave and sP depth-phase data. (d) East–west and (e) north–south vertical cross-sections of the earthquakes shown in (c).

  • Here is a map from the recent update of the Japan National Seismic Hazard Maps, resulting from knowledge gained following the 2011 M 9.1 earthquake (Fujiwara et al., 2012). The color represents the chance that a region will experience ground shaking at or greater that Japan Meteorological Agency (JMA) seismic intensity 6 in the next 30 years. JMA intensity is a scale of shaking intensity similar to the Modified Mercalli Intensity (MMI) Scale. The numbers are different, so they are difficult to compare. The JMA intensity 6 is similar to MMI X. Today’s earthquakes are in a region of slightly elevated chance of ground shaking (between 6-26%). Today’s M 6.6 earthquake may have reached


  • This is a map from the National Research Institute for Earth Science and Disaster Resilience, where data from the Strong-motion Seismograph Networks in Japan are located. This shows measurements of JMA intensity. It appears that a site near the epicenter (red star) reached JMA intensity 7.


  • This is an animation from the same source showing observations of JMA intensity recorded at the surface throughout Japan. h/t to Jascha Polet for sharing this on twitter.

  • Here is the upper figure showing the tectonic setting (Kurikami et al., 2009). Note how the Okhotsk plate has a strike-slip fault that terminates near the Hidaka Collision Zone (called a forearc-sliver fault, formed because the plate convergence is oblique to the subduction zone fault). I include their figure caption as a blockquote.


    • Tectonic setting of Kyushu within the Japanese island arc. The locations of active faults and volcanoes that have been active in the last 10,000 years are also shown.

  • This is a fantastic educational video from IRIS that discusses the plate tectonics and mentions some earthquakes in the region of Japan.

  • Here is a USGS poster than summarizes the earthquake history and plate geometry for this region. This is the USGS Open File Report 2010-1083-D (Rhea et al., 2010).


    •

Earthquake Triggered Landslides

  • Here is the aerial video from NHK that shows some of the landslides triggered by this sequence of earthquakes today. This comes from a tweet below.
  • Well, here is a great figure from Keefer (1984) that shows that the larger the magnitude of an earthquake, the larger an area can be subject to triggering of landslides from the ground shaking from that earthquake.


    • Area affected by landslides in earthquakes of different magnitudes. Numbers beside data points are earthquakes listed in Table 1. Dots = onshore earthquakes; x = offshore earthquakes. Horizontal bars indicate range in reported magnitudes. Solid line is approximate upper bound enclosing all data.

  • In 2008 there was an earthquake in China with a magnitude M 7.9. Unfortunately this earthquake caused many deaths. Using satellite imagery, geologists identified about 60,000 individual landslides (Gorum et al., 2011). Below is a map that shows the faults in the region, as well as epicenters from the earthquakes from this sequence.


    • Location and 12May 2008Wenchuan earthquake fault surface rupturemap, and focalmechanisms of the main earthquake (12May) and two of the major aftershocks (13 May and 25 May). Also the epicenters of historic earthquakes are indicated. The following faults are indicated: WMF: Wenchuan–Maowen fault; BF: Beichuan–Yingxiu fault; PF: Pengguan fault; JGF: Jiangyou–Guanxian fault; QCF: Qingchuan fault; HYF: Huya fault;MJF:Minjian fault based on the following sources: (Surface rupture: Xu et al., 2009a,b; Epicenter and aftershocks: USGS 2008; Historic earthquakes: Kirby et al., 2000; Li et al., 2008; Xu et al., 2009a,b).

  • This map shows the region where there was a high density of landslides (Fan et al., 2012). Note how the majority of landslides are located near the larger earthquakes (the larger circles in the above map).


    • Distribution of landslide dams triggered by the Wenchuan earthquake, China. The high landslide density zone is defined by a landslide area density >0.1 km−2; also shown are epicenters of historical earthquakes (USGS, 2008) and the historical Diexi landslide dams (Dahaizi, Xiaohaizi and Diexi). White polygons are unmapped due to the presence of clouds and shadows in post-earthquake imagery. WMF: Wenchuan–Maowen fault; YBF: Yingxiu–Beichuan fault; PF: Pengguan fault; JGF: Jiangyou–Guanxian fault; QCF: Qingchuan fault; HYF: Huya fault; MJF: Minjiang fault (after X. Xu et al., 2009). MJR: Minjiang River; MYR: Mianyuan River; JJR: Jianjiang River; QR: Qingjiang River.

  • Many of these landslides dammed rivers, causing an additional hazard. These earthen dams block rivers, leading to a large lake forming upstream of these dams. The dams can be overtopped when the lakes fill with water. once the water reaches the top of the dam, they can overflow and rapidly down cut back to the level of the river prior to the dam emplacement. If this happens too rapidly, a flood can occur, putting those downstream at risk of flooding.


    • Comparison of densities of blocking and non-blocking landslides. (a) Landslide density. (b) Landslide dam point density. White dashed lines are 240-km by 25-km swath profiles. (c). Mean normalized landslide and landslide dam densities along the SW–NE profile. Red lines are Yingxiu-Beichuan fault (YBF) and Pengguan fault (PF). Yellow dash lines are the boundary of the P1–P7 watersheds in the Pengguan Massif. YX, WC, HW, BC, and QC are the cities of Yingxiu, Wenchuan, Hanwang, Beichuan and Qingchuan, respectively. MJR, JJR, FJR, and QR represent Minjiang, Jianjiang, Fujiang and Qingjiang rivers, respectively.

  • In 1959, there was an earthquake in southwestern Montana, the M 7.2 Hebgen Lake Earthquake. This earthquake triggered a landslide that dammed the Madison River. This dam created a lake now called “Earthquake Lake.” I was actually driving on a road trip following my graduation from Oregon State University in 2014. I drove to this area and arrived the day that the Earthquake Lake Visitor Center opened. Pretty cool.
  • Here is a view of the lake as it was in May, 2014. Note the dead trees. The landslide is the bare looking mountainside in the distance on the left. We are looking to the West.


  • Here is a view of the landslide from my truck.


  • Here are all the people waiting to go into the visitor center on opening day.


  • Here is another cool view of the ghost forest.


  • Here is an educational display near the lake. Click on the image and one may zoom in within their browser, or save the image and zoom in that way. The text is readable if one wants to follow along.


  • This is from the poster and shows the landslide dam after it formed.


    •

Geologic Fundamentals

  • For more on the graphical representation of moment tensors and focal mechnisms, check this IRIS video out:

  • Here is a fantastic infographic from Frisch et al. (2011). This figure shows some examples of earthquakes in different plate tectonic settings, and what their fault plane solutions are. There is a cross section showing these focal mechanisms for a thrust or reverse earthquake. The upper right corner includes my favorite figure of all time. This shows the first motion (up or down) for each of the four quadrants. This figure also shows how the amplitude of the seismic waves are greatest (generally) in the middle of the quadrant and decrease to zero at the nodal planes (the boundary of each quadrant).


    •

  • Here is another way to look at these beach balls.

  • There are three types of earthquakes, strike-slip, compressional (reverse or thrust, depending upon the dip of the fault), and extensional (normal). Here is are some animations of these three types of earthquake faults. The following three animations are from IRIS.

    • Strike Slip:

    Compressional:

    Extensional:

  • This is an image from the USGS that shows how, when an oceanic plate moves over a hotspot, the volcanoes formed over the hotspot form a series of volcanoes that increase in age in the direction of plate motion. The presumption is that the hotspot is stable and stays in one location. Torsvik et al. (2017) use various methods to evaluate why this is a false presumption for the Hawaii Hotspot.


    • A cutaway view along the Hawaiian island chain showing the inferred mantle plume that has fed the Hawaiian hot spot on the overriding Pacific Plate. The geologic ages of the oldest volcano on each island (Ma = millions of years ago) are progressively older to the northwest, consistent with the hot spot model for the origin of the Hawaiian Ridge-Emperor Seamount Chain. (Modified from image of Joel E. Robinson, USGS, in “This Dynamic Planet” map of Simkin and others, 2006.)

  • Here is a map from Torsvik et al. (2017) that shows the age of volcanic rocks at different locations along the Hawaii-Emperor Seamount Chain.


    • Hawaiian-Emperor Chain. White dots are the locations of radiometrically dated seamounts, atolls and islands, based on compilations of Doubrovine et al. and O’Connor et al. Features encircled with larger white circles are discussed in the text and Fig. 2. Marine gravity anomaly map is from Sandwell and Smith.

    References:

  • Chapman et al., 2009. Development of Methodologies for the Identification of Volcanic and Tectonic Hazards to Potential HLW Repository Sites in Japan –The Kyushu Case Study-, NUMO-TR-09-02, NOv. 2009, 192 pp.
  • Fan, X., et al., 2012. Transient water and sediment storage of the decaying landslide dams induced by the 2008 Wenchuan earthquake, China in Geomorphology, v. 171-172, p. 58-68, doi:10.1016/j.geomorph.2012.05.003
  • Fujiwara, H., Morikawa, N., Okumura, T., Ishikawa, Y., and Nojima, N., 2012. Revision of Probabilistic Seismic Hazard Assessment for Japan after the 2011 Tohoku-oki Mega-thrust Earthquake (M9.0) in Proceedings of the 15th World Conference on Earthquake Engineering, 15th World Conference on Earthquake Engineering, Lisbon.
  • Gorum, T., Fan, X., van Westen, C.J., Huang, R., Xu, Q., Tang, C., Wang, G., 2011. Distribution pattern of earthquake-induced landslides triggered by the 12 May 2008 Wenchuan earthquake in Geomorphology, v. 133, p. 152-167, doi:10.1016/j.geomorph.2010.12.030
  • Hayes, G. P., D. J. Wald, and R. L. Johnson, 2012. Slab1.0: A three-dimensional model of global subduction zone geometries in J. Geophys. Res., 117, B01302, doi:10.1029/2011JB008524.
  • Itoh, Y., Ishiuyama, T., and Nagasaki, Y., 2005. Deformation mode in the frontal edge of an arc–arc collision zone: subsurface geology, active faults and paleomagnetism in southern central Hokkaido, Japan in Tectonophysics, v. 395, p. 81-97 doi:10.1016/j.tecto.2004.09.003
  • Iwasaki, T., et al., 2004. Upper and middle crustal deformation of an arc–arc collision across Hokkaido, Japan, inferred from seismic refraction/wide-angle reflection experiments in Tectonophysics, v. 388, p. 59-73, doi:10.1016/j.tecto.2004.03.025
  • Keefer, D.K., 1984. Landslides caused by earthquakes in Geological Society of America Bulletin, v. 95, p. 406-421, doi: 10.1130/0016-7606(1984)95<406:LCBE>2.0.CO;2
  • Kurikami et al., 2009. Study on strategy and methodology for repository concept development for the Japanese geological disposal project, NUMO-TR-09-04, Sept. 20-09, 101 pp.
  • Lay, T., and Kanamori, H., 1980, Earthquake doublets in the Solomon Islands: Physics of the Earth and Planetary Interiors, v. 21, p. 283-304.
  • Lay, T., Ammon, C.J., Kanamori, H., Kim, M.J., and Xue, L., 2011. Outer trench-slope faulting and the 2011 Mw 9.0 off the Pacific coast of Tohoku Earthquake in Earth Planets Space, v. 63, p. 713-718.
  • Lay, T., H. Kanamori, C. J. Ammon, A. R. Hutko, K. Furlong, and L. Rivera, 2009. The 2006 – 2007 Kuril Islands great earthquake sequence in J. Geophys. Res., 114, B11308, doi:10.1029/2008JB006280.
  • Liu, X., Zhao, D., and Li, DS., 2013. Seismic heterogeneity and anisotropy of the southern Kuril arc: insight into megathrust earthquakes in Geophysical Journal International, Volume 194, Issue 2, 1 August 2013, Pages 1069–1090, https://doi.org/10.1093/gji/ggt150
  • Meyer, B., Saltus, R., Chulliat, a., 2017. EMAG2: Earth Magnetic Anomaly Grid (2-arc-minute resolution) Version 3. National Centers for Environmental Information, NOAA. Model. doi:10.7289/V5H70CVX
  • Rhea, S., Tarr, A.C., Hayes, G., Villaseñor, A., Furlong, K.P., and Benz, H.M., 2010. Seismicity of the Earth 1900-2007, Kuril-Kamchatka arc and vicinity: U.S. Geological Survey Open-File Report 2010-1083-C, 1 map sheet, scale 1:5,000,000.
  • Van Horne, A., Sato, H., Ishiyama, T., 2017. Evolution of the Sea of Japan back-arc and some unsolved issues in Tectonophysics, v. 710-711, p. 6-20, http://dx.doi.org/10.1016/j.tecto.2016.08.020

Return to the Earthquake Reports page.


°
≥
ñ

2 thoughts on “Earthquake Report: Hokkaido, Japan”

  1. Pingback: Deadly earthquake below Hokkaido, Japan - September 5, 2018
  2. Pingback: Deadly earthquake below Hokkaido, Japan - September 5, 2018 | Best Science News

Leave a Reply

Your email address will not be published. Required fields are marked *