(2017) Stress accumulation in the Marmara Sea estimated through ground-motion

57 Pages·2017·2.78 MB·English

1 Stress accumulation in the Marmara Sea estimated through ground-
2 motion simulations from dynamic rupture scenarios
3
4 Hideo Aochi
5 Bureau de Recherches Géologiques et Minières, Orléans, France
6
7 John Douglas
8 University of Strathclyde, Department of Civil and Environmental Engineering, UK
9
10 Thomas Ulrich
11 Ludwig-Maximilians-Universität, München, Germany
12
13
14 Submitted to Journal of Geophysical Research (25th November 2016)
15 In revision, on 10th February 2017
16
17 Keypoints
18 • Earthquake ground motions are simulated in the Marmara Sea region based on
19 dynamic rupture scenarios.
20 • Ground motion prediction equations (GMPEs) allow the revision of the probabilities
21 of model parameters used in a logic tree framework.
22 • The revised probabilities thereby support the hypothesis of the low stress conditions
23 for this part of the North Anatolian fault.
24
1
25 ABSTRACT
26
27 We compare ground motions simulated from dynamic rupture scenarios, for the seismic gap
28 along the North Anatolian Fault under the Marmara Sea (Turkey), to estimates from empirical
29 ground motion prediction equations (GMPEs). Ground motions are simulated using a finite
30 difference method and a 3D model of the local crustal structure. They are analyzed at more
31 than a thousand locations in terms of horizontal peak ground velocity. Characteristics of
32 probable earthquake scenarios are strongly dependent on the hypothesized level of
33 accumulated stress, in terms of a normalized stress parameter T (Aochi and Ulrich, 2015).
34 With respect to the GMPEs, it is found that simulations for many scenarios systematically
35 overestimate the ground motions at all distances. Simulations for only some scenarios,
36 corresponding to moderate stress accumulation, match the estimates from the GMPEs. The
37 difference between the simulations and the GMPEs is used to quantify the relative
38 probabilities of each scenario and, therefore, to revise the probability of the stress field. A
39 magnitude Mw7+ operating at moderate prestress field (0.6 < T ≤ 0.7) is statistically more
40 probable, as previously assumed in the logic tree of probabilistic assessment of rupture
41 scenarios. This approach of revising the mechanical hypothesis by means of comparison to an
42 empirical statistical model (e.g. a GMPE) is useful not only for practical seismic hazard
43 assessments but also to understand crustal dynamics.
44 1. INTRODUCTION
45 The part of the North Anatolian Fault under the Marmara Sea (Turkey) is recognized
46 as a major seismic gap, where a large earthquake of magnitude larger than 7 is expected in the
47 coming years (Atakan et al., 2002; Erdik et al., 2004; Parsons 2004). Seismic risk mitigation
48 plans have been developed around this gap, in particular for the megacity of Istanbul, in terms
49 of, for example, early warning systems using dense observation networks and structural
50 retrofitting. Regardless of the high seismic hazard/risk, the description of the potential source
51 area remains imprecise, principally because: the faults are running under the Marmara Sea;
52 the largest earthquakes were not recorded instrumentally as they occurred in the 18th century
53 or earlier (e.g. Ergintav et al., 2014 and references therein); and the current seismicity rate is
54 relatively low. Many studies have been carried out to assess potential earthquake scenarios in
55 this area (e.g. Erdik et al., 2004; Pulido et al., 2004; Oglesby et al., 2008; Oglesby and Mai,
56 2012). Among them, Aochi and Ulrich (2015), hereafter referred as AU15, have discussed
2
57 statistically the occurrence probability of various earthquake scenarios, which are dynamically
58 simulated along different non-planar fault systems for various levels of hypothesized stress
59 accumulation and hypocentral locations. Although each model parameter is assigned a
60 probability, this process is based on expert judgement through a logic-tree approach. It is
61 found that some scenarios correspond to extreme conditions, such as a very fast rupture
62 velocity (super-shear rupture), which leads to high-amplitude ground motions (e.g. Akinci et
63 al., 2016). The purpose of this article is to quantify the ground-motion variations as a function
64 of the hypotheses, especially with respect to stress accumulation, through comparisons with
65 empirical ground motion prediction equations (GMPEs) (e.g. Douglas, 2003), which are
66 simple models used extensively in engineering seismology to estimate ground motions at a
67 site given the occurrence of an earthquake of a particular magnitude and location.
68 There have been several attempts to validate numerical simulations of earthquake
69 ground motions using GMPEs (e.x. SCEC broadband platform validation exercise; Goulet et
70 al., 2015), which had been previously derived from curve fitting to observations. Aochi and
71 Douglas (2006) compare ground motions that were dynamically simulated from simple
72 earthquake scenarios (double-couple sources) to GMPEs for various intensity measures: peak
73 ground acceleration (PGA), peak ground velocity (PGV), response spectral acceleration, Arias
74 intensity and duration. Although there were no complex propagation or site effects in the
75 simulations, the results implied a strong dispersion of the ground motion due only to source
76 effects. Baumann and Dalguer (2014) show more systematically the consistency of ground-
77 motion simulations from dynamic rupture modelling with GMPEs for frequencies up to 1 Hz,
78 in particular for intermediate distances (between 10 and 45 km). They also observed
79 supersaturation (a large variation both in positive and negative senses) very close (less than 1
80 km) to the source. Although most of the studies are for planar faults, ground motion from
81 dynamic rupture simulations based on non-planar fault systems can also be coherent in terms
82 of different metrics such as PGV, duration and response spectral ordinates (Aochi and
83 Yoshimi, 2016).
84 In the following sections, we first explain the procedure of our analysis, in which
85 model parameters are “inverted” from the difference between the ground motion simulated for
86 different dynamic rupture scenarios and those predicted by GMPEs. Then we introduce the
87 simulations both in terms of the earthquake source and the 3D regional structure, such that the
88 calculated ground motions are comparable with available observations for this region.
89 Finally, we discuss the implications of our analysis for this area and more generally for future
90 seismic hazard assessments.
3
91
92
93 2. METHOD AND MODEL
94 In this section we present the method followed in this article and the models used to make the
95 calculations.
96 2.1 METHOD
97 Figure 1 illustrates schematically the forward modelling of earthquake scenarios and
98 ground-motion simulations in a logic-tree formulation (top) and an inversion with respect to
99 the GMPE for calibrating the model-parameter probabilities (bottom), from steps (1) to (4).
100 Each forward modelling is a deterministic process based on fracture mechanics and elasto-
101 dynamics under given conditions. Various model parameters are considered together with a
102 weighting probability based on expert judgment, as in many recent probabilistic seismic
103 hazard assessments (e.g. Douglas et al., 2014). In the forward modelling of AU15, three
104 variable model parameters are considered and for each of which three branches are given:
105 namely 3 × 3 × 3 = 27 simulations are carried out, with probabilities between 0.5 % and 7.5
106 %. The model parameters are: (1) fault geometry, (2) stress accumulation level and (3)
107 hypocentral location. This study focuses only on the variation in stress accumulation level,
108 which strongly affects the simulated ground motions. Previously the three increasing levels of
109 stress accumulation were attributed decreasing probabilities, i.e. higher probability for a lower
110 stress accumulation. This assumption is reasonable in the sense that a characteristic
111 earthquake may be triggered once the regional accumulated stress is just sufficient, without
112 requiring a high and unstable stress. The assigned probabilities, however, remain uncertain.
113 We previously remarked that some simulated scenarios show severe or extreme features, such
114 as very large slip velocity and/or very fast rupture propagation (i.e. super-shear rupture), and,
115 thus, lead to extreme ground motions (AU15). Such scenarios are most commonly associated
116 to a high accumulated stress and are given a low probability. However, they appear quite
117 frequently so that, in the sum, they become statistically important. The motivation of this
118 study is to quantify the resultant ground motions so as to better constrain the probabilities of
119 different stress accumulation levels.
120 For convenience, we briefly explain the stress conditions and friction that are essential
121 in dynamic rupture simulations. AU15 introduced a stress parameter T indicating the
4
122 closeness of a given Mohr circle to Coulomb friction stick threshold in a Mohr-Coulomb
123 diagram (Figure 2). Parameter T is a global, non-dimensional parameter that describes the
124 externally loaded principal stresses, namely the initial level of stress loading, and the
125 frictional parameters, and is uniform everywhere in a simulation. The condition of T = 1
126 means that the Mohr circle tangentially meets the static friction threshold τ :
p
127 τ =σ +μ⋅σ (1)
p 0 s n
128 whereσ is the cohesive force, μ the static frictional coefficient, and σ the applied normal
0 s n
129 stress. Following simulations of the 1999 Izmit earthquake (Aochi and Madariaga, 2003),
130 these parameters are set as σ= 5 MPa and μ = 0.3. Mapping the regional stress tensor on
0 s
131 the complex fault leads to spatially heterogeneous shear and normal stress acting on the fault.
132 In the simulations, the external principal stress direction is set according to the tectonic
133 movement (Le Pichon et al., 2003), and its absolute level varies with depth. The condition of
134 T = 0 means, in contrast, that the Mohr circle tangentially meets the dynamic friction
τ
135 threshold :
r
136 τ =μ⋅σ (2)
r d n
137 where the dynamic frictional coefficient μ is set equal to 0.24, 20 % less than μ , which
d s
138 allows modeling the released fracture energy during the earthquake(e.g. McGarr, 1999) (thus,
139 always τ >τ ). T is given by:
p r
τ −τ
140 T ≡ 0 r for an optimal fault orientation (Φ = 36.7°; see Figure 2). (3)
τ −τ
p r
141 We assume that the seismogenic zone extends down to 12 km depth and that its shallow part
142 is less favorable to rupture propagation. As a consequence, we introduce a depth-dependent
143 critical displacement [D of the slip weakening law (Figure 2c)], identical to the depth
c
144 dependent D used by Aochi and Madariaga (2003) to simulate the 1999 Izmit earthquake.
c
145 For a positive stress drop during an earthquake, T must be between 0 and 1. If a planar fault is
146 ideally oriented with respect to the external principal stress, the condition T = 0.5 easily
147 allows spontaneous rupture propagation (this corresponds to S = 2 of Das and Aki, 1977);
148 however this is not the case for a non-planar fault system. Figure 3 shows the principal
149 stresses as a function of stress parameter. Although the rupture process is sensitive to the
150 parameter T (or S), the external principal stress does not change significantly. Note that S is
151 only defined for a planar fault, while stress parameter T is defined in a Mohr-Coulomb
152 diagram such that stress conditions can be calculated for faults of any orientation. It should
5
153 also be noted that S is strongly variable along the fault traces for a given T. In this study, we
154 vary T when discussing in detail the probability distribution.
155 The purpose of this study is to constrain the probable stress levels by analyzing
156 statistically the calculated ground motions of various scenarios by comparing them to
157 estimates from GMPEs (final step of Figure 1). GMPEs are generally obtained by regression
158 analysis on intensity measures, such as PGA and PGV, computed from observed strong-
159 motion data. GMPEs are usually simple functions of earthquake magnitude, faulting
160 mechanism, source-to-site distance and local site conditions. The equations provide a median
161 value as well as its associated standard deviation. Therefore, the deviation of the simulated
162 ground motions with respect to the GMPEs indicates a likelihood of each case, and the degree
163 of this deviation can be used to update the probabilities of model parameters of interest among
164 the bins.
165 Scherbaum et al. (2009) and Delavaud et al. (2012) propose and demonstrate the
166 model selection of GMPEs using an information-theoretic approach. For every sampling point
167 ( x(i=1,...,N) ), function g(x) evaluates the likelihood of observing x (ground-motion
i
168 estimate) for a given reference model (GMPE, here). Negative average sample log-likelihood
169 (LLH) is, therefore, defined as
1 N
170 LLH(g)=− log (g(x)), (4)
N 2 i
i=1
171 which is a good estimate of the relative distance between a model and a given reference. A
172 small LLH means a close match between the model and the observation. Then, for residuals
173 from different models g (x),k =1,...K , the likely weight of each model p is given as
k j
2−LLH(gj)
174 p = (5)
j K 2−LLH(gk)
k=1
175 This criterion provides a simple expression to rank the models statistically. Hence, this can be
176 regarded as a probabilistic inversion process. Equation (5) provides a relative probability
177 among all the bins, which is used in logic tree (Figure 1). The bins may be biased for some
178 reasons but the probability provided by Equation (5) is valid in the framework of a logic tree.
179 In this study, we calculate the ground motions in the 3D structure using a finite
180 difference method (Aochi and Madariaga, 2003; Aochi et al., 2013) according to different
181 finite source scenarios dynamically simulated. For the purpose mentioned above, we adopt
182 PGV as the primary intensity measure as we find that our tests on moderate earthquakes (M of
183 about 5) show a consistency among the simulation, the observations and the GMPEs (section
6
184 2.3) and also because PGV is sensible at lower frequencies (around 1 Hz) which we calculate
185 than PGA. Below we also demonstrate that other metrics lead to similar findings.
186
187 2.2 MODELS
188 Various earthquake scenarios along the North Anatolian Fault under the Marmara Sea were
189 dynamically simulated using the 3D boundary integral equation method (AU15). Each of the
190 27 simulated scenarios has a probability between 0.5 and 7.5 % and results in earthquakes
191 with moment magnitudes (M ) between 5.49 and 7.56. Each scenario is given a weight
w
192 corresponding to the selected model parameters for fault geometry, stress accumulation level
193 and hypocentral location in a logic-tree formulation as explained in the previous section. The
194 assigned probabilities remain, however, qualitative and uncertain, in particular, for the stress
195 accumulation level, which significantly controls the consequent ground motions. Figure 4
196 presents three examples of dynamic rupture simulations based on three different stress levels
197 (T = 0.66, 0.75 and 0.97) along the same fault geometry and assuming the same hypocenter
198 (AU15). After forcing rupture on an initial crack in the same way, the subsequent dynamic
199 rupture propagation is controlled by energy balance due to the static and dynamic stress field
200 and a frictional law. The external principal stress determined by parameter T is assumed
201 homogeneous along the tectonic rotational movement of the Marmara block, but the applied
202 normal/shear stress field along the fault trace is heterogeneous with respect to the strike
203 direction. Irregularity in fault geometry perturbs the stress field. If the assumed accumulated
204 stress (T) is high enough, the rupture propagation overcomes the stress heterogeneity to
205 become a large earthquake. This tendency is what we generally observe in similar numerical
206 simulations of dynamic rupture process along non-planar faults such as the 1999 Izmit
207 earthquake nearby, for example (Harris et al., 2002; Aochi and Madariaga, 2003). As a result,
208 the ruptured area, the amount of slip and the magnitude of an earthquake, hence, become
209 larger (M 7.04, 7.27 and 7.46) as T increases (Figure 4) and these changes should influence
w
210 the ground motions (AU15).
211 The ground-motion simulations are carried out in a 3D structural model using a finite
212 difference method, as in AU15. The finite difference is based on approximations that are
213 fourth order in space and second order in time for a staggered grid scheme (Aochi and
214 Madariaga, 2003; Aochi et al., 2013, and references therein). In the model domain of 200 km
215 (EW) by 120 km (NS) by 40 km (UD), the 3D structural model combines a 3D tomographic
216 model (Bayrakci et al., 2013), bathymetry from the General Bathymetric Chart of the Ocean
7
217 and a 1D layered model (Karabulut, personal communication; See Figure 5). The 3D model
218 describes the basin structure of the Marmara Sea but it assumes a flat surface. The sea layer is
219 given by V = 1.5 km/s, V = 0 km/s and ρ = 1000 kg/m3. The 3D model of Bayrakci et al.
p s
220 (2013) gives only the structure of V , so we define V and the material density ρ simply by the
p s
221 following equations: V =V / 3and ρ=1000×(−0.045V 2 +0.432V +1.711), where V is in
s p s s s
222 km/s and ρ is in kg/m3 (Ludwig et al., 1970: Horikawa et al., 2009). Figure 5 shows a 1D
223 profile extracted from the constructed 3D model. The basin structure is pronounced down to
224 about 12 km. The minimum velocity (V ) included is about 0.8 km/s (V beneath the sea
min s
225 floor). We set the grid size of the finite difference scheme as Δs = 160 m. The maximum
226 frequency f that is reliable numerically in the fourth order approximation is estimated as
max
227 f =V /(5⋅Δs) = 1.0 Hz. Our simulations do not take into account any site effects beyond
max min
228 those accounted for directly by the simulations. However, we can assume that at the relatively
229 low frequencies considered, site effects may have minimal impact on the ground motions. The
230 maximum velocity included in this model is V = 7.94 km/s (V ). Therefore, a time step of
max p
231 the finite difference should be Δt≤0.5⋅Δs/V =0.010 s to ensure stability; and, hence, we
max
232 fix Δt = 0.008 s. Each calculation is run for a duration of 120 s. The quality factor Q for
233 anelastic attenuation is introduced as a damping term in the finite difference formulation
234 (Graves, 1996) with central frequency about 1 Hz. There are other methods to take Q into
235 account in visco-elastic (e.g. Day and Bradley, 2001; Kristek and Moczo, 2003) or visco-
236 plastic (e.g. Roten et al., 2014) formulations. However, as shown in the next section, our
237 current knowledge of this region is still limited, even in terms of velocity structure, and a
238 detailed description of attenuation would not be useful yet.
239
240 2.3. CURRENT MODEL PERFORMANCE
241 We first simulate a moderate earthquake to check the reliability of the 3D model and
242 the accuracy of the ground-motion simulations. Figure 6 shows maps of the 2011/07/25 M
L
243 5.2 earthquake (40.823°N, 27.743°E, 10.8 km depth) and the 2012/06/07 M 5.1 earthquake
L
244 (40.853°N, 27.920°E, 7.1 km depth) (Karabulut, personal communication; who also provided
245 their focal mechanisms) and available seismograph stations (Broadband stations from Kandilli
246 Observatory of Earthquake Research Institute, Turkey). For the comparison, we calculate the
247 ground motions using the 1D layered model with a fixed Q = 300 and the 3D model assuming
248 Q = V/10 (Figure 5). Moment magnitude M of 4.85 and 4.9 are given to obtain the best fit to
s w
8
249 the spectral accelerations at 20 seconds. The source time function is assumed as a spline
250 function of durations of 2.0 s and 1.8 s, respectively. Figure 7 shows the waveforms filtered
251 up to 1 Hz at the common stations for the two earthquakes, aligned at the origin time (t = 0).
252 The arrivals of the waves are generally good except for SLVT where the observed wave
253 arrives slightly later than the simulated one. The peak values of the waveforms (here, PGV)
254 are mostly caused by the direct S-wave, but they are associated with the later phases at station
255 SLVT. In our models, no 3D structure exists beneath SLVT. This indicates that the model
256 could be further improved locally to reflect the geological context. As well as considering a
257 single ground-motion parameter, such as PGV, we also calculate the Goodness-of-Fit (GoF)
258 score between simulations and observations based on several metrics from each time series
259 (Olsen and Mayhew, 2010; Aochi and Yoshimi, 2016): averaged over long-period SA (10, 20
260 and 50 s), short-period SA (1, 2 and 5 s) and duration defined as interval between 5 and 95%
261 of Arias Intensity. GoF for two positive values (x, y) to compare is calculated with
2 x−y
262 GoF =erfc( ) (6)
x+ y
263 where erfc is the complementary error function. Figure 8 shows the scores averaged over the
264 defined metrics for each component for the two earthquakes. We do not aim to discuss the
265 scores of any particular station; we intend to show that the 3D model is better than the 1D-
266 layered model in this area. In general, GoF scores larger than 0.8 and 0.6 are considered
267 “excellent” and “good”, respectively, so that the current 3D model could be further improved.
268 Note that a factor-of-two difference between values leads to GoF = 0.35, although it would be
269 better to compare spectra not by its linear amplitude but by its logarithm. Hereafter we use
270 the 3D model assuming Q = Vs/10 for ground motion simulations of large scenario
271 earthquakes. Assumption of much stronger attenuation of Q = Vs/20 does not give a better
272 result in terms of GoF, as the duration becomes very short.
273 We also compare the horizontal PGVs with GMPEs in Figure 9. The selected GMPEs
274 are Bindi et al. (2014), Akkar et al. (2014a, b), Boore et al. (2014) and Cauzzi et al. (2015).
275 We show the non-filtered values of PGV from both simulations and observations because the
276 duration of the source time function of a smooth B-spline function is given for 1 s in the
277 simulations. The comparisons show a good coherency between the simulations and the
278 observations for both the moderate earthquakes. Although the used source parameter is simple
279 (a smooth source time function at a point source) in these tests, this agreement is supportive of
280 applying any of these GMPEs at a frequency of around 1 Hz for large scenario earthquakes
281 for which no observations exist. For the nearby 1999 Izmit earthquake, Bouchon et al. (2002)
9
282 showed the near-field ground velocities with no filtering and used them directly for
283 seismological analyses. The waveforms were dominated by low frequency content: strong
284 pulses of about a few seconds are visible and are closely related to the earthquake source
285 process (e.g. Aochi et al., 2011). In the following section, we use the GMPE of Bindi et al.
286 (2014), which is found roughly in the middle of the considered GMPEs and the simulations
287 and was found by Douglas and Aochi (2016) to be suitable for ground-motion simulations
288 from moderate earthquakes in this region. Similar results would be found for any of the other
289 three GMPEs as they also predict similar ground motions.
290 3. Results
291
292 Let us first explore the PGV maps simulated for an M 7.27 earthquake scenario (central
w
293 hypocenter position with T = 0.75, shown in Figure 4B). We calculate the ground motions
294 assuming a constant Q = 300 or a variable Q = V/10. Figure 10 shows the PGV maps and
s
295 PGVs as a function of distance, evaluated at 2851 distributed receivers. Both the calculations
296 tend to overestimate the ground motion levels predicted by the GMPE of Bindi et al. (2014).
297 For comparison, the three other recent GMPEs considered above are also shown. Although
298 some differences are found at distances shorter than 5 km, all these GMPEs show similar
299 overall features. The simulated PGVs are high particularly in the near field compared with the
300 predictions of GMPEs, but we find that strong attenuation (Q = V/10) shows a smaller
s
301 average residual (in terms of common logarithms): 0.402 compared to 0.467. Hereafter, our
302 discussion is based on the calculations assuming Q = V/10. Ground motions at very short
s
303 distances of a few numerical grids of FDM should not be taken into account because each
304 fault element of 500 m size within the BIEM simulation is represented as point source in the
305 ground-motion simulations, namely seismic moment release on a finite sub-element of finite
306 surface is concentrated at a point. The residuals are calculated for distances larger than 1 km.
307 For further distances, the receiver coverage is no longer azimuthally uniform as our physical
308 model dimension is 200 km × 120 km. Thus we use distances up to 60 km, which is roughly
309 the distance from the fault to the outer model boundary. We recall that the purpose of this
310 study is not to adjust the model parameters to fit the GMPEs but to assess the probabilities of
311 the model parameters using the relative discrepancy between the simulations and the GMPEs.
312 We keep the same fault geometry as in Figure 4 [Model LP, based on the fault trace
313 mapping of Le Pichon et al. (2003)] and the three hypocentral locations (West, Center and
10