Earthquake sequence simulations with measured properties for JFAST core samples

1. Introduction

After the devastating 2011 Tohoku-Oki earthquake, many numerical models were published for the dynamic rupture (e.g. [1,2]) and the sequence of earthquakes (e.g. [3–7]) to explain the behaviour of the Japan Trench subduction fault, raising potential key mechanisms (e.g. a hierarchical asperity concept, a free surface effect, dynamic weakening, and its combined effect with rate strengthening at low slip rates). Meanwhile, new data have been accumulating in both observational and experimental fields, which are useful in reducing the model parameter space significantly. Here, we would like to present dynamic earthquake sequence simulations accounting for experimental outcomes and explaining observations. We believe that updating numerical models with new data is an important effort towards realistic modelling of natural fault behaviour.

It was known prior to the earthquake that the coastal region near Sendai, one of the main cities in Tohoku region, Japan, has been struck by devastating tsunamis at intervals of approximately 1000 years, as revealed by dating of tsunami deposits [8], with the predecessor event being the 869 Jogan Sanriku earthquake. Further investigation published after the Tohoku-Oki earthquake [9,10] identified more tsunami deposits and inferred a shorter recurrence interval between 500 and 800 years.

The Japan Trench Fast Drilling Project (JFAST) was conducted shortly after the earthquake. In this project, sea floor drilling to the coseismic slip surface was conducted shortly after the Tohoku-Oki earthquake about 6.4 km from the trench [11]. For detailed information, please see the previous publication [11]. Temperature measurements by the slip surface in this project revealed that frictional heat of 27 MJ m⁻² (with 90% confidence interval 19 to 51) was generated [12] for the coseismic slip of about 50 m [13].

On a larger scale, data of heat flow measurements in the Japan Trench subduction zone within about 250 km from the trench were explained by long-term heat generation on the subduction interface comparable to the effective friction coefficient (long-term heat generation area-density per slip normalized by overburden pressure) of about 0.025 [14].

Laboratory-measured permeability of the smectite-rich fault material picked up from drilled core samples in JFAST was quite low, below 10⁻²⁰ m² at the in situ condition (822 mbsf (metres below sea floor), ‘very likely from the plate boundary’) [15], indicating that frictional heat probably caused significant build-up of pore fluid pressure during coseismic slip (thermal pressurization, TP (e.g. [16])) there. The frictional properties of the sample collected from the same core were measured under hydrothermal conditions [17], showing that the rate dependences change with slip rate and temperature. In addition, frictional properties of blueschist, which presumably exists along the subduction interface in the seismogenic depth range [18], were measured under hydrothermal conditions [19]. The consequence of the experimental data should be investigated by quantitative methods such as numerical simulations. In doing so, we need to revise a currently available friction law (e.g. logarithmic rate- and state-dependent friction (RSF) law [20,21]) so that it can express such complex experimental results.

The present study aims to investigate whether it is possible to develop models of earthquake sequences on the Japan Trench subduction fault which account for the experimental outcomes and explain the observations raised above. Note that the existence of such a model is not obvious, since direct applicability of experimental results for centimetre-scale samples is not guaranteed because of potential scale effects. Section 2 is devoted to a description of the friction law used in the present study. The model setting is explained in §3, and the results and comparison with observations are presented in §4. Discussion and conclusions follow in §§5 and 6.

2. Friction law

In the shallow part of the fault, we adopted experimental results for the JFAST core sample [17], smectite-rich material collected from 822 mbsf [11]. We modelled the data by a slightly generalized rate- and state-dependent friction law. In addition, TP (e.g. [16]) was implemented there with measured hydraulic properties for samples from the same core section [15]. Since smectite dehydrates as it goes deep, we used the model for JFAST friction up to 120°C ambient temperature, and experimental results for blueschist [19] for the deeper part. The detailed distribution of the frictional properties is described in the next section.

(a) Quadratic rate- and state-dependent friction law for JFAST samples

In conventional RSF laws, the friction coefficient f logarithmically depends on slip rate V and on a state variable φ which represents recent slowness. For convenience, let us define

2.1

where V_* is an arbitrary reference slip rate. We have L_V= L_W at steady state for commonly used RSF laws. In addition, f is usually expressed as

2.2

where a and b represent the amount of direct and evolution effects, respectively, and f_* is the friction coefficient at the reference state L_V= L_W= 0. The friction coefficient at steady state f_ss is given by

2.3

For JFAST samples, experimentally determined a and a–b from velocity step tests depend on temperature and slip rate [17], as shown in figure 1. Specifically, with increasing slip rate, a − b changes from negative to positive at low temperature (approx. 20°C) but from positive to negative at high temperature (≥150°C). In addition, the velocity step tests from 30 to 100 µm s⁻¹ yield positive a − b values at 20 and 200°C, but typically negative values in a temperature range from 50 to 150°C. We aim to implement those dependences in the friction law. The rate-dependent a − b cannot be expressed by writing a and b as functions of the ambient temperature T_a, as is often done in previous studies (e.g. [22]). We read equation (2.2) as the separation of three effects: the absolute level, the direct effect, and the evolution effect. From this viewpoint, we generalize equation (2.2) to

2.4

where F and G represent the direct and evolution effects, respectively. They satisfy

2.5

In the present simulations, f refers to the ambient temperature T_a which does not change in time. The temperature in the fault core changes significantly during rapid slip which is much faster than the experimental range. At such high slip rates, other effects such as TP may dominate the change in f. Still, the experimental data at low slip rates can be used in constraining fault motion at low slip rates where the temperature of the fault core is not much different from T_a.

• Download figure
• Open in new tab
• Download PowerPoint

We have F = aL_V and G = −bL_W in conventional RSF laws. Here we assume that they are quadratic functions

2.6

and

2.7

where T_* is the reference ambient temperature. Then f_ss is written as

2.8

At least one quadratic term with respect to L_V is required to express the change in a − b with L_V. Since the slope changes with temperature (figure 1), a term with the coefficient (F₂₁+ G₂₁) is necessary. At least one quadratic term with respect to (T_a− T_*) is necessary to express the non-monotonic change in a − b as a function of temperature. The biquadratic term (F₂₂ + G₂₂)L_V²(T_a− T_*)² represents the change in the second derivative of a − b with respect to temperature with the change in the slip rate (the right column in figure 1b).

Because of the quadratic form, we can determine the parameters by fitting linear functions of L_V to the experimentally determined rate dependences [17] plotted against the average value of pre- and post-step L_V (figure 1). The best-fit parameters are listed in the electronic supplementary material, table S1. The plotted data points were equally weighted in determining the parameters. We selected the central values of the logarithmic velocities and the temperatures of the experimental condition as the reference values, and T_* = 104°C. The best-fit model successfully captured the tendencies recognized in the experimental study.

We can estimate f_* from the absolute value of f during the experiments [17] after determining the functions F and G. Although f_* may depend on T_a in principle, we cannot recognize a clear tendency (electronic supplementary material, figure S1). Therefore, we adopted the average value 0.369.

We used the aging law for the state-evolution law,

2.9

where t is time and L is state-evolution distance.

During earthquake sequences, V changes from orders of magnitude smaller than the plate convergence rate (approx. 10⁻⁹ m s⁻¹) to the coseismic slip rate (approx. 1 m s⁻¹). This range is wider than the range of experimental conditions. We did not simply extrapolate the quadratic functions, since we would obtain negative direct rate dependence (i.e. ∂f/∂V < 0), which could lead to an ill-posed problem [23]. We assumed that the derivatives of F and G were continuous and they were linear functions outside the experimental condition (i.e. V < 0.3 µm s⁻¹ or 100 µm s⁻¹ < V) (figure 2). At an ambient temperature below 20°C, we used T_a = 20°C to calculate f (figure 2a).

• Download figure
• Open in new tab
• Download PowerPoint

(c) Hydrothermal effect

We assume that TP operates at the shallow part of the fault [5–7]. We considered frictional heat and one-dimensional diffusion of temperature T and excess pore pressure p normal to the fault:

2.10

and

2.11

where y is the coordinate normal to the fault, α_th and α_hy are diffusion coefficients of temperature and pore pressure, respectively, ρc is specific heat capacity, Λ is pore pressure change per temperature change under undrained conditions and ω is heat generation power density. We assume a Gaussian heat source:

2.12

where w is the half-width of the heat source distribution and τ is shear resistance given by

2.13

with σ_ea the ambient effective normal stress and p₀= p|_y=0. We assumed w = 1 cm.

The solution technique for TP in the present study [24] is suitable for earthquake sequence simulations owing to unconditional stability, ability to take a logarithmically wide range of time steps, and affordable numerical cost. But it cannot deal with nonlinearity such as temporal changes in the transport properties. We used the measured hydraulic properties for the JFAST sample corrected from the slip surface [15] after interpolating to the in situ condition where the sample was collected. The parameters are listed in the electronic supplementary material, table S3.

At a greater depth where the rock type is different, the hydraulic properties should be different. However, there is no measurement available for materials from seismogenic depths. In addition, fracturing of indurated materials such as cataclasites presumably causes significant dilatancy associated with a decrease in p and an increase in α_hy. TP may become less efficient there. We express this effect by decreasing Λ to 0 linearly with increasing depth. Investigation of the effect of dilatancy for indurated fault materials is an important future study.

3. Model setting

We conducted two-dimensional earthquake sequence simulations with a spectral boundary integral equation method with inertial effects (e.g. [22]). Previous earthquakes affect the stress distribution by static stress change due to cumulative slip distribution. We assumed the shear wave speed c_s = 3000 m s⁻¹, Poisson's ratio ν = 0.25 and density ρ = 2700 kg m⁻³. The fault was discretized in the along-dip direction (x-direction) (figure 3a), and slip in this direction was considered. Spatial periodicity was assumed for efficient computation. The free surface was not considered, but an axis of rotational symmetry of order 2 was assumed at x = 0 as an approximation of the trench axis. The periodicity causes another axis at the deep end of the calculation domain x = 600 km. We prescribed a constant slip rate V_pl = 8.5 cm yr⁻¹ for x > 300 km, which loads the remaining part of the fault.

• Download figure
• Open in new tab
• Download PowerPoint

The dipping angle of the fault and the slope of the sea floor were assumed as 8.8° and 2.6°, respectively, which were used to estimate overburden stress on the fault. The distribution T_a was assumed based on thermal modelling in a previous study [14]. We assume that T_a is 0°C at x = 0, increases linearly with x, and reaches 400°C at a horizontal distance of 250 km from the trench (figure 3b).

We assumed the frictional and hydraulic properties of JFAST samples for T_a < 120°C, the frictional properties of blueschist and Λ = 0 for 150°C < T_a < 300°C, and a = F₁₀= 0.02, b = −G₁₀= 0 and f_*= 0.84 (f_ss at V_pl is the same as that for blueschist) for 350°C < T_a. Between these regimes, parameters f_*, F_ij, G_ij and Λ are linearly interpolated. The resulting apparent a − b value (i.e. ∂f_ss/∂ln(V)), which is a function of x and V, is plotted in figure 3c. At V = V_pl, the fault shows rate strengthening in about 60 km < x < 80 km. The shallowest part is rate weakening at V_pl, but ∂f_ss/∂ln(V) increases with V and becomes positive at about 10 µm s⁻¹.

Previous studies (e.g. [25]) suggested the existence of excess pore pressure at depth. We tested different distributions of the ambient effective normal stress σ_ea after parametrization. We considered two linear distributions and connected them smoothly:

3.1

A shorter state-evolution distance L causes a shorter recurrence of earthquakes [26]. Experimentally determined L is a fraction of 1 mm [17,19], which is too short for restricted numerical resources. We used a uniform distribution of L for simplicity, and changed its value from 10 mm to 40 mm to understand the dependence of the system behaviour on L.

The length scale of a process zone at a rupture front is given by (e.g. [27])

3.2

At the initial condition, τ was given by steady-state frictional resistance at V = V_pl. The initial value of φ was assumed to be the steady-state value at V = V_pl/2.

4. Model results and comparison with observations

(a) Recurrence intervals

As stated in §1, recent studies [9,10] suggested a recurrence interval of devastating tsunamis and tsunamigenic earthquakes between 500 and 800 years. In addition, there are sites of M7+ class earthquakes in the deeper part of the ruptured area of the 2011 Tohoku-Oki earthquake (e.g. [28]) at about 100 to 200 km from the trench. In the Miyagi-Oki area, for example, a sequence of large earthquakes occurred with an interval of about 37 years. A favourable model should show such a supercycle behaviour.

In the present simulations, we obtained many earthquakes defined by a threshold in V of 0.1 m s⁻¹. They can be classified into two categories: (i) catastrophic earthquakes, which reached x = 0 and produced more than 10 m of slip there, and (ii) the others, most of which produced significant slip within a region about 100 km < x < 200 km. We compare the earthquake in the first category to the tsunamigenic earthquakes and those in the second category to the deeper earthquakes. Note that in two-dimensional problems 10 times difference in the length scale of an earthquake corresponds to two orders difference in two-dimensional potency, while it corresponds to three orders difference in potency and to a difference in moment magnitude by 2 (e.g. Mw 9 and Mw 7).

(i) Effect of shallow effective normal stress

Figure 4 shows magnitude–time (M-T) diagrams and cumulative slip distributions for various σ_ea50, i.e. 7.5 MPa, 20 MPa and 30 MPa, with a fixed value of L = 20 mm. Note that the time scales of the M-T diagrams and the slip scales of the slip distributions are different for different cases. In addition to the catastrophic earthquakes defined in the previous section (red markers in figure 4a–c), smaller earthquakes (blue markers in figure 4a–c) took place typically in the deeper part of the seismogenic zone approximately from x = 100 km to 200 km.

• Download figure
• Open in new tab
• Download PowerPoint

With increasing σ_ea50, the interval of the catastrophic earthquakes increases, and so does their magnitude. They are sometimes nucleated at the bottom of the seismogenic zone (x ∼ 200 km), but sometimes at the point where the smaller earthquakes are arrested (x ∼ 100 km). The number of earthquakes between one catastrophic earthquake and another (a supercycle) increases with increasing σ_ea50. With σ_ea50 = 7.5 km, the earthquakes marked in blue in the M-T diagram sometimes produce significant coseismic slip at x = 0 (figure 4d) though it is smaller than 10 m.

The average recurrence intervals of all the earthquakes and of the catastrophic earthquakes are plotted in figure 5a. Error bars represent standard deviation. The fault behaviour is calculated at least until the ninth catastrophic earthquake, and the earthquakes during the first supercycle were excluded in the averaging. The interval of all earthquakes does not change significantly on average, with its variance increasing.

• Download figure
• Open in new tab
• Download PowerPoint

Figure 5b shows histograms of the intervals of all earthquakes for the cases plotted in figure 4. The increase in the variance in the earthquake interval from the case σ_ea50 = 7.5 MPa is primarily due to the appearance of very short intervals. Note that the afterslip is graphically represented in, for example, figure 4f as regions of sparse lines in x > 200 km. After a catastrophic earthquake, smaller earthquakes near the deeper limit of the seismogenic zone occurred frequently during an afterslip [29].

Since the favourable recurrence interval of the catastrophic earthquakes is between 500 and 800 years [9,10], we selected σ_ea50 = 25 MPa and looked into the effect of L as presented in the next subsection.

(ii) Effect of state-evolution distance L

Figure 6 is similar to figure 4 for the cases with σ_ea50 = 25 MPa and with different values of L. We used common time and slip scales here. With decreasing L, the intervals of catastrophic earthquakes and all earthquakes both decrease by a similar amount (figure 7a). This means that, although the nucleation frequency increases, catastrophic earthquakes have to wait for the shallow part of the fault to get ready.

• Download figure
• Open in new tab
• Download PowerPoint

• Download figure
• Open in new tab
• Download PowerPoint

With L = 40 mm, there are two significant ranges of concentrations in the histogram (figure 7b). One is below 30 years for earthquakes during afterslips of catastrophic earthquakes, and the other is above 200 years for earthquakes generated by loading due to deep fault slip at a plate rate. With decreasing L by a factor of 4, the latter concentration moves to shorter intervals by a factor slightly larger than 2. This is roughly consistent with the scaling of L^1/2 suggested by a previous study [26]. It should be pointed out that the smaller L caused the larger variation in the size of the earthquakes [23].

In order to have non-catastrophic earthquakes about every 40 years, which is the case in the Miyagi-Oki area, we have to use a smaller value of L. We did not conduct simulations with such a small L owing to restricted numerical resources, but the obtained dependences of the recurrence intervals on σ_ea50 and L are useful in further tuning of the model. It is conceivable that, keeping the value of σ_ea50 = 25 MPa and decreasing L to 2–3 mm, the interval of the earthquakes in the deep part of the seismogenic zone would become smaller than 50 years, while that of the catastrophic earthquakes would decrease by only tens of years.

(b) Long-term frictional heat

For all cases shown in the last section, long-term frictional heat per slip (dissipative stress) was calculated by dividing the distribution of cumulative dissipation (time integral of τV) by that of cumulative slip. The first supercycle was not counted. The dissipative stress increases in x < 200 km with increasing σ_ea50 (figure 8a), while it is hardly affected by L (figure 8b).

• Download figure
• Open in new tab
• Download PowerPoint

As introduced in §1, the data of heat flow measurements within 250 km from the trench axis were favourably explained by an effective friction coefficient of 0.025 in a previous study [14] with significant data scattering, overall being covered by a range of friction coefficient from 0 to 0.06. Black dashed lines in figure 8 indicate these estimates. All of the simulated cases are well within the range, and consistent with the observation.

(c) Slip and frictional heat during catastrophic earthquakes

Coseismic slip and frictional heat during catastrophic earthquakes at x = 6.4 km, which is comparable to the location of the JFAST drilling site [11], are plotted in figure 9. The first catastrophic earthquake in a simulated sequence is excluded, since it is affected significantly by the initial condition. The cases with σ_ea50 = 25 MPa are consistent with the 50 m of coseismic slip during the Tohoku-Oki earthquake (upper panels in figure 9), inferred from direct measurement of seafloor topography [13]. The constraint of coseismic frictional heat (the optimum value 27 MJ m⁻² and the 90% confidence interval 19–51 MJ m⁻² [12]) is satisfied by the cases with σ_ea50 ≥ 20 MPa (lower panels in figure 9). For more detailed description of a catastrophic earthquake, please see supplementary information (electronic supplementary material, figure S3).

• Download figure
• Open in new tab
• Download PowerPoint

5. Discussion

During the simulations, we obtained sub-seismic (i.e. maximum slip rate less than 0.1 m s⁻¹) events at the very shallow part of the system. Figure 10a represents a spatio-temporal distribution of V during such an event in a case with σ_ea50 = 25 MPa and L = 10 mm. This event started from around 11 km < x < 21 km, but downward propagation was not significant. Although the source dimension is about 30 km (= c_s × 10 s), the event duration is about 30 s (figure 10b). Figure 10c shows a trajectory of (V, τ) at x = 15 km. After a long-term increase in τ (arrow in figure 10c), state evolution is associated with acceleration from the plate rate. The acceleration was halted after the trajectory reached the plot of steady-state frictional resistance. The change in the rate dependence from negative to positive for the JFAST sample suppressed further acceleration of the fault motion. Stress drop due to this event was a fraction of 1 MPa (figure 10d). The focus of the present paper is on the comparison of models and observations, and we do not report the further detailed description of the system behaviour due to restricted page length.

• Download figure
• Open in new tab
• Download PowerPoint

Forward models of sequences of earthquakes have so many degrees of freedom if we consider non-uniform distributions of fault properties. Using conventional RSF laws, some parameter studies have been published (e.g. [7]) in order to understand how the fault behaviour depends on the frictional parameters. It is true that those studies are useful in understanding the underlying physics based on a simplified friction law, but it is difficult to obtain a convincing model that can be compared to a famous earthquake because of the many unconstrained degrees of freedom (e.g. distributions of aσ_e, bσ_e, L, f₀σ_e, α_hy, Λ and so on). In addition, real experimental results are not so simple as what is expected from the conventional RSF laws as shown in figure 1. We think it is important to constrain the parameter space using available information before conducting parameter studies.

Friction experimental outcomes for granite [30] have long been used as typical properties in numerical simulations (e.g. [22]). Those for gabbro [31], which may have more similar chemical composition to the seafloor rocks, were shown to yield much better simulated fault motion in explaining GPS observations in Northern Cascadia than those for granite [32]. Setting assumptions of frictional properties guided by frictional experiments is a reasonable strategy in exploring the huge parameter space.

In this paper, we demonstrated the consequence of the measured physical frictional properties [15,17,19]. It was shown that the model can be tuned to agree with different kinds of observations listed in §1, and the tuned model yielded sub-seismic slow events very close to the trench. The direct applicability of laboratory experiments is, however, under debate. The state-evolution distance L may be scaled by the shear zone thickness if strain is not localized [33]. In this case, the slip rate in the friction law may have to be scaled in the opposite way if the shear strain rate matters in the friction law. It is noteworthy in considering this problem that the sign of rate dependence changes with V for the JFAST sample. Rate weakening would prefer strain localization, while rate strengthening would prefer pervasive deformation. Then the internal structure and the shear zone thickness may change for deformation conditions. Experimental investigation of the scaling of the friction law is very important, and establishment of a physics-based, rather than phenomenological, friction law consistent with the experimental results is a crucially important future task. But, given the small amount of the core sample available, such basic investigation should be done for analogous material for which we can check the repeatability of experiments.

In the present study, we adopted a simplified geometry for the subduction fault. Implementation of realistic geometry may cause larger coseismic slip near the trench (e.g. [2]) and thus the optimum distribution of σ_e may change significantly. However, the way we tuned the system to fit the observation is probably applicable to such a system: decreasing L to decrease the interval of the earthquakes and increasing shallow effective normal stress to increase the interval of catastrophic earthquakes. How the shallowest part of the fault behaves in a tuned model accounting for the realistic geometry is an important future study.

6. Conclusion

The five years after the 2011 Tohoku-Oki earthquake have illuminated the behaviour of the Japan Trench subduction fault much more clearly than just after the earthquake (e.g. [9–14]). In addition, suggestions about the frictional properties of the fault were made in experimental studies [15,17,19]. The present study showed that it is possible to develop numerical models of earthquake sequences accounting for the experimental suggestions and explaining the fault behaviour by performing parameter studies for a few parameters.

We used experimental results for core samples recovered during the Japan Trench Fast Drilling Project (JFAST) [15,17] and blueschist [19]. For the JFAST sample, rate dependences change with slip rate, as well as temperature. In order to implement such properties, we generalized a rate- and state-dependent friction law. Thermal pressurization of pore fluid (e.g. [16]) was also implemented in the shallow part using the measured hydraulic properties for the core sample. We parametrized the distribution of ambient effective normal stress by σ_ea50 and conducted a parameter study for it. Experimentally determined state-evolution distance L was too short for us to conduct earthquake sequence simulations for restricted numerical resources. Thus we decreased L from 40 mm to 10 mm and investigated the effect of L.

Roughly speaking, there are two types of earthquakes. One is a catastrophic earthquake which reaches the trench axis and produces significant coseismic slip there, and the other is a smaller earthquake confined in the blueschist region. By controlling the two parameters, the average interval of the two types of earthquakes can be controlled independently. In the case with σ_ea50 = 25 MPa and L = 10 mm, we obtained catastrophic earthquakes every about 500 years and smaller deep earthquakes every tens of years. The smaller earthquakes took place more frequently than average during afterslip of the catastrophic earthquakes [29]. The preferred model is consistent with a direct temperature measurement during JFAST [12] and large-scale heat flow measurements [14].

Sub-seismic events sometimes occurred in our numerical simulations close to the trench. Acceleration of the fault motion was suppressed due to the complex rate dependence of JFAST samples. This indicates the importance of frictional properties of the shallow part of the subduction interface on the generation of slow earthquakes there. Whether lab experimentally determined frictional properties can be directly used in large-scale simulations is a very important problem given potential scale effects. But we would like to emphasize that the consequence of experimental outcomes in fault behaviour should be investigated with quantitative methods such as numerical simulations, and experimental results could be used as a guide in exploring the huge model parameter space.

Data accessibility

The experimental data used to constrain the model were either published in journals or are in a PhD thesis as listed in the references. The data are available on request to the author, M.S. The software for computation of earthquake sequences was given by Prof. N. Lapusta.

Authors' contributions

H.N. formulated the friction law, conducted simulations and took the lead in composition. M.S. provided experimental data, and contributed in discussions. B.S. gave the basic idea, arranged collaboration and contributed in discussions. All authors gave final approval for publication.

Competing interests

The authors have no competing interests.

Funding

H.N. is supported by MEXT KAKENHI Grant Number 26109007. M.S. is supported by a grant-in-aid for JSPS Fellows (24.7181). B.S. is supported by MEXT KAKENHI Grant Number 26109007, and by the Ministry of Education, Culture, Sports, Science and Technology (MEXT) of Japan, under its Earthquake and Volcano Hazards Observation and Research Program.