IPhO 2017 · Theory problem 2

Earthquake, Volcano and Tsunami

Download PDFIPhO 2017

Indonesia is the supermarket of natural hazards. Almost all of the natural hazards have occurred in Indonesia, such as volcano eruptions, earthquakes and tsunami.

A. Merapi Volcano Eruption

Merapi volcano in Yogyakarta is one of the most active volcano in Java. Pyroclastic flows are well-known eruption characteristics of the volcano. The pyroclastic flow is a hot mixture of gas and rock which travels away from a volcano. In October 26th, 2010, Merapi showed his explosive character by producing an ash plume that reached 12 km altitude (Fig. 1) and pyroclastic currents displacing more than 20,000 people around the volcano.

Figure 1

(Fig. 1: Pyroclastic cloud during Merapi eruption, Courtesy of Volcanological Office of Yogyakarta, BPPTKG)

Let us look at the causes of the largest eruption of Merapi in 2010. It is known by geophysicists that the influence of the external water into the magma plays an important role to the explosive behavior of volcanic eruptions (hydro magmatic eruptions). Let’s assume that we dealt with a volcano as a system that consists of mixture of magmatic particles and water. The volcano vents structures and atmosphere are being boundary of the system. The explosive eruption is considered to be happening in two stages, (1) an instantaneous magma-water interaction, and (2) a system expansion. In the first stage, a mass of magma (\(m_m\)) at an absolute temperature (\(T_m\)) mixes with a mass of external water (\(m_w\)) at an absolute temperature (\(T_w\)). The thermal equilibrium is reached almost instantaneously. This interaction can be perceived as a nearly-constant volume process. Latent heat of evaporation of water and latent heat of melting of magma can be neglected.

A.1 Find the equilibrium temperature at the first stage in terms of the masses and heat capacity per unit mass of water \(c_{Vw}\) and magma \(c_{Vm}\). 0.5 pt.
A.2 Determine its equilibrium pressure at the first stage by assuming that the mixture can be modeled as ideal gas. Assume that the volume per unit mole of the mixture is \(v_e\). 0.3 pt.

The system expansion (the second stage) can occur through several possibilities, one of which is thermal detonation. Although such process is quite complicated, we can empirically measure the relative velocity of the erupted mixture. The velocity of gas during the eruption depend on the pressure \(p\), the total mass \(m\) and the volume \(V\) of the mixture in the conduit of a volcano.

A.3 Express the velocity of gas during the eruption in terms of \(p\), \(m\), and \(V\) up to a proportional constant \(\kappa\). 0.5 pt.

The observed pressure is around the order of 100 MPa. This makes the eruption (relative) velocity can be as large as ballistic velocity.

B. The Yogyakarta Earthquake

The 2006 Yogyakarta earthquake of magnitude \(M_w=6.4\), which destroyed many buildings in the Bantul and Yogyakarta area, occurred at 05:54:00.00 local time or 22:54:00.00 UTC. The earthquake was caused by a sudden displacement of the Opak fault segment (see Fig. 2). The hypocenter was located 15 km below the surface.

The seismic wave that propagates on the earth crust can be recorded using seismometer. The diagram from seismometer is called seismogram (Fig. 2 and 3, Lower graph). These seismograms represent the vertical ground velocity as a function of time recorded by the seismic station at Gamping Station Yogyakarta (YOGI) (Fig. 2) and Denpasar, Bali (DNP) (Fig. 3). In general, seismic wave consists of three wave types: the longitudinal or primary (\(P\)-wave), the transversal or secondary (\(S\)-wave), and the surface wave. The \(P\)-wave and \(S\)-wave travel in the subsurface while the surface wave travels along the Earth surface. Seismic waves traveling through subsurfaces to the stations can be divided into those that propagate in a straight line, those that are reflected by a layer’s boundary, and those that are refracted into the next layer. The longitudinal wave or the primary wave has the highest velocity, while the surface wave has the lowest velocity, around 60% of the \(P\)-wave.

Figure 2

(Fig. 2: The maps location of YOGI)

Figure 3

(Fig. 3: The maps location of DNP (Denpasar))

The distance between the epicenter (the projection of hypocenter on the Earth surface) and the YOGI and DNP stations respectively are 22.5 km and 500 km. The depth of the Earth crust layer in Java, Indonesia, is 30 km. Beneath the Earth’s crust is the Earth’s mantle layer. Just like other wave phenomena, seismic wave also satisfies the Snell’s law. The seismic wave can also be reflected by the mantle layer. In this problem we assume that the earth curvature is neglected.

B.1 Fig. 2 shows the seismogram at the YOGI station. Use the data to find the velocity of the \(P\)-wave in the Earth crust. 0.5 pt.
B.2 Find the travel time of the direct \(P\) wave and reflected wave due to the Yogyakarta earthquake that arrived at the DNP station in Denpasar. 0.6 pt.

By assuming that the Earth is composed of only two layers: the crust and the mantle, the primary wave propagates in the crust and in the mantle with different constant velocities. The velocity in the mantle is faster than in the crust. Note that \(P\)-wave refracted into the mantle at the right angle (\(90^\circ\)) is being partly refracted back into the crust along its entire path of propagation along the crust-mantle boundary.

B.3 Find the velocity of the \(P\)-wave in the mantle. 1.2 pt.

For a more realistic Earth structure, the crust can be divided into a number of thin layers so that the velocity of seismic wave is a function of the depth \(z\) according to \(v(z)=v_0+az\) where \(a\) is a constant and the hypocenter is approximated on the surface. In this model, the wave ray is curving.

B.4 Let us define the ray parameter \(p=\sin\theta(z)/v(z)\), where \(\theta(z)\) is the angle between the ray and the normal. Suppose a seismic wave arrives to the station with ray parameter \(p\); express the distance to the hypocenter in terms of \(p\), \(v_0\), and \(a\). Assume that the hypocenter is very close to the ground surface. 1.4 pt.
B.5 Find the travel time \(T\) from hypocenter to any station, in form of integral over \(z\). 1.0 pt.

The earth consists of a stack of homogeneous layers with the velocity of each layer is \(v_i\) and the thickness of \(\delta z_i\).

B.6 From the result of the previous problem, approximate the travel time (\(T\)) from the hypocenter to DNP station by assuming that the crust consists of only three discretized layers, \((i=1,2,3)\), characterized by \(v_1=6.65\) km/sec, \(v_2=6.97\) km/sec, \(v_3=6.99\) km/sec, \(p=0.143\) sec/km, \(\delta z_1=6.0\) km, \(\delta z_2=9.0\) km, \(\delta z_3=15\) km. 1.0 pt.
Figure 4

Fig. 4: A simplified model of earth’s layers.

C. Java Tsunami

The 2006 Pangandaran earthquake and tsunami occurred on July 17 at 15:19:27 local time off the coast of west and central Java. During the earthquake where the epicenter fault is on the ocean floor, the fault may be displaced producing a remarkably large water wave called tsunami. In other words, a tsunami is a shallow-water wave which is initiated by a tiny amplitude, but with an extremely large wavelength. Consider a sudden fault causing a lifting of some ocean floor as shown in Fig. 5. Assume that the energy of the earthquake is transformed to the potential energy of this raised ocean water. For simple model we approximate that the raised water has a geometry of a box with its area of \(\lambda L/2\) (where \(L\gg\lambda\)) and height of \(h\).

Figure 5

Fig. 5: Illustration for the tsunami wave \(d\) is the depth of the ocean.

C.1 Find the potential energy stored in the raised ocean water due to the earthquake with respect to ocean surface. Assume that the density of sea water is \(\rho\). 0.5 pt.
C.2 Find the speed of tsunami wave up to dimensionless factor. 1.2 pt.
C.3 Using energy argument, determine the amplitude of the tsunami wave as a function of the depth, assuming that the depth varies slowly and also knowing that at a depth of \(d_0\) the amplitude is \(A_0\). 1.3 pt.

Official English text of the problem.

Answers

The seismogram readings used below are approximately \(5\,\mathrm{s}\) after the earthquake at YOGI and \(72\,\mathrm{s}\) at DNP. Let \(v_c\) and \(v_m\) denote the crust and mantle speeds.

  • A.1: \(\displaystyle T_e=\frac{m_m c_{Vm}T_m+m_w c_{Vw}T_w}{m_m c_{Vm}+m_w c_{Vw}}\).
  • A.2: \(\displaystyle p_e=\frac{RT_e}{v_e}\).
  • A.3: \(\displaystyle u=\kappa\sqrt{\frac{pV}{m}}\).
  • B.1: \(\displaystyle v_c\simeq5.4\,\mathrm{km\,s^{-1}}\).
  • B.2: \(\displaystyle t_{\rm direct}\simeq92.5\,\mathrm{s},\qquad t_{\rm reflected}\simeq92.8\,\mathrm{s}\), both approximately \(93\,\mathrm{s}\) at the precision of the graph.
  • B.3: \(\displaystyle v_m\simeq7.6\,\mathrm{km\,s^{-1}}\).
  • B.4: \(\displaystyle X=\frac{2}{ap}\sqrt{1-p^2v_0^2}\), assuming \(a>0\) and \(0<pv_0<1\).
  • B.5: \(\displaystyle T=2\int_0^{z_t}\frac{dz}{v(z)\sqrt{1-p^2v(z)^2}},\qquad z_t=\frac{p^{-1}-v_0}{a}\).
  • B.6: Direct discretization gives \(\displaystyle T\simeq1.84\times10^2\,\mathrm{s}\). This layer sum has a geometric inconsistency with the stated \(500\,\mathrm{km}\) offset, explained below.
  • C.1: \(\displaystyle U=\frac14\rho g\lambda Lh^2\).
  • C.2: \(\displaystyle c=C\sqrt{gd}\), where \(C\) is dimensionless; linear shallow-water theory gives \(C=1\).
  • C.3: \(\displaystyle A(d)=A_0\left(\frac{d_0}{d}\right)^{1/4}\).

Solution

A.1 — Equilibrium temperature

During the rapid mixing, take heat exchange with the surroundings and boundary work to be negligible. Neglecting the specified latent heats and treating the heat capacities as constant, conservation of internal energy gives

\[ m_m c_{Vm}(T_e-T_m)+m_w c_{Vw}(T_e-T_w)=0. \]

Thus

\[ \boxed{ T_e=\frac{m_m c_{Vm}T_m+m_w c_{Vw}T_w} {m_m c_{Vm}+m_w c_{Vw}} }. \]

This is the initial-temperature average weighted by the two heat capacities.

A.2 — Equilibrium pressure

For \(n\) moles of the ideal-gas mixture,

\[ p_eV=nRT_e. \]

Since its molar volume is \(v_e=V/n\),

\[ \boxed{ p_e=\frac{RT_e}{v_e} =\frac{R}{v_e} \frac{m_m c_{Vm}T_m+m_w c_{Vw}T_w} {m_m c_{Vm}+m_w c_{Vw}} }. \]

Here \(R\) is the universal gas constant.

A.3 — Eruption velocity

Assume, as specified, that the velocity depends only on \(p,m,V\). Write

\[ u=\kappa p^\alpha m^\beta V^\gamma. \]

Using

\[ [p]=ML^{-1}T^{-2},\qquad [m]=M,\qquad [V]=L^3, \]

dimensional consistency requires

\[ \alpha+\beta=0,\qquad -\alpha+3\gamma=1,\qquad -2\alpha=-1. \]

Therefore

\[ \alpha=\frac12,\qquad \beta=-\frac12,\qquad \gamma=\frac12, \]

and

\[ \boxed{u=\kappa\sqrt{\frac{pV}{m}}}. \]

Equivalently, \(u=\kappa\sqrt{p/\rho_{\rm mix}}\), where \(\rho_{\rm mix}=m/V\). Dimensional analysis does not determine \(\kappa\).

B.1 — \(P\)-wave speed in the crust

Let the source depth be \(z_s=15\,\mathrm{km}\). The horizontal distance to YOGI is \(22.5\,\mathrm{km}\), so the direct source-to-station distance is

\[ s_Y=\sqrt{22.5^2+15^2} =27.04\,\mathrm{km}. \]

The first visible departure from the baseline in Fig. 2 is approximately at 22:54:05 UTC. Relative to the earthquake time 22:54:00, this gives

\[ t_Y\simeq5\,\mathrm{s}. \]

Hence

\[ \boxed{ v_c=\frac{s_Y}{t_Y} \simeq\frac{27.04}{5} \simeq5.4\,\mathrm{km\,s^{-1}} }. \]

The numerical calculations below retain \(5.4083\,\mathrm{km\,s^{-1}}\) internally, although the graph does not justify that many significant figures.

B.2 — Direct and reflected travel times to DNP

Write the horizontal source-to-DNP distance as \(X=500\,\mathrm{km}\), and the crust thickness as \(H=30\,\mathrm{km}\).

For the direct wave,

\[ s_{\rm direct}=\sqrt{X^2+z_s^2}, \]

so

\[ \boxed{ t_{\rm direct} =\frac{\sqrt{500^2+15^2}}{v_c} \simeq92.5\,\mathrm{s} }. \]

For reflection at the horizontal crust–mantle boundary, reflect the source across the plane \(z=H\). Its image is at depth

\[ z_{\rm image}=2H-z_s=45\,\mathrm{km}. \]

Unfolding the reflected ray gives a straight line from this image to DNP. Consequently,

\[ s_{\rm reflected}=\sqrt{X^2+(2H-z_s)^2}, \]

and

\[ \boxed{ t_{\rm reflected} =\frac{\sqrt{500^2+45^2}}{v_c} \simeq92.8\,\mathrm{s} }. \]

Both travel times are about \(93\,\mathrm{s}\); their calculated difference is approximately \(0.33\,\mathrm{s}\).

B.3 — \(P\)-wave speed in the mantle

The first departure from the flat baseline in Fig. 3 is approximately at 22:55:12 UTC, giving

\[ t_D\simeq72\,\mathrm{s}. \]

This precedes the direct and reflected crustal arrivals calculated above. We identify it with the critically refracted wave that travels along the mantle boundary and radiates back into the crust.

Let \(\theta_c\) be the critical angle measured from the vertical. Snell’s law gives

\[ \frac{\sin\theta_c}{v_c}=\frac{1}{v_m}, \qquad \sin\theta_c=\frac{v_c}{v_m}. \]

The downward crustal leg spans depth \(H-z_s\), and the upward leg spans depth \(H\). Their combined vertical displacement in magnitude is

\[ K=(H-z_s)+H=2H-z_s=45\,\mathrm{km}. \]

Their combined length and horizontal projection are therefore

\[ s_{\rm crust}=\frac{K}{\cos\theta_c}, \qquad X_{\rm crust}=K\tan\theta_c. \]

The distance travelled along the mantle boundary is

\[ X_{\rm mantle}=X-K\tan\theta_c. \]

Thus

\[ t_D= \frac{K}{v_c\cos\theta_c} +\frac{X-K\tan\theta_c}{v_m}. \]

Using \(\sin\theta_c=v_c/v_m\), this simplifies to

\[ \boxed{ t_D=\frac{X}{v_m} +\frac{K}{v_c}\sqrt{1-\frac{v_c^2}{v_m^2}} }. \]

Substituting the graph reading,

\[ 72= \frac{500}{v_m} +\frac{45}{5.4083} \sqrt{1-\left(\frac{5.4083}{v_m}\right)^2}, \]

with distances in kilometres and speeds in kilometres per second. Solving gives

\[ v_m\simeq7.554\,\mathrm{km\,s^{-1}}, \]

or, at appropriate precision,

\[ \boxed{v_m\simeq7.6\,\mathrm{km\,s^{-1}}}. \]

As a check, \(\theta_c\simeq45.7^\circ\), and the mantle segment is approximately \(454\,\mathrm{km}\) long. The crustal legs take about \(11.9\,\mathrm{s}\), and the mantle segment about \(60.1\,\mathrm{s}\).

B.4 — Distance for the continuously varying velocity

Take \(z\) positive downward and \(x\) horizontal. Snell’s law across the thin horizontal layers implies that

\[ p=\frac{\sin\theta(z)}{v(z)} \]

is constant along a ray. Therefore

\[ \sin\theta(z)=pv(z), \qquad \cos\theta(z)=\sqrt{1-p^2v(z)^2}. \]

On the downward leg,

\[ \frac{dx}{dz} =\tan\theta(z) =\frac{pv(z)}{\sqrt{1-p^2v(z)^2}}. \]

For a returning ray in a medium with \(a>0\), the deepest point occurs where the ray becomes horizontal:

\[ pv(z_t)=1. \]

Since \(v(z)=v_0+az\),

\[ \boxed{z_t=\frac{p^{-1}-v_0}{a}}. \]

With both source and receiver approximated as lying on the surface, the downward and upward portions are symmetric. The horizontal distance is

\[ X=2\int_0^{z_t} \frac{p(v_0+az)} {\sqrt{1-p^2(v_0+az)^2}}\,dz. \]

Set \(u=p(v_0+az)\), so \(du=pa\,dz\). Then

\[ X=\frac{2}{ap}\int_{pv_0}^{1} \frac{u}{\sqrt{1-u^2}}\,du. \]

Since the antiderivative is \(-\sqrt{1-u^2}\),

\[ \boxed{ X=\frac{2}{ap}\sqrt{1-p^2v_0^2} =\frac{2}{a}\sqrt{\frac{1}{p^2}-v_0^2} }. \]

The returning-ray solution requires \(0<pv_0<1\).

B.5 — Travel-time integral

A ray segment of length \(ds\) takes time

\[ dt=\frac{ds}{v(z)}. \]

Since \(\lvert dz\rvert=ds\cos\theta(z)\),

\[ dt= \frac{\lvert dz\rvert} {v(z)\sqrt{1-p^2v(z)^2}}. \]

Including both the downward and upward legs gives

\[ \boxed{ T=2\int_0^{z_t} \frac{dz}{v(z)\sqrt{1-p^2v(z)^2}} }. \]

For the linear velocity profile, this is

\[ \boxed{ T= 2\int_0^{(p^{-1}-v_0)/a} \frac{dz} {(v_0+az)\sqrt{1-p^2(v_0+az)^2}} }. \]

The singularity at the turning point is integrable. Evaluating the integral also gives

\[ T=\frac{2}{a} \ln\!\left( \frac{1+\sqrt{1-p^2v_0^2}}{pv_0} \right). \]

An equivalent form, useful for examining the numerical approximation in B.6, follows by subtracting \(pX\):

\[ \begin{aligned} T-pX &= 2\int_0^{z_t} \left[ \frac{1}{v\sqrt{1-p^2v^2}} -\frac{p^2v}{\sqrt{1-p^2v^2}} \right]dz\\ &= 2\int_0^{z_t}\sqrt{\frac{1}{v(z)^2}-p^2}\,dz. \end{aligned} \]

Hence, for a geometrically consistent ray,

\[ \boxed{ T=pX+2\int_0^{z_t}\sqrt{\frac{1}{v(z)^2}-p^2}\,dz }. \]

B.6 — Three-layer approximation

Retaining the surface-source approximation from B.4–B.5, discretize the first travel-time integral:

\[ \boxed{ T\approx 2\sum_{i=1}^{3} \frac{\delta z_i} {v_i\sqrt{1-p^2v_i^2}} }. \]

The factor of two accounts for passage through each layer on the downward and upward legs.

Using the numbers exactly as printed gives:

Layer \(i\) \(pv_i\) Two-leg contribution \(\displaystyle \frac{2\delta z_i}{v_i\sqrt{1-p^2v_i^2}}\)
1 0.95095 \(5.833\,\mathrm{s}\)
2 0.99671 \(31.863\,\mathrm{s}\)
3 0.99957 \(146.366\,\mathrm{s}\)

Thus the direct layer quadrature gives

\[ \boxed{ T\approx184.06\,\mathrm{s} \simeq1.84\times10^2\,\mathrm{s} }. \]

There is a consistency limitation in interpreting this number as the travel time to DNP. The horizontal distance of the same discrete ray would be

\[ X_{\rm layers} = 2\sum_{i=1}^{3} \frac{\delta z_i\,pv_i} {\sqrt{1-p^2v_i^2}} \simeq1280.9\,\mathrm{km}, \]

which differs substantially from \(500\,\mathrm{km}\).

If instead one inserts the prescribed \(500\,\mathrm{km}\) offset into the second form of B.5 and approximates its integral, one obtains the formal estimate

\[ \begin{aligned} T &\approx p(500) +2\sum_{i=1}^{3}\delta z_i \sqrt{\frac{1}{v_i^2}-p^2}\\ &=71.5+0.8934\\ &\approx72.39\,\mathrm{s}. \end{aligned} \]

The two continuum formulas are equivalent, but these coarse approximations fail to preserve their required distance relation. Thus \(184\,\mathrm{s}\) is the direct travel-time layer sum; \(72.4\,\mathrm{s}\) is a fixed-offset estimate. They cannot both describe the same literal three-layer ray with the supplied parameters.

C.1 — Potential energy of the raised water

Let the horizontal area of the raised box be

\[ S=\frac{\lambda L}{2}. \]

Taking the undisturbed ocean surface as zero gravitational potential, a horizontal slab of thickness \(dz\) at height \(z\) has mass \(\rho S\,dz\). Consequently,

\[ U=\int_0^h \rho Sgz\,dz =\frac12\rho Sgh^2. \]

Therefore

\[ \boxed{U=\frac14\rho g\lambda Lh^2}. \]

Equivalently, the water above the original surface has mass \(\rho Sh\), with its centre of mass at height \(h/2\).

C.2 — Tsunami speed

Use the shallow-water, small-amplitude approximations

\[ \lambda\gg d,\qquad h\ll d. \]

The horizontal fluid velocity is then approximately uniform through the depth, and the vertical pressure distribution is approximately hydrostatic.

Let \(\eta(x,t)\) be the surface displacement and \(u(x,t)\) the horizontal fluid velocity. Conservation of water volume, linearized about depth \(d\), gives

\[ \frac{\partial\eta}{\partial t} +d\frac{\partial u}{\partial x}=0. \]

The horizontal pressure gradient is \(\rho g\,\partial\eta/\partial x\). Newton’s second law therefore gives

\[ \frac{\partial u}{\partial t} =-g\frac{\partial\eta}{\partial x}. \]

Differentiating the first equation with respect to time and substituting the second,

\[ \frac{\partial^2\eta}{\partial t^2} =gd\frac{\partial^2\eta}{\partial x^2}. \]

This is the wave equation, so

\[ \boxed{c=\sqrt{gd}}. \]

At the level requested in the question,

\[ \boxed{c=C\sqrt{gd}}, \]

with a dimensionless factor \(C\); the stated linear shallow-water assumptions give \(C=1\). The density cancels because it multiplies both inertia and the gravitational restoring force.

C.3 — Amplitude as the depth changes

Assume slow depth variation, negligible reflection and dissipation, and constant wavefront width. Since \(L\gg\lambda\), the wave can be treated locally as a plane wave.

For a sinusoidal surface displacement of amplitude \(A\),

\[ \eta=A\cos(kx-\omega t), \qquad \langle\eta^2\rangle=\frac{A^2}{2}. \]

The mean gravitational potential energy per unit horizontal area is

\[ \overline{\mathcal E}_{\rm pot} =\frac12\rho g\langle\eta^2\rangle =\frac14\rho gA^2. \]

For a progressive shallow-water wave, continuity gives

\[ u=\frac{c}{d}\eta. \]

Its mean kinetic energy per unit horizontal area is therefore

\[ \overline{\mathcal E}_{\rm kin} =\frac12\rho d\langle u^2\rangle =\frac{\rho c^2A^2}{4d} =\frac14\rho gA^2. \]

Thus the total mean energy density is

\[ \overline{\mathcal E}=\frac12\rho gA^2. \]

Shallow-water waves are nondispersive, so their group velocity equals \(c=\sqrt{gd}\). Conservation of energy flux per unit wavefront width gives

\[ \overline{\mathcal E}\,c=\text{constant}, \]

or

\[ A^2\sqrt d=\text{constant}. \]

Using the reference amplitude \(A_0\) at depth \(d_0\),

\[ A^2\sqrt d=A_0^2\sqrt{d_0}. \]

Therefore

\[ \boxed{ A(d)=A_0\left(\frac{d_0}{d}\right)^{1/4} }. \]

This is Green’s law. Its validity requires small amplitude relative to depth and sufficiently slow depth variation; it does not describe the eventual breaking of the wave.

The solution was checked against the official solution.