# Unsupervised machine learning reveals slab hydration variations from deep earthquake distributions beneath the northwest Pacific

Gilbert Mao, Thomas P Ferrand, Jiaqi Li, Brian Zhu, Ziyi Xi, Min Chen

► To cite this version:

Gilbert Mao, Thomas P Ferrand, Jiaqi Li, Brian Zhu, Ziyi Xi, et al. Unsupervised machine learning reveals slab hydration variations from deep earthquake distributions beneath the northwest Pacific. Communications Earth & Environment, 2022, 3 (1), pp.56. 10.1038/s43247-022-00377-x. insu-03690222

HAL Id: insu-03690222

https://insu.hal.science/insu-03690222

Submitted on 8 Jun 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Unsupervised machine learning reveals slab hydration variations from deep earthquake distributions beneath the northwest Pacific

Gilbert L. Mao<fnref n="1" />,², Thomas P. Ferrand³,⁴✉, Jiaqi Li<fnref n="1" />,⁵, Brian Zhu<fnref n="1" />,⁶, Ziyi Xi<fnref n="1" /> & Min Chen<fnref n="1" />,⁷

Although transformational faulting in the rim of the metastable olivine wedge is hypothesized as a triggering mechanism of deep-focus earthquakes, there is no direct evidence of such rim. Variations of the *b* value - slope of the Gutenberg-Richter distribution - have been used to decipher triggering and rupture mechanisms of deep earthquakes. However, detection limits prevent full understanding of these mechanisms. Using the Japan Meteorological Agency catalog, we estimate *b* values of deep earthquakes in the northwestern Pacific Plate, clustered in four regions with unsupervised machine learning. The *b*-value analysis of Honshu and Izu deep seismicity reveals a kink at magnitude 3.7-3.8, where the *b* value abruptly changes from 1.4-1.7 to 0.6-0.7. The anomalously high *b* values for small earthquakes highlight enhanced transformational faulting, likely catalyzed by deep hydrous defects coinciding with the unstable rim of the metastable olivine wedge, the thickness of which we estimate at ~1 km.

<footnote n="1">¹Department of Computational Mathematics, Science and Engineering, Michigan State University, East Lansing, MI 48824, USA. ²William Mason High School, Mason, OH 45040, USA. ³Institute of Geological Sciences, Freie Universität Berlin, 12249 Berlin, Germany. ⁴Institut des Sciences de la Terre d'Orléans, UMR-7327, Université d'Orléans-CNRS, Orléans, France. ⁵Department of Earth, Planetary, and Space Sciences, University of California, Los Angeles, CA 90095, USA. ⁶Department of Statistics and Data Science, Yale University, New Haven, CT 06520, USA. ⁷Department of Earth and Environmental Sciences, Michigan State University, East Lansing, MI 48824, USA. ✉email: thomas.ferrand@fu-berlin.de</footnote>
Deep earthquakes, within subducting slabs at depths greater than 300 km, occur at conditions of high temperatures and pressures, which should prohibit brittle failure¹,². Their mechanism remains enigmatic due to the limited experimental observations, unresolved thermomechanical structural properties, and unclear source rupture processes. Three main hypotheses have been proposed in the literature to explain the occurrence of deep earthquakes: thermal shear instability³,⁴, dehydration embrittlement⁵,⁶, and transformational faulting⁷⁻¹⁰. While these three mechanisms are still up for debate, recent studies argue that the “dehydration embrittlement” model, based on fluid overpressure, is not a suitable mechanism for deep intraslab earthquakes, especially concerning the sinking lithospheric mantle¹¹⁻¹⁴. There has been an increasing number of seismic observations of metastable olivine wedges (MOW)¹⁵⁻¹⁸ and laboratory experiments⁹,¹⁹,²⁰ that support the triggering of deep earthquakes via transformational faulting from α-olivine to β-spinel (wadsleyite) in the rim of the MOW. The geometry of the MOW is controlled by the thermal parameter of the subducting slab²¹⁻²³. At the same depth, colder slabs such as Tonga are qualitatively inferred to have thicker MOW than warmer slabs, such as Japan and South America²⁴. Hydrous defects should also be considered as water activates metamorphic transformations and thus could prevent olivine metastability²⁵⁻²⁷. Consequently, besides temperature, the structure of the MOW is also influenced by the degree of slab hydration, which can account for the intermediate-depth seismicity distribution and is controlled by the structural inheritance and fault distribution of the oceanic plate prior to subduction¹²,²⁸. At intermediate-depths (30–300 km), the earthquake nucleation can originate from a dehydration-driven stress transfer (DDST), involving grain size reduction instead of dehydration-induced fluid overpressure¹¹,²⁷,²⁹,³⁰. However, it is unclear if slab dehydration and DDST can still be a viable explanation for deep-focus earthquakes.

The variability of the *b* value—slope of the Gutenberg-Richter distribution—provides important insight into the nature of seismic ruptures²⁴,³¹⁻³³. Especially for deep earthquakes, the *b*-value variations with slab thermal state and earthquake magnitude have been used to infer deep earthquake mechanisms: the nucleation can be due to transformational faulting within a metastable olivine wedge (MOW) and the rupture can propagate outside the MOW due to thermal runaway for larger earthquakes²⁴,³³. Although it is speculated that deep earthquakes should nucleate within the rim of the MOW¹⁵,¹⁶, there has been no seismic constraints on the dimensions of the rim, largely due to limited detection of small earthquakes. This prevents full understanding of both nucleation and dynamic propagation of deep ruptures.

Here we analyzed one of the most complete catalogs, the Japan Meteorological Agency (JMA) catalog, to estimate the *b* values of four clusters of deep earthquakes in the northwestern Pacific subduction zones, which host abundant deep earthquakes (Fig. 1; Fig. 2). We performed the earthquake clustering based on an unsupervised machine learning approach, i.e., the K-means algorithm, using the location of hypocenters. The *b*-value analysis of the deep seismicity in the Honshu and Izu segments of the Pacific slab reveals a kink at a magnitude of ~3.7–3.8 for earthquakes deeper than 300 km. The *b* value abruptly changes from 1.4–1.7 to 0.6–0.7 at the kink. This indicates a fractal dimension reduction for earthquakes above $M_w$ 3.7–3.8 in a distinct rupture domain within the MOW. The elevated *b* values greater than 1.5 are likely related to a highly hydrated rim of the MOW within the Pacific slab of the Honshu and Izu trenches. In comparison, no *b*-value kink is detected at the Kuril and Bonin trenches. Our unique observation of the *b*-value kink at $M_w$ 3.7–3.8 supports the existence of an unstable rim of the MOW with a characteristic

rim thickness of ~1 km (between ≈ 1 and ≈ 4 km), in which many low-magnitude earthquakes are triggered by transformational faulting. Additionally, earthquakes with intermediate magnitudes (>3.7–3.8) can rupture inward from the unstable rim into the (meta)stable interior of the MOW, with a characteristic thickness of 10–20 km. We propose that the anomalously high *b* values for small earthquakes ($M_w$ < 3.7–3.8) highlight enhanced transformational faulting catalyzed by deep hydrous defects. The coincidence between enhanced seismicity within the MOW unstable rim and a denser distribution of pre-subduction hydrated faults would lead to the observed frequency-magnitude distribution anomaly at lower magnitudes.

## Results

The *b* values calculated using three different methods are listed in Table 1. It has been shown that 2000 events can suffice to calculate the *b* value to within 0.05 at 98% confidence and 1000 events can still provide a robust *b* value to within 0.1 at 98% confidence³⁴. Another study also concluded that 1000 events yield small variation in the *b* value (roughly 0.1)³⁵. There are 2108, 5338, 1571, and 1643 events in the four defined regions (Table S7). Thus, we argue each cluster has enough events for a robust *b*-value analysis, as supported by the error analysis (Table 1). For the Kuril and Bonin clusters, the GR distributions show a linear trend with calculated *b* values of around 1.0 (Fig. 3a, d). However, clear kinks (breaks in the slopes of the GR distribution) are observed at a $M_w$ of 3.7–3.8, for the Honshu and Izu clusters (Fig. 3b, c), with *b* values decreased by half from 1.4 and 1.7 ($M_w$ < 3.7–3.8) to 0.6 and 0.7 ($M_w$ > 3.7–3.8), respectively. This drastic *b*-values variation across the kink is much larger than the standard deviations (<0.03).

To further validate the robustness of the above analysis, we resampled the earthquakes in each cluster using the Monte Carlo approach to conduct bootstrapping²⁴ (Supplementary Note 1; Fig. S1). The Monte Carlo simulations of the *b* values are plotted on histograms (Fig. S1). Although the *b* values from the Monte Carlo simulations show some difference from those calculated via the maximum likelihood estimation in the Gutenberg-Richter (GR) plot ($b_1$ in Table 1), the bimodal distribution of the *b* values in cluster 1 (Izu) and cluster 2 (Honshu) is obvious. Moreover, we confirmed that the kinks for the Izu and Honshu regions were not affected by the occurrence of specific large events affecting the regional stress field (e.g., Tohoku-Oki earthquake) (Supplementary Note 2; Table S1). These results further corroborate the robustness of the *b*-value analysis in each cluster and thus, the kinks observed in clusters 1 and 2 are significant and reflect structural characteristics.

Compared to the kink at $M_w$ 6.5, previously reported from the global CMT catalog²⁴, the kinks observed in this study are at much lower magnitudes (i.e., $M_w$ 3.7–3.8) for earthquakes deeper than 300 km. Because of the limited number of large deep earthquakes in each cluster from the JMA catalog covered by a short time-period between 2002 and 2016, the kink at $M_w$ 6.5 for earthquakes deeper than 500 km, is not observed. To confirm the possible existence of a kink at larger magnitude for earthquakes deeper than 300 km, we revisit the GR distribution of the CMT catalog with three more years of data compared to the previous study by Zhan²⁴. The deep earthquakes in “warm” subduction zones combined (Japan-Kuril (JK), Izu-Bonin-Mariana (IBM), South America (SA), and Philippine (PH)) show a clear kink at $M_w$ 6.7 (Fig. S5a). This kink occurs at a slightly larger $M_w$ likely due to a thicker MOW at 300 km than at 500 km. The calculated *b* values are almost doubled, from 0.6 below to 1.1 above $M_w$ 6.7, consistent with previous observations²⁴. The *b* value of 1.1 for the Tonga (cold) subduction without obvious kink also confirms previous results. The IRIS catalog with mixed magnitudes is also
<figure><img src="image_3.png" /><figcaption>**Fig. 1 Tectonic map and deep-focus seismicity distribution for the period 2002-2016 from the JMA catalog clustered by the K-means method.**</figcaption></figure>

Earthquakes in four clusters (clusters 0-3, i.e., result of the K-means algorithm) are colored corresponding to the Bonin (red), Izu (cyan), Honshu (purple), and Kuril (yellow) regions. Each cluster is correlated with different slab hydration state controlled by oceanic plate features and fabrics orientation and distribution. Major plate boundaries are outlined by the double thin black lines and microplate boundaries by double dashed lines. The key features of the oceanic floor prior to subduction are annotated accordingly with the legends for faults and fabrics. KFZ Kashima Fault Zone, NFZ Nosappu Fault Zone.

<figure><img src="image_4.png" /><figcaption>**Fig. 2 Depth distribution of the earthquakes of the JMA catalog for the 300-700 km depth range.** Mw > 2.4 (yellow); Mw > 3.4 (red); Mw > 4.4 (black).</figcaption></figure>
analyzed. The GR distribution shows a similar kink around $M_w$ 6.9, where the *b* value approximately doubles (Fig. S5b) for “warm” subduction zones, but remains constant (i.e., *b* value of 1.1) for the deep earthquake distribution in the colder Tonga subducting slab.

## Discussion

Constant *b* values have been viewed as an indication of self-similarity for earthquakes at all sizes. However, changes in *b* values (kinks) have been reported for both shallow and deep earthquakes<sup>24,36,37</sup>. The departure from a constant *b* value in a GR distribution has been attributed to either the saturation of the fault width<sup>37</sup> or a change of rupture mechanism at large magnitudes from transformational faulting to thermal shear instability<sup>24,33</sup>. Our observation of anomalously high *b* values (e.g., 1.4 and 1.7) below $M_w$ 3.7–3.8 in Honshu and Izu likely requires either high pore-pressure—unlikely at such depths—or strong heterogeneity that facilitates the occurrence of low-magnitude earthquakes<sup>31,38</sup>. Here the anomalously high *b* values and the existence of the kink at $M_w$ 3.7–3.8 happen to be

correlated to the highly hydrated lithospheric segments prior to subduction (Figs. 1, 2). For the regions of Honshu and Izu, the oblique subduction of spreading faults and oceanic fracture zones (e.g., Kashima Fault Zone) relative to the trench leads to the formation of numerous bending faults at the outer rise parallel to the trench. These bending faults enhance water percolation and deep serpentinization<sup>12,39–41</sup> resulting in highly hydrated slabs. In contrast, the slabs in Kuril and Bonin are relatively dry due to different reasons. In the Kuril region (Cluster 3), additional bending faults are not favored due to the pre-existing faults (i.e., spreading faults and Nozappu Fracture Zone) which are either already parallel or normal to the Kuril trench, which limits the opening of new bending faults and causes limited slab hydration. For the Bonin region, the dry slab is due to the warmer state of the Pacific lithosphere associated to the Cretaceous-Paleogene magmatism within the Ogasawara Plateau and its related ridges as results of the regional thermal anomaly (Fig. 1)<sup>42,43</sup>. We argue that the thermal anomaly associated with hotspot volcanism necessarily implies a warmer temperature profile in the lithosphere, and thus a shallower brittle-ductile transition as the bending faults penetrate only down to the brittle-ductile transition<sup>39</sup>. This implies that the Pacific Plate in Bonin is much drier prior to subduction compared to the three other lithospheric segments analyzed in this study. The dry MOW can be considered as relatively homogeneous, while its H-bearing transforming rim has strong intrinsic heterogeneity due to scattered transformation nuclei.

For the highly hydrated slabs in the regions of Honshu and Izu, the kinks observed at lower magnitudes ($M_w$ 3.7–3.8), with *b* values reduced by half from small to large earthquakes (Fig. 3b, c), suggest that the size of the $M_w$-3.7–3.8 events corresponds to the characteristic thickness of a seismogenic zone of less than 1 km<sup>44,45</sup>. This seismogenic zone at depths greater than 300 km is very likely the rim of the metastable olivine wedge (MOW), where olivine would transform to its high-pressure polymorph wadsleyite through transformational faulting, in a narrow unstable

<table><caption>**Table 1** *b* values and associated standard deviations using three different measurement methods.</caption><thead><tr><th>Region</th><th>b<sub>0</sub></th><th>σ<sub>b0</sub></th><th>b<sub>1</sub></th><th>σ<sub>b1</sub></th><th>b<sub>2</sub></th><th>σ<sub>b2</sub></th></tr></thead><tbody><tr><td>Kuril</td><td>0.97</td><td>0.01</td><td>0.87</td><td>0.01</td><td>0.92</td><td>0.01</td></tr><tr><td>Honshu<sup>a</sup></td><td>1.72</td><td>0.02</td><td>1.43</td><td>0.01</td><td>1.57</td><td>0.01</td></tr><tr><td>Honshu<sup>b</sup></td><td>0.65</td><td>0.02</td><td>0.60</td><td>0.02</td><td>0.62</td><td>0.02</td></tr><tr><td>Izu<sup>a</sup></td><td>2.11</td><td>0.01</td><td>1.70</td><td>0.01</td><td>1.89</td><td>0.01</td></tr><tr><td>Izu<sup>b</sup></td><td>0.80</td><td>0.03</td><td>0.73</td><td>0.02</td><td>0.77</td><td>0.02</td></tr><tr><td>Bonin</td><td>1.09</td><td>0.01</td><td>1.03</td><td>0.01</td><td>1.03</td><td>0.01</td></tr></tbody></table>

The *b*<sub>0</sub>, *b*<sub>1</sub>, and *b*<sub>2</sub> values are described in section Data and methods (b-value analysis).

The σ<sub>b</sub> are the standard deviations for each distribution.

<sup>a</sup>Mw < 3.8.

<sup>b</sup>Mw > 3.8.

<figure><img src="image_5.png" /><figcaption>**Fig. 3** Gutenberg-Richter distribution for earthquakes occurring between 2002 and 2016 deeper than 300 km in the four clusters. **a** Bonin (cluster 0), **b** Izu (cluster 1), **c** Honshu (cluster 2), and **d** Kuril (cluster 3). In **b** and **c**, very different slopes are observed below and above $M_w$ 3.7–3.8. The two straight lines *b* = 1.7 and *b* = 0.7 are the best fitting lines for 3.0 < $M_w$ < 3.8 (red) and 3.8 < $M_w$ < 5.5 (black) for cluster 1 (Izu), while *b* = 1.4 and *b* = 0.6 for 3.0 < $M_w$ < 3.7 (red) and 3.7 < $M_w$ < 5.8 (black) for cluster 2 (Honshu).</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Fig. 4** Schematic illustration of the hypothesis of transformational faulting enhanced by hydrous defects in the unstable rim of a metastable olivine wedge (MOW). **a** For a highly hydrated slab mantle, hydrous defects would coincide with the MOW transforming rim and make it a highly reactive halo that enhances the nucleation of small earthquakes, which leads to high *b* values (e.g., *b* > 1.5) for the Honshu and Izu deep-focus seismicity. **b** For a relatively dry slab mantle, there is limited seismicity inside the transforming rim and the *b* values are normal (e.g., *b* ≈ 1) for the Bonin and Kuril slabs.</figcaption></figure>

zone composed of a mixture of different olivine phases characterized by a small grain size<sup>1,19</sup>, which is key for seismic rupture nucleation and propagation<sup>29,46–48</sup>. While the MOW should not be expected in wet slabs due to the water-enhanced kinetics of the olivine-wadsleyite transition<sup>25,26,49</sup>, experimental evidence shows that the MOW can still coexist with hydrous phases (e.g., Phase A) in water-undersaturated conditions, which allows transformational faulting within dry MOW<sup>50</sup>. Considering water-undersaturated conditions at greater depth for highly hydrated slabs<sup>51</sup>, it is possible that hydrous defects in the rim will locally trigger fast transformation and generate stress transfers toward drier rock volumes inside the MOW. We want to emphasize again that the mechanism proposed here is different from the “dehydration embrittlement” model often mentioned in the literature. At depths greater than 300 km as in this study, we only consider the remaining water defects within the nominally anhydrous minerals (e.g., olivine) and/or at grain boundaries, i.e., long after the dehydration of serpentine minerals within the subducting lithospheric mantle.

We propose that the highly hydrated slabs in Honshu and Izu are likely to preserve substantial hydrous defects at greater depths spatially coinciding with the MOW rim closer to the slab top (Fig. 4a). These hydrous defects can be nominal without producing fluid overpressure<sup>27,30</sup> but should catalyze the transformational faulting by either minor water precipitation at grain boundaries<sup>52</sup> or via transformation-driven stress transfers<sup>11,27</sup>. When $M_w$< 3.7–3.8, the intrinsic heterogeneity introduced by hydrous defects in the rim can account for the observed high *b* values (1.4 and 1.7); in contrast, in dry slabs only negligible hydrous defects should spatially coexist with the MOW rim (Fig. 4b). This leads to limited low-magnitude seismicity nucleated in the MOW rim and normal *b* values close to 1.0 in the Bonin and Kuril slabs. When $M_w$ > 3.7–3.8, transformation-driven stress transfers (TDST)<sup>11,27</sup> can trigger rupture propagation beyond the MOW rim and into the metastable olivine wedge (fault lengths > 1 km Figs. 3, 4). We want to point out that small events can well occur inside of the MOW. Only the abnormally numerous small events observed in some slabs (like the Izu and

Honshu regions) in addition to the normal trend would be located within the rim.

The break in the GR distribution observed at $M_w$ 6.7 from the global CMT catalog (Fig. S5) for deep-focus earthquakes suggests a fault size of 10–20 km<sup>44,45</sup>. The increase of *b* values at $M_w$ 6.5 for earthquakes at depths > 500 km was previously attributed to a change in the fractal dimension of earthquake size distribution from 1 to 2 due to a thin MOW (with a width smaller than 3–5 km) and the change in controlling rupture mechanism<sup>24</sup>. We argue here that changes in *b* values at different sizes should be determined by changes of both the rupture domain heterogeneity and the physical mechanisms controlling the rupture, which operate at different stages from rupture nucleation to dynamic propagation and arrest<sup>53</sup>. For large earthquakes ($M_w$ > 6.7) at depths > 300 km, it is possible for the earthquake to rupture through the entire width of the MOW and potentially propagate outside the MOW (but still within the slab) into the neighboring wadsleyite through thermal runaway instability<sup>46,54</sup> (Fig. 5d). Assuming circular ruptures, the fault size of a $M_w$ 6.7 event would be around 10–20 km<sup>44,45</sup>. In such scenario, the change of controlling mechanism and the increase of heterogeneity (and thus fractal dimension outside the MOW) can explain the increase of *b* values to 1–1.5 for $M_w$ > 6.7 (Fig. S5). Furthermore, enhanced transformational faulting in the water-bearing MOW rim would induce a heterogeneity level high enough to increase the *b* values for small events ($M_w$< 3.7–3.8). Depending on the slab hydration state, the density of hydrous defects will modulate the *b* value at low magnitudes from 1 (dry slab with limited hydrous defects) to 2 (initially highly hydrated slab with abundant remaining hydrous defects; Fig. 5b).

The development of ruptures of different sizes can all be related to three distinct rupture domains: the unstable rim of the MOW, the interior of the MOW, and the warmer wadsleyite domain around the MOW including slab material that already endured the olivine-wadsleyite transformation. These three distinct materials are the locations of three different rupture mechanisms, which can operate at three different stages for large ruptures (Fig. 5a, c, d). Considering a large event, the initial stage involves
<figure><img src="image_7.png" /><figcaption>Fig. 5 Mechanisms proposed to explain the *b*-value kinks observed around $M_w$ 3.8 from the JMA catalog and $M_w$ 6.7 from the CMT catalog.</figcaption></figure>

a Schematic representation of the unstable rim (with an example thickness of ≈ 1 km) around the MOW (with an example thickness of ≈ 20 km), with deep earthquakes of different sizes: $M_w$ < 3.8 within the transforming rim, 3.8 < $M_w$ < 6.7 in case of stress transfer rupturing through the MOW, and $M_w$ > 6.7 for ruptures propagating beyond the MOW toward the stable slab material already transformed into wadsleyite peridotites. b Range of *b* values corresponding to the different segments of the Gutenberg-Richter distribution defined by the observed two kinks at $M_w$ 3.8 and $M_w$ 6.7. The slope of each segment is controlled by the hydration state of the slab and the thermal state of the MOW. c Illustration of the link between rupture size and ruptured materials for deep-focus earthquake: within the MOW rim (red), inside the MOW (green), and outside the MOW (blue). d Corresponding mechanisms required in each material: primary transformational faulting in the unstable rim of the MOW (red), transformation-driven stress transfer inducing secondary transformational faulting rupturing into the MOW (green), and the additional rupture propagation outside the MOW via thermal runaway instability (blue), with increasing rupture width and earthquake magnitude.

low-magnitude earthquakes nucleating within the MOW rim via transformational faulting when (meta)stable olivine becomes unstable and transforms into wadsleyite. As the fault grows larger, at the second stage, the rupture can propagate into the interior of the dry MOW via a transformation-driven stress transfer. If the rupture passes the center of the MOW (i.e., half width), the shearing instability likely continues via thermal runaway in the slab-strike direction within the MOW and slightly beyond the MOW further into the wadsleyite within the slab (Fig. 5c, d). The variations of the *b* values with all earthquake sizes can be explained by differences in fractal dimension between these three distinct rupture domains controlled by the slab hydrous state and the thicknesses of the MOW and its unstable rim, determined at first order by the slab thermal parameter. It should be noted that the *b*-value changes across different subduction zones for deep-focus earthquakes may indicate changes in physical rupture mechanisms. For example, for the cold Tonga subduction zone, the *b* value shows no significant change at larger magnitudes. This would suggest that the MOW in the Tonga slab is thick enough to accommodate large earthquakes with magnitude 7 or above with the same mechanism of transformation-driven stress transfer. However, due to the earthquake detection limit, the kink at smaller magnitudes still needs to be verified in the future with more complete earthquake catalogs for other subduction zones such as Tonga.

Furthermore, the observation of the *b*-value kink in the subducting slab can provide additional insights regarding features that characterized the oceanic lithosphere before subduction. For instance, we have enough data in the Bonin region in the 300–550 km depth range to report the absence of kink in the *b* value. As shown in the recent paper of Hirano et al.⁴², the hotspot volcanism that overprinted the Pacific lithosphere in the region of

the Ogasawara Plateau at ≈ 55 Ma also affected the Minamitorishima area at ≈ 35 Ma. Considering the relatively steep slab geometry in the Bonin subduction zone⁵⁵ and a subduction velocity of ≈ 80 mm.yr⁻¹⁵⁶, and assuming a fixed location of the hotspot volcanism⁴², a depth of 550 km would correspond to a subduction period of ≈ 15 Ma. This would suggest that the activity of the Ogasawara hotspot magmatism started at least around 80 million years ago, consistent with the literature⁴².

## Conclusions

We have analyzed the Gutenberg-Richter distribution for deep earthquakes from one of the most complete regional earthquake catalogs: the JMA catalog of 15 years of data, in the northwestern Pacific subduction zones. An unsupervised machine learning approach, the K-means algorithm, divides all the deep-focus earthquakes (depths ≥ 300 km) into four clusters corresponding to subducting slabs with different hydration state prior to subduction. For the first time, we observe kinks of *b* values reducing from 1.4–1.7 down to 0.6–0.7 at $M_w$ 3.7–3.8 for the Honshu and Izu clusters, while normal *b* values (1.0) for the Bonin and Kuril clusters are observed at low and intermediate magnitudes 3.3 < $M_w$ < 6.0. The variation of *b* values at low magnitudes ($M_w$ < 3.7–3.8) indicates their dependence on the slab hydration state, with higher *b* values (1.4–1.7) correlated with highly hydrated slab coinciding with the expected location of the MOW rim. The hydrous defects favor the nucleation of small earthquakes via the mechanism of transformational faulting within the rim. Such mechanism operates for small events ($M_w$ < 3.7–3.8), with rupture lengths ≤ 1 km, which tells us the thickness of the MOW rim. Through transformation-driven stress transfer, the rupture can propagate into the MOW, which is more homogeneous with
smaller fractal dimension, corresponding to b values between 0.5 and 1.0.

Combining with the b-value analysis from the latest CMT catalog, the kink at $M_w$ 6.7 suggests that the mechanism of thermal runaway instability operates for larger earthquakes rupturing through and propagating outside the MOW with increased heterogeneity in the new rupture domain. The change of controlling mechanism and rupture domain heterogeneity can explain the spatially varying b values due to the slab hydrous state and thermal state. The $M_w$ threshold of 3.7–3.8 provides a seismic constraint on the MOW rim thickness ~ 1 km, a narrow zone where hypocenters of low-magnitude earthquake should concentrate in, while the centroids of larger earthquakes below $M_w$ 6.7 should coincide with the geometry of the MOW with a thickness of 10–20 km. Furthermore, considering cold subduction zones such as Tonga and moderate hydration, it is possible that a more complete earthquake catalog reveals a similar kink and high b values at low moment magnitudes in case of the existence of a thin unstable rim.

## Data and methods

The previous b-value analysis of deep earthquakes involves using the global CMT catalog with moment magnitude of completeness ($M_C$) around 5.3$^{24}$. The high $M_C$ prevents the full understanding of earthquake nucleation at lower magnitudes. The depth distribution of earthquakes recorded in the JMA catalog is shown in Fig. 2. To study the distribution of low-magnitude earthquakes and its implications on deep rupture mechanism, we performed clustering and b-value analyses of deep-focus earthquakes using the JMA catalog, i.e., the most complete catalog of the Kuril, Japan, and Izu-Bonin subduction zones (Fig. 1) with a much lower $M_C$ around 3.0 (Fig. 3). There are a total of 10,660 earthquakes recorded at depths greater than 300 km in Northwestern Pacific subduction zones between 2002 and 2016. The reason to select the 300–700 km depth range for this study is detailed in the Supplementary Note 3 (Tables S2 to S6). To cluster these deep earthquakes more objectively, we employed the K-means clustering, an unsupervised machine learning method based on clustering analysis. The b-value analysis was then performed in each of the four clusters, which turned out to also correspond to different segments of the Pacific Plate with different slab hydration degrees (Fig. 1).

**Clustering analysis.** Cluster analysis or clustering is the process of finding similarities in the dataset and grouping similar data points together. It is one of the important tools for the pattern recognition$^{57}$ and generally regarded as unsupervised learning$^{58}$. Because the K-means algorithm is easy to implement and does not require high computational power, it is the most widely used algorithm for cluster analysis$^{59–61}$. Cluster analyses of earthquakes provide insights into seismic hazard in different geological or tectonic settings$^{62–64}$. In this study, we used the earthquake hypocenter locations as the features for the K-means clustering, with the intention of minimizing human interference. The optimal number of clusters is selected through silhouette analysis (Supplementary Note 4; Fig. S2). It evaluates how well the resulting clusters are separated, with higher silhouette scores (closer to 1) indicating better clustering. In this study, the optimal cluster number is four because it leads to not only a local maximum silhouette score (Fig. S2) but also clusters in excellent agreement with the existing segments of the subducting Northwestern Pacific Slab, i.e., Kuril, Honshu, Izu, and Bonin (Fig. 1). The robustness of the K-means clustering algorithm has been further corroborated by the spectral clustering and Gaussian mixture models (GMM) clustering algorithms. The clustering analysis results are similar for all three different algorithms and

included in the SI (Supplementary Notes 5-7; Figs. S3 and S4. Tables S8 to S11).

It is worth noting that the Bonin segment is characterized by a remarkable magmatic anomaly that overprints the oceanic lithosphere in the region during the Cretaceous and the Paleogene$^{42,43}$. The Ogasawara Plateau corresponds to the area affected by this magmatic activity, which necessarily has a prolongation within the lithosphere subducting at the Bonin trench. As the magmatic activity mobilizes water, we argue that it necessarily left the lithospheric mantle in a drier state than without additional magmatic activity (portion without overprinting subducting at the Izu trench). The lower water content in such a drier lithosphere (Bonin) is expected to reduce the probability of nucleation of metamorphic transformations$^{13,25–27}$; therefore, we argue that the separation of the Izu-Bonin Arc seismic events into two individual clusters, as done by the unsupervised machine learning in this study, is appropriate.

**b-value analysis.** To carry out the Gutenberg-Richter (GR) distribution analysis for each cluster, we first converted the JMA magnitude scale $M_j$ to the physics-based moment magnitude $M_w$ scale using the following empirical relationship$^{65}$:

$$M_w = aM_j^2 + bM_j + c \quad (1)$$

where $a = 0.053 \pm 0.003$, $b = 0.33 \pm 0.02$, and $c = 1.68 \pm 0.03$. The magnitude of completeness ($M_c$) is approximately estimated from the non-cumulative (binned) distribution (Supplementary Note 8; Table S12). Since the $M_c$ determined from all approaches are similar, the completeness determination in this study is reliable. Due to the seismic network coverage variations, $M_c$ changes across the different clusters, e.g., $M_c$ 3.3 for the Bonin cluster with fewer stations in its vicinity, while around 3.0 for the rest of the clusters with more station coverage (Fig. 3).

We then calculated b values ($b_0$) by means of the maximum likelihood estimation$^{66}$. We have also estimated the more accurate b values ($b_1$) following the modified formula$^{67–69}$. We note that Marzocchi and Sandri stated that either the original estimation ($b_0$) or the modified estimation ($b_1$) will be influenced by the possible error on the magnitude$^{35}$. Therefore, we also applied another b-value estimation ($b_2$) from Tinti and Mulargia$^{70}$, which is much less affected by the magnitude error$^{35}$:

$$b_0 = \frac{1}{\ln(10) * (\mu - M_c)} \quad (2)$$

$$b_1 = \frac{1}{\ln(10) * (\mu - (M_c - \Delta M/2))} \quad (3)$$

$$b_2 = \frac{1}{\ln(10) * \Delta M} * \ln \left( 1 + \frac{\Delta M}{\mu - M_c} \right) \quad (4)$$

where $\mu$ is the average of the magnitudes, $M_c$ is the magnitude of completeness, and $\Delta M$ is the magnitude interval (0.05 in this study).

We calculated the corresponding standard deviations of all three b-value estimation methods:

$$\sigma_{b_0} = 2.3 b_0^2 \sqrt{\frac{\sum_{i=1}^{N} (M_i - \mu)^2}{N(N-1)}} \quad (5)$$

$$\sigma_{b_1} = 2.3 b_1^2 \sqrt{\frac{\sum_{i=1}^{N} (M_i - \mu)^2}{N(N-1)}} \quad (6)$$
$$ \sigma_{b2} = \frac{\frac{\Delta M}{\mu - M_c}}{\ln(10) * \Delta M * \sqrt{N * (1 + \frac{\Delta M}{\mu - M_c})}} \quad (7) $$

where N is the total number of earthquakes.

### Data availability

Seismic datasets used for the clustering analyses in this study can be accessed from: ftp://ftp.eri.utokyo.ac.jp/pub/data/jma/hypo/ for the Japan Meteorological Agency (JMA) catalog; https://ds.iris.edu/wilber3/find_stations/ for the IRIS catalog; and http://www.globalcmt.org/ for the CMT catalog.

### Code availability

The pre-existing Python codes from the *Scikit Learn* cluster library were used for the clustering analyses. The map-view figure is produced using the Generic Mapping Tools (GMT). The scripts for data analysis and generating figures were programmed by the authors and can be requested from the first author.

Received: 18 June 2021; Accepted: 3 February 2022;
Published online: 10 March 2022

### References

1. Green, H. W. & Houston, H. The mechanics of deep earthquakes. *Annu. Rev. Earth Planet. Sci.* **23**, 169–213 (1995).

2. Griggs, D. & Handin, J. *Rock Deformation (A Symposium)*. (Geological Society of America, 1960).

3. Kanamori, H. Frictional melting during the rupture of the 1994 Bolivian earthquake. *Science* **279**, 839–842 (1998).

4. Ogawa, M. Shear instability in a viscoelastic material as the cause of deep focus earthquakes. *J. Geophys. Res.* **92**, 13801–13810 (1987).

5. Meade, C. & Jeanloz, R. Deep-focus earthquakes and recycling of water into the earth's mantle. *Science* **252**, 68–72 (1991).

6. Raleigh, C. B. & Paterson, M. S. Experimental deformation of serpentinite and its tectonic implications. *J. Geophys. Res.* **70**, 3965–3985 (1965).

7. Kirby, S. H., Durham, W. B. & Stern, L. A. Mantle phase changes and deep earthquake faulting in subducting lithosphere. *Science* **252**, 216–225 (1991).

8. Kao, H. & Liu, L. G. A hypothesis for the seismogenesis of a double seismic zone. *Geophys. J. Int* **123**, 71–84 (1995).

9. Schubnel, A. et al. Deep-focus earthquake analogs recorded at high pressure and temperature in the laboratory. *Science* **341**, 1377–1380 (2013).

10. Green, H. W. Shearing instabilities accompanying high-pressure phase transformations and the mechanics of deep earthquakes. *Proc. Natl Acad. Sci. USA* **104**, 9133–9138 (2007).

11. Ferrand, T. P. et al. Dehydration-driven stress transfer triggers intermediate-depth earthquakes. *Nat. Commun.* **8**, 1–11 (2017).

12. Kita, S. & Ferrand, T. P. Physical mechanisms of oceanic mantle earthquakes: comparison of natural and experimental events. *Sci. Rep.* **8**, 1–11 (2018).

13. Ferrand, T. P. Neither antigorite nor its dehydration is “metastable”. *Am. Mineral.* **104**, 788–790 (2019).

14. Ferrand, T. P. & Manea, E. F. Dehydration-induced earthquakes identified in a subducted oceanic slab beneath Vrancea, Romania. *Sci. Rep.* **11**, 1–9 (2021).

15. Wiens, D. A., McGuire, J. J. & Shore, P. J. Evidence for transformational faulting from a deep double seismic zone in Tonga. *Nature* **364**, 790–793 (1993).

16. Iidaka, T. & Furukawa, Y. Double seismic zone for deep earthquakes in the Izu-Bonin subduction zone. *Science* **263**, 1116–1118 (1994).

17. Kawakatsu, H. & Yoshioka, S. Metastable olivine wedge and deep dry cold slab beneath southwest Japan. *Earth Planet. Sci. Lett.* **303**, 1–10 (2011).

18. Shen, Z. & Zhan, Z. Metastable olivine wedge beneath the Japan sea imaged by seismic interferometry. *Geophys. Res. Lett.* **47**, 1–8 (2020).

19. Wang, Y. et al. A laboratory nanoseismological study on deep-focus earthquake micromechanics. *Sci. Adv.* **3**, e1601896 (2017).

20. Shi, F., Wang, Y., Zhang, J., Yu, T. & Zhu, L. Instability induced by orthopyroxene phase transformation and implications for deep earthquakes below 300 km depth. *AGU Fall Meeting Abstracts* **2017**, S43B–S40841 (2017).

21. Kirby, S. H., Stein, S., Okal, E. A. & Rubie, D. C. Metastable mantle phase transformations and deep earthquakes in subducting oceanic lithosphere. *Rev. Geophys.* **34**, 261–306 (1996).

22. Däßler, R. & Yuen, D. A. The metastable olivine wedge in fast subducting slabs: constraints from thermo-kinetic coupling. *Earth Planet. Sci. Lett.* **137**, 109–118 (1996).

23. Devaux, J. P., Schubert, G. & Anderson, C. Formation of a metastable olivine wedge in a descending slab. *J. Geophys. Res. Solid Earth* **102**, 24627–24637 (1997).

24. Zhan, Z. Gutenberg–Richter law for deep earthquakes revisited: a dual-mechanism hypothesis. *Earth Planet. Sci. Lett.* **461**, 1–7 (2017).

25. Fyfe, W. S., Turner, F. J. & Verhoogen, J. *Metamorphic reactions and metamorphic facies*. Vol. 73 (Geological Society of America MEMOIRS, 1959).

26. Rubie, D. C. Role of kinetics in the formation and preservation of eclogites. In *Eclogite facies Rocks, Chapter 5*, 111–140 (1990).

27. Ferrand, T. P. Seismicity and mineral destabilizations in the subducting mantle up to 6 GPa, 200 km depth. *Lithos* **334–335**, 205–230 (2019).

28. Boneh, Y. et al. Intermediate-depth earthquakes controlled by incoming plate hydration along bending-related faults. *Geophys. Res. Lett.* **46**, 3688–3697 (2019).

29. Incel, S. et al. Laboratory earthquakes triggered during eclogitization of lawsonite-bearing blueschist. *Earth Planet. Sci. Lett.* **459**, 320–331 (2017).

30. Gasc, J. et al. Faulting of natural serpentinite: Implications for intermediate-depth seismicity. *Earth Planet. Sci. Lett.* **474**, 138–147 (2017).

31. Wyss, M., Klein, F., Nagamine, K. & Wiemer, S. Anomalously high b-values in the South Flank of Kilauea volcano, Hawaii: Evidence for the distribution of magma below Kilauea's East rift zone. *J. Volcanol. Geotherm. Res.* **106**, 23–37 (2001).

32. Mori, J. & Abercrombie, R. E. Depth dependence of earthquake frequency-magnitude distributions in California: Implications for rupture initiation. *J. Geophys. Res. Solid Earth* **102**, 15081–15090 (1997).

33. Wiens, D. A. & Gilbert, H. J. Effect of slab temperature on deep-earthquake aftershock productivity and magnitude-frequency relations. *Nature* **384**, 153–156 (1996).

34. Felzer, K. R. Calculating the Gutenberg-Richter b value. *AGU Fall Meeting abstracts* **2006**, S42C–08 (2006).

35. Marzocchi, W. & Sandri, L. A review and new insights on the estimation of the b-value and its uncertainty. *Ann. Geophys.* **46**, 1271–1282 (2003).

36. Okal, E. A. & Kirby, S. H. Frequency-moment distribution of deep earthquakes; implications for the seismogenic zone at the bottom of slabs. *Phys. Earth Planet. Inter.* **92**, 169–187 (1995).

37. Pacheco, J. F., Scholz, C. H. & Sykes, L. R. Changes in frequency-size relationship from small to large earthquakes. *Nature* **355**, 71–73 (1992).

38. Wiemer, S. & Benoit, J. P. Mapping the b-value anomaly at 100 km depth in the Alaska and New Zealand subduction zones. *Geophys. Res. Lett.* **23**, 1557–1560 (1996).

39. Ranero, C. R., Phipps Morgan, J., McIntosh, K. & Relchert, C. Bending-related faulting and mantle serpentinization at the Middle America trench. *Nature* **425**, 367–373 (2003).

40. Shillington, D. J. et al. Link between plate fabric, hydration and subduction zone seismicity in Alaska. *Nat. Geosci.* **8**, 961–964 (2015).

41. Fujie, G. et al. Controlling factor of incoming plate hydration at the north-western Pacific margin. *Nat. Commun.* **9**, 1–7 (2018).

42. Hirano, N. et al. An Paleogene magmatic overprint on Cretaceous seamounts of the western Pacific. *Island Arc* **30**, 1–16 (2021).

43. Tsuji, T., Nakamura, Y., Tokuyama, H., Coffin, M. F. & Koda, K. Oceanic crust and Moho of the Pacific Plate in the eastern Ogasawara Plateau region. *Island Arc* **16**, 361–373 (2007).

44. Kanamori, H. Quantification of Earthquakes. *Nature* **271**, 411–414 (1978).

45. Kanamori, H. The energy release in Great Earthquakes. *J. Geophys. Res.* **82**, 2981–2987 (1977).

46. Thielmann, M., Rozel, A., Kaus, B. J. P. & Ricard, Y. Intermediate-depth earthquake generation and shear zone formation caused by grain size reduction and shear heating. *Geology* **43**, 791–794 (2015).

47. Green, H. W. II, Shi, F., Bozhilov, K., Xia, G. & Reches, A. Z. Phase transformation and nanometric flow cause extreme weakening during fault slip. *Nat. Geosci.* **8**, 484–489 (2015).

48. Thielmann, M. Grain size assisted thermal runaway as a nucleation mechanism for continental mantle earthquakes: Impact of complex rheologies. *Tectonophysics* **746**, 611–623 (2018).

49. Ohtani, E., Litasov, K., Hosoya, T., Kubo, T. & Kondo, T. Water transport into the deep mantle and formation of a hydrous transition zone. *Phys. Earth Planet. Inter.* **143**, 255–269 (2004).

50. Ishii, T. & Ohtani, E. Dry metastable olivine and slab deformation in a wet subducting slab. *Nat. Geosci.* **14**, 526–530 (2021).

51. Green, H. W., Chen, W. P. & Brudzinski, M. R. Seismic evidence of negligible water carried below 400-km depth in subducting lithosphere. *Nature* **467**, 828–831 (2010).

52. Zhang, J., Green, H. W., Bozhilov, K. & Jin, Z. Faulting induced by precipitation of water at grain boundaries in hot subducting oceanic crust. *Nature* **428**, 633–636 (2004).

53. Ferrand, T. P., Nielsen, S., Labrousse, L. & Schubnel, A. Scaling seismic fault thickness from the laboratory to the field. *J. Geophys. Res. Solid Earth* **126**, e2020JB020694 (2021).
54. John, T. et al. Generation of intermediate-depth earthquakes by self-localizing thermal runaway. *Nat. Geosci.* **2**, 137–140 (2009).

55. Zhang, H., Wang, F., Myhill, R. & Guo, H. Slab morphology and deformation beneath Izu-Bonin. *Nat. Commun.* **10**, 1–8 (2019).

56. Argus, D. F., Gordon, R. G. & DeMets, C. Geologically current motion of 56 plates relative to the no-net-rotation reference frame. *Geochem. Geophys. Geosyst.* **12**, Q11001 (2011).

57. Jain, A. K., Duin, R. P. W. & Mao, J. Statistical pattern recognition: a review. *IEEE Trans. Pattern Anal. Mach. Intell.* **22**, 4–37 (2000).

58. Pan, W., Shen, X. & Liu, B. Cluster analysis: unsupervised learning via supervised learning with a non-convex penalty. *J. Mach. Learn. Res.* **14**, 1865–1889 (2013).

59. Steinley, D. K-means clustering: a half-century synthesis. *Br. J. Math. Stat. Psychol.* **59**, 1–34 (2006).

60. Likas, A., Vlassis, N. & J. Verbeek, J. The global k-means clustering algorithm. *Pattern Recognit.* **36**, 451–461 (2003).

61. Ahmed, M., Seraj, R. & Islam, S. M. S. The k-means algorithm: a comprehensive survey and performance evaluation. *Electronics* **9**, 1–12 (2020).

62. Dzwinel, W. et al. Nonlinear multidimensional scaling and visualization of earthquake clusters over space, time and feature space. *Nonlinear Process. Geophys.* **12**, 117–128 (2005).

63. Rehman, K., Burton, P. W. & Weatherill, G. A. K-means cluster analysis and seismicity partitioning for Pakistan. *J. Seismol.* **18**, 401–419 (2014).

64. Weatherill, G. & Burton, P. W. Delineation of shallow seismic source zones using K-means cluster analysis, with application to the Aegean region. *Geophys. J. Int.* **176**, 565–588 (2009).

65. Uchide, T. & Imanishi, K. Underestimation of microearthquake size by the magnitude scale of the Japan Meteorological Agency: influence on earthquake statistics. *J. Geophys. Res. Solid Earth* **123**, 606–620 (2018).

66. Aki, K. Maximum likelihood estimate of b in the formula log N= a·bM and its confidence limits. *Bull. Earthq. Res. Inst.* **43**, 237–239 (1965).

67. Utsu, T. A statistical significance test of the difference in b-value between two earthquake groups. *J. Phys. Earth* **14**, 37–40 (1966).

68. Bender, B. Maximum likelihood estimation of b values for magnitude grouped data. *Bull. Seismol. Soc. Am.* **73**, 831–851 (1983).

69. Guo, Z. Q. & Ogata, Y. Statistical relations between the parameters of aftershocks in time, space, and magnitude. *J. Geophys. Res. Solid Earth* **102**, 2857–2873 (1997).

70. Tinti, S. & Mulargia, F. Confidence intervals of b values for grouped magnitudes. *Bull. Seismol. Soc. Am.* **77**, 2125–2134 (1987).

## Acknowledgements

We thank the Japan Meteorological Agency for the earthquake catalog and Prof. Zhongwen Zhan for the constructive discussion. We acknowledge the Institute for Cyber-Enabled Research (ICER) at Michigan State University and the Extreme Science and Engineering Discovery Environment (XSEDE supported by the NSF grant ACI-1053575) for providing the high-performance computing resources. The map-view figure is produced using the Generic Mapping Tools (GMT). T.P.F. acknowledges the *Labex VOLTAIRE* (Orléans University, France) and the *Alexander von Humboldt*

Foundation (Germany). This research was supported by the NSF grant 1802247 (MinChen's grant).

## Author contributions

M.C. triggered the international collaboration. G.L.M. performed the data analysis. T.P.F. proposed interpretations based on a review of the literature, provided advice on the presentation of the data, and drew the explanatory diagrams. G.L.M., T.P.F., and M.C. wrote the original version of the paper. G.L.M., T.P.F., and J.L. implemented the reviewer's comments and suggestions. B.Z. helped G.L.M. on the coding and Z.X. did the calculation of the eigenvalues for the spectral clustering analysis. We dedicate the final version of this work to our late colleague M.C.

## Funding

Open Access funding enabled and organized by Projekt DEAL.

## Competing interests

The authors declare no competing interests.

## Additional information

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s43247-022-00377-x.

**Correspondence** and requests for materials should be addressed to Thomas P. Ferrand.

**Peer review information** *Communications Earth & Environment* thanks the anonymous reviewers for their contribution to the peer review of this work. Primary Handling Editors: Luca Dal Zilio, Joe Aslin. Peer reviewer reports are available.

**Reprints and permission information** is available at http://www.nature.com/reprints

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_8.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2022