Variations in the crustal structure and strength of plate coupling along the Ryukyu subduction zone

Abstract

The Ryukyu trench-arc-back arc system is part of the subduction margins of the Philippine Sea plate. Previous studies have indicated that several geophysical and geological characteristics reveal significant variations (including convergent rate, topography, subducting slab angle etc.) along this subduction system. In addition, the strength of plate coupling and the potential of large earthquake occurrence in the Ryukyu subduction zone have been major subjects of debate for decades. To gain new insights into the spatial variations in the crustal structure and strength of plate coupling along the Ryukyu subduction zone, in the present study, based on three P-wave seismic velocity profiles, we construct density models for 2-D gravity modeling. Then, we estimate the mantle lithosphere buoyancy (H_m) using these three density models to determine the strength of plate coupling between the subducting Philippine Sea plate and the overriding Eurasian plate, which could provide information for evaluating large earthquakes potential. 2-D gravity modeling results reveal that oceanic plateaus and/or submarine ridges with obviously less dense and thick oceanic crust are subducting in the northern and central parts of the Ryukyu Trench, which could increase the slab buoyancy in these regions. The H_m results indicate that the strength of plate coupling is almost weak in the north and is relatively strong in the central Ryukyu subduction zone.

Introduction

It is well-known that major devastating earthquakes generally occur along the plate interface in subduction zones, which are boundaries, where two plates converge and one plate dives beneath the other. According to the historical seismic catalog, the observed maximum magnitudes of subduction zone earthquakes are highly variable worldwide. For decades, scientists have made great efforts on better understanding the structural characteristics in each subduction zone to evaluate its seismic potential. However, subduction zone earthquakes are still the most impactful natural hazards on Earth and often cause great harm to human life and property. The Ryukyu Trench is part of the circum-Pacific seismic belt (Fig. 1), which is the world's greatest earthquake belt. Historically, no large earthquakes (M > 8) have been recorded in the last 300 years in this region (Ando et al. 2009); thus, the seismic coupling has been interpreted to be weak (Peterson and Seno 1984; Scholz and Campos 2012). However, Lin et al. (2014) suggested a high possibility of megathrust earthquakes in the Ryukyu subduction zone because of the tectonic similarities between the Ryukyu and Sumatra subduction zones. In addition, based on the dating of tsunami deposits, Ando et al. (2018) suggested that large tsunamis (such as the 1771 event) may occur at an average interval of approximately 600 years in the south Ryukyu subduction zone. Compared with the recurrence interval of large earthquakes, seismic observation alone is obviously limited to assessing the occurrence of great subduction zone earthquakes (McCaffrey 2008). Consequently, other independent approaches are worthy and necessary.

Fig. 1

a Map showing the locations of major subduction zones of the Circum-Pacific belt. The earthquake data, which show events with magnitudes larger than 8.0 (yellow dots indicate 8 < M < 9), were obtained from the GCMT project (1976 ~ 2017) data catalog. The red dots indicate M9.0 or larger earthquakes that occurred between 1700 and 2020. The thick green lines indicate the locations of subduction zones. b Regional map of the Ryukyu subduction zone. Black arrows illustrate the motion of the Philippine Sea plate with respect to the Eurasian plate. Red diamonds indicate the 1911 earthquake near Amami Island and the 1771 large tsunami off Ishigaki Island. Thick black dashed line indicates volcanic front (identified from Arai et al. (2018)). Earthquakes data were obtained from the Global CMT project (1976 ~ 2017) data catalog. EU Eurasian plate; GR Gagua Ridge; LOFZ Luzon–Okinawa fracture zone; PSP Philippine Sea plate; T Taiwan; WPB West Philippine Basin

Several studies have focused on the plate coupling conditions of the Ryukyu subduction zone and have further evaluated the potential of large earthquake occurrence. For example, Igarashi (2010) proposed weak plate coupling and a low potential for large interplate earthquake occurrence in the Ryukyu subduction zone by estimating slip rates of small repeating earthquakes. Arai et al. (2016) considered that the plate interface in the southern Ryukyu Trench is dominated by slow earthquakes (including low-frequency earthquakes and slow slip events), which thus lacks a typical locked zone. In conflict with the inference of a weakly coupled condition, Ando et al. (2012) suggested the strong likelihood of a coupled portion in the Ryukyu subduction zone (by very low-frequency earthquake observations); Hsu et al. (2012) considered that the plate interface could be fully locked in the southern Ryukyu subduction zone (using GPS data); Hsu et al. (2013) proposed that a strong plate coupling condition due to the subducted Luzon–Okinawa fracture zone (LOFZ) could resist subduction; Tadokoro et al. (2018) proposed a strongly coupled state in the central portion (seafloor crustal deformation measurement results). Kano et al. (2021) proposed that Mw of 7.5 or larger earthquakes could occur in the southern Ryukyu Trench (from GNSS data analyzing). In brief, the strength of plate coupling and the potential for large earthquake occurrence in the Ryukyu subduction zone are still controversial.

As an independent method from seismic observations, considering the relationship between the topography and the interaction between overriding and subducting plates, Gvirtzman and Nur (1999a) proposed that variations in mantle lithosphere buoyancy (H_m) across the subduction zone could reveal changes in the strength of plate coupling. Intuitively, in subduction zones, if the subducting and overriding plates become locked together (strong coupling), stress is more easily accumulated, and more energy can be released to produce earthquakes. Doo et al. (2020) have discussed their relationships in several subduction zones. Their conclusion indicates that strong plate coupling correlates well with the occurrence of large earthquakes, whereas weak plate coupling can constrain the potential occurrence of large earthquakes. Using this method, Hsu (2001) and Doo et al. (2018) obtained a relatively strong plate coupling status in the southern Ryukyu subduction zone. However, the degree of coupling along a long subduction zone interface can vary (Lamb and Davis 2003; Lamb 2006). In addition, previous studies pointed out that several geophysical and geological characteristics vary significantly along the entire Ryukyu subduction zone. To understand the lateral variation in the strength of plate coupling along the whole Ryukyu Trench, in this study, we process 2-D gravity modeling along three velocity transects cutting across the northern and central Ryukyu subduction zone (Fig. 2) and then calculate their corresponding H_m. 2-D gravity modeling results can provide new constraints on deep crustal structures, where the resolution of the velocity model is poor; then, we can better understand the variations in subsurface structures. H_m results can provide the information necessary for understanding the strength of plate coupling and assessing the large earthquake potential in the Ryukyu subduction zone.

Fig. 2

Free-air gravity anomaly map (Sandwell et al. 2014) of the Ryukyu subduction zone. Thick black lines represent the profiles that we used to process 2-D gravity modeling and estimate the buoyancy of mantle lithosphere in this study. Black arrows illustrate the motion of the Philippine Sea plate with respect to the Eurasian plate. EU Eurasian plate; GR Gagua Ridge; LOFZ Luzon–Okinawa fracture zone; PSP Philippine Sea plate; T Taiwan

Tectonic background

Along the Ryukyu Trench, the Philippine Sea plate (PSP) subducts northwestward beneath the Eurasian plate (EU), forming a back-arc basin, namely, the Okinawa Trough (OT). The convergence rate is approximately 7 cm/yr in the north and progressively increases to approximately 13 cm/yr in the south (Argus et al. 2011). The Ryukyu Trench trends NE–SW in the north and turns approximately E–W at its southwestern end. Traditionally, the 1200 km-long Ryukyu arc-trench system can be divided into three parts (northern, central, and southern) by two significant bathymetric depressions running in the arc-perpendicular direction: the Tokara Channel in the north and the Kerama Gap in the south (Shiono et al. 1980; Kuramoto and Konishi 1989). The incoming PSP contains several significant topographic features along the Ryukyu Trench: in the northern and central parts, where the subduction direction is almost perpendicular to the trench orientation, several volcanic ridges (including the Kyushu-Palau Ridge, Oki-Daito Ridge, Daito Ridge, and Amami Plateau) with a large crustal thickness (Nishizawa et al. 2017) are subducting. A clear lineament feature (trending approximately NE–SW) cuts across the West Philippine Basin, the LOFZ, whose northern part has already been subducted beneath the Ryukyu Trench (Hsu et al. 2013). In the southern part of the subduction zone, the N–S trending Gagua Ridge extends northward and subducts beneath the Ryukyu Trench along 123°E (Dominguez et al. 1998; Schnürle et al. 1998; Doo et al. 2021), which separates the seafloor spreading direction and the age of the Huatung basin from that of the oceanic crust of the PSP (Hilde and Lee 1984; Deschamps et al. 2000; Sibuet et al. 2002; Doo et al. 2015a). Overall, the topographic features are flat in the south and rough in the central and northern regions. The maximum depths in the three segments along the trench vary, increasing from north (< 5000 m) to south (> 7000 m). To the deeper part, the dip angles of the subducting PSP also reveal different features; in the northern part of the Ryukyu Trench, the slab angle is steeper at deeper depths (> 50 km). In contrast, the slab dip angles in the central and southern parts are relatively gentle (Iwasaki et al. 1990; Kodaira et al. 1996).

In the overriding EU, the crustal thickness is flat along the northern and central Ryukyu Arc and OT and changes abruptly along the southern Ryukyu Arc and OT (Nakamura and Umedu 2009). The OT shows an obvious variation in seafloor topography: the northern OT is a gentle depression with shallow bathymetry (< 1000 m); water depths are within 1000–2000 m in the central OT, while the southern OT consists of a relatively narrow basin with a maximum depth greater than 2000 m (Fig. 1). For such significant characteristics, the formation mechanism of the Tokara Channel and Kerama Gap (geographic boundaries) and the relations to the evolution of the OT are interesting but are still not well-understood. In addition, similar to the variation in the convergence rate, the rifting rate along the OT increases southward from ~ 2 cm/yr in the north to ~ 5 cm/yr in the south (Argus et al. 2011). Better imaging of subsurface structures would be a critical step for understanding how such complex tectonic processes work in the Ryukyu subduction zone.

Buoyancy of mantle lithosphere (H_m)

Different from previous studies that used seismic and/or GPS data to evaluate plate coupling conditions in subduction zones, Gvirtzman and Nur (1999a) provided a method by estimating the buoyancy of mantle lithosphere to determine the strength of plate coupling between the subducting and overriding plates. We provide a brief introduction to this method in the Appendix. For details on this method, readers can refer to the literature (Gvirtzman and Nur 1999a; 1999b; 2001).

P-wave seismic velocity models and 2-D gravity modeling

To determine H_m, the lithospheric density model is needed. We applied a 2-D gravity modeling technique to establish the subsurface density model (e.g., Doo et al. 2015b; 2016; 2018; 2021). To reduce the nonunique problem, the seismic velocity structure is crucial and necessary for 2-D gravity modeling processing. Thus, three approximately NW–SE trending P-wave seismic velocity profiles (from a joint ocean bottom seismometer (OBS) and multichannel reflection seismic experiments) were used in this study (Fig. 2). Profiles ECr11 and 2-A&2-B cut across the northern and central Ryukyu subduction zone, respectively. P-wave seismic velocity models of profiles ECr11, RK02, and 2-A were presented by Nishizawa et al. (2017), Arai et al. (2017), and Kodaira et al. (1996), respectively. The resolution of the velocity model of profile ECr11 (Fig. 3a) is relatively good. The locations of the subducted slab and the Moho depths can be well-identified. For profile RK02, which is roughly across the Tokara Channel, the locations of the deep slab and Moho depths are relatively uncertain (Fig. 4a). In the central Ryukyu subduction zone, the OBS data were acquired along profiles 2-A and 2-B in 1988. Kodaira et al. (1996) only provided the velocity model of only profile 2-A (Fig. 5a). To extend the velocity image northwestward (profile 2-B), we applied the Slab 2.0 model (Hayes et al. 2018) to define the location of the subducted slab. Because Kodaira et al. (1996) provided only layer structures (which is different from profiles ECr11 and RK02), we thus extended layer boundaries directly to construct the initial velocity model of profile 2-B (Fig. 5a). Through processing forward gravity modeling, we modified the initial velocity model and then provided a possible velocity model of profile 2-B.

Fig. 3.

2-D gravity modeling result of profile ECr11. a P-wave seismic velocity model of Nishizawa et al. (2017) for profile ECr11. b Observed (Sandwell et al. 2014) and synthetic gravity anomalies. RMS: root mean square. Synthetic gravity anomaly is the result of the density model shown in (c)

Fig. 4.

2-D gravity modeling results of profile RK02. a P-wave seismic velocity model of Arai et al. (2017) for profile RK02. Black dots indicate the location of OBSs. b Observed (Sandwell et al. 2014) and synthetic gravity anomalies. Synthetic I is the result of the initial density model. To fit the observed gravity data, we modify the Moho depth, increase density of the arc crust, and add a subducted plateau (blue dashed lines). Then, we obtain a better fitting gravity anomaly (Synthetic II). RMS: root mean square. c Initial and final density models

Fig. 5.

2-D gravity modeling results of profile 2A&2B. a P-wave seismic velocity model of Kodaira et al. (1996). Dashed red line indicates the location of the subducted slab (Hayes et al. 2018). b Observed (Sandwell et al. 2014) and synthetic gravity anomalies. Synthetic I is the result of the initial density model. To fit the observed gravity data, we modify the Moho depth (within 40–100 km), and reduce the densities of layer 2 (150–210 km) and the oceanic crust (270–296 km). Synthetic II is the result of the final density model. RMS: root mean square. c Initial and final density models

To process forward gravity modeling, first, we constructed subsurface structural geometries based on P-wave seismic velocity models. In areas, where the velocity contours are unclear, we extended the velocity contours smoothly as the layer boundaries. P-wave velocities were converted to densities using the velocity–density function (Eq. 1) of Brocher (2005):

$$ \rho (g/cm^{3} ) = 1.6612V_{p} - 0.4721V_{p}^{2} + 0.067V_{p}^{3} - 0.0043V_{p}^{4} + 0.000106V_{p}^{5} $$

This equation is valid for V_p values between 1.5 and 8.5 km/s. Subsequently, the initial 2-D density model was constructed for each profile (Figs. 3c, 4, 5c). Then, the gravity anomaly was calculated using the formula of Blakely (1995). If the synthetic anomaly could not adequately fit the observed gravity data (Sandwell et al. 2014), we slightly adjusted the densities (within a reasonable range) until the synthetic anomaly was fitted to the observed gravity anomaly. For details on the process of forward gravity modeling, readers are referred to the literature (Doo et al. 2015b; 2016). Readers may be concerned about the nonunique problem of 2-D gravity modeling, P-wave seismic velocity models could provide good constraints on subsurface structural geometry and density information that help to greatly reduce this problem.

Results

2-D gravity modeling

For profile ECr11, the synthetic gravity anomaly (green dashed line) fits the observed data well (Fig. 3b). Consequently, we do not modify the density and geometries of the velocity model. For profile RK02, a relatively large difference between synthetic I (green dashed line) and observed data is revealed in three parts: 20–120, 120–190, and 230–290 km (Fig. 4b). We modify the Moho depth (considering the resolution of velocity model and wavelength of gravity signal), increase density of the arc crust, and add a subducted plateau (Amami Plateau) for these three segments. Then, we obtain a better fitting gravity anomaly (blue dashed line shown in Fig. 4b). For profile 2-A&2-B (Fig. 5), to fit the observed gravity data well, we modify the Moho depth (within 40–100 km) according to the results of past studies in this area (Kodaira et al. 1996; Nishizawa et al. 2017) and reduce the densities of layer 2 (150–210 km) and the oceanic crust (270–296 km). From 2-D gravity modeling, even though the final result is a compendium map of the subsurface structural geometry, we find that the variation in each stratum is small along profile 2-B (Fig. 5).

In profile ECr11 (Fig. 3c), we find a typical thickness of oceanic crust subducting; in contrast, the thickness of the subducted oceanic crust of profile RK02 is thicker (~ 12 km), and the density is smaller than that of typical oceanic crust (Fig. 4c), which may be because the crust is influenced by the Amami Plateau. Arai et al. (2017) also interpreted this part to be a plateau. Basis on the modeling results of profiles ECr11 and RK02, we find that plateau could obviously result in large influence in crustal thickness and density. In the central part of the Ryukyu subduction system, topographically, the Oki–Daito Ridge does not extend to the trench (Fig. 1). However, in profile 2-A&2-B, the southeast end (270–296 km) shows thick oceanic crust with a low density (Fig. 5c), which may be influenced by the Oki–Daito Ridge. This result implies that the deep part of the plateaus and/or ridges is larger than that revealed in the seafloor and further indicates that the Oki–Daito Ridge is close to the trench. Gravity modeling results show that densities of submarine plateaus and/or ridges are obviously smaller than those of the lower crust of the overriding plate. Not only the conspicuously low density but also the thick crustal thickness of the incoming plateaus is expected to increase the slab buoyancy.

Strength of plate coupling

Gvirtzman and Nur (1999a) provided four typical cases to illustrate the variation characteristics of H_m and its corresponding plate coupling status. For the decoupled plate case (Calabria case shown in Fig. 6), the subducted slab leaves the overriding plate and deep material upwells into the corner between the plates; thus, the curve of H_m across the subduction zone would depict a sharp undulation over a short distance. In contrast, strong plate coupling (Andes case in Fig. 6), based on a large portion of the overriding plate is coupled to the subducted slab, the H_m curve across the subduction zone would reveal small variations. For the Kurile and Izu–Bonin cases, H_m variations are large over a wide distance, giving the definition of the intermediate plate coupling status. For this study, the H_m estimation results from the three profiles are shown in Fig. 6. Compared with the four cases provided by Gvirtzman and Nur (1999a), we find that the H_m variation pattern of profile ECr11 is similar but smaller in scale to the Calabria case, which indicates that the coupling status between the PSP and the EU is close to a weak state in this region. Profile RK02 cuts across the Tokara Channel, the boundary between the northern and central Ryukyu Arc; thus, it is difficult to distinguish whether this profile belongs to the north or central part of the Ryukyu subduction zone. Its minimum H_m value is larger than that of ECr11 at the plate boundary (the subducting plate suffer a strong suction), which indicates that the plate coupling status is stronger than that in profile ECr11. For profile 2-A&2-B, the variation of the H_m curve across the central Ryukyu subduction zone is relatively small, suggesting that the strength of plate coupling is relatively strong, and this result is consistent with the idea by Tadokoro et al. (2018).

Fig. 6

Map showing H_m curves. The curves of ECr11, RK02, and 2-A&2-B represent the results of estimating H_m in this study. The curves of the Andes, Calabria, Izu–Bonin and Kurile are modified from Gvirtzman and Nur (1999a). The Andes curve represents a case of strong plate coupling. The Calabria curve represents a case featuring weak plate coupling. The Izu–Bonin and Kurile cases display intermediate plate coupling status. The gray band zone (2.0 ± 0.5) indicates the variations in H_m along the North America passive margin, which represents the normal contribution of the mantle lithosphere to the Earth’s topography

Discussion

Behind the northern and central parts of the Ryukyu Arc, based on bathymetric data, a series of active volcanoes (volcanic front) were identified (thick black dashed line shown in Fig. 1; Arai et al. 2018). This linear distribution of active volcanoes starts from southern Kyushu and ends at the central Ryukyu Arc. Coincidently, approximately in the corresponding area, on the other side of the trench, large-scale plateaus and/or ridges are present. Arai et al. (2018) proposed that oceanic plateau subduction may increase the density of active volcanoes on the arc. Nakamura et al. (2003) considered that the subduction angle of the PSP may have caused differences in volcanism in the Ryukyu subduction zone. A similar case is that Yang et al. (1996) proposed a ridge subduction model to explain the double arc structure observed in the Bashi segment of the Luzon Arc (Manila subduction zone). Considering the mechanism and tectonic features, the model of Yang et al. (1996) may be comparable to that for explaining the phenomenon appearing in the Ryukyu subduction zone. We thus proposed that distinct differences in geophysical and geological features displayed in the Ryukyu subduction zone may result from large-scale oceanic plateaus and/or ridges subduction.

On the other hand, Doo et al. (2018) obtained a density profile across the southern Ryukyu Trench (profile AA' shown in Fig. 2) through 2-D gravity modeling. Comparing these density profiles (Fig. 7), we observed a thick sediment wedge with a high density (~ 2.55 g/cm³) in the north and thin thicknesses of the sediment wedge with a low density (~ 2.4 g/cm³) in the south. This result may imply that the amount of sediment supplied from the surrounding areas is larger in the north than that in the south. Correspondingly, the maximum water depth of the trench is shallower in the north and deeper in the south (Arai et al. 2017). Combining the bathymetric data with gravity modeling results, we considered that most northern and central parts of the Ryukyu subduction zone could be influenced by plateau and ridge subduction. Previous studies (Ballance et al. 1989; Yamazaki and Okamura 1989; von Huene and Scholl 1991) suggested that submarine ridge and/or plateau subduction could enhance subduction-related erosion at the base of the fore-arc wedge and cause subsidence of fore-arc slopes. In addition, according to OBS refraction data (at almost the same location as profile ECr11), Iwasaki et al. (1990) proposed that the origin of the sediment wedge may be oceanic. Their conclusions could well explain the high density of thick sediments that we observed in the northern and central parts of the fore-arc region.

Fig. 7

Density model of four profiles across the Ryukyu subduction zone. a Profile ECr11, b profile RK02, c profile 2-A&2-B and d profile AA' (modify from Doo et al. (2018))

According to these four density profiles (Fig. 7), overall, in the study area, the PSP oceanic crust shows normal oceanic crustal thickness (5–7 km) except for some regions that are covered by plateaus and/or ridges (Doo et al. 2015a; Eakin et al. 2015; Arai et al. 2016; Nishizawa et al. 2017). Gravity modeling results demonstrate that large-scale submarine ridges and plateaus could greatly increase the thickness and decrease the density of the PSP oceanic crust. In addition to the variation feature of the bathymetry, the crustal thickness of the OT has been significantly thinned from the central to southern segments. This along-trough crustal thickness variation is consistent with the southward increase in rifting rates (Argus et al. 2011).

Two other H_m profile analyses were performed by Hsu (2001) and Doo et al. (2018), and these analyses show a relatively strong plate coupling status in the southern Ryukyu subduction zone. Summarizing all the results, we can obtain that plate coupling status is relatively strong in the southern and central Ryukyu subduction zone, while plate coupling status is relatively weak in the northern segment. As mentioned above, potential large earthquake occurrence could be positively related to the strength of plate coupling. We thus conclude that the potential of large earthquake occurrence in the northern Ryukyu subduction zone may be lower. This result seems consistent with the historical seismic catalog. However, according to the topographic features and gravity modeling results, we find that large-scale plateaus and/or submarine ridges are or prior to subduction. Although the relationship between subducting features (such as fracture zones, plateaus, and ridges) and the occurrence of large earthquakes is complex and is not consistent among various subduction zones (Cloos 1992; Kato and Ando 1997; Schole and Small 1997; Kodaira et al. 2000; 2003; Abercrombie et al. 2001; Wang and Bilek 2011; Wang and Lin 2022), obviously, thick crust with less density could physically increase slab buoyancy and prevent slab subduction. In addition, Nishikawa and Ide (2014) suggested that slab buoyancy is an important control on the stress state and the earthquake size distribution in subduction zones. In the northern and central parts of the Ryukyu subduction zone, most of the incoming PSP is covered with large-scale plateaus and submarine ridges (Fig. 1). The effect of ridge subduction could apparently be more critical. Further tectonic investigations and earthquake monitoring in the northern Ryukyu subduction zone are still essential.

Conclusions

On the basis of three P-wave seismic velocity models, we obtained subsurface density structures across the northern and central Ryukyu subduction zone using 2-D gravity modeling. Combined with an additional density profile (across the southern Ryukyu Trench), we find large thicknesses with high-density sedimentary wedge in the north and small thicknesses with low-density sedimentary wedge in the south Ryukyu subduction zone. The origin of these high-density materials may come from the subducting plateaus and/or ridges. 2-D gravity modeling results also reveal that plateaus and/or ridges could greatly increase the thickness and decrease the density of the PSP oceanic crust. Several distinct differences in geophysical and geological features displayed in the Ryukyu subduction zone may result from large-scale oceanic plateaus and/or ridges subduction. The H_m results indicate that the strength of plate coupling in the northern Ryukyu subduction zone is relatively weak compared to other segments. For this result, the potential of large earthquake occurrence in the northern Ryukyu subduction zone may be lower.

Appendix: introduction of buoyancy of mantle lithosphere (Hm)

Under the condition of isostatic equilibrium, the mean elevation (ɛ) of a region is the sum of the buoyancies of the crust (${H}_{c}$) and the mantle (${H}_{m}$) lithosphere (Lachenbruch and Morgan 1990) as follows:

$$\upvarepsilon =\mathrm{a}({H}_{c}+{H}_{m}-{H}_{0})$$

(1)

a = 1 for $\upvarepsilon \ge 0$

$\mathrm{a}=\frac{{\rho }_{a}}{{\rho }_{a}-{\rho }_{w}}$ for $\upvarepsilon <0$

where ${\rho }_{a}$ and ${\rho }_{w}$ are the densities of the asthenosphere and sea water, respectively. ${H}_{c}$ and ${H}_{m}$ are defined as

$${H}_{c}=\frac{1}{{\rho }_{a}}\left({\rho }_{a}-{\rho }_{c}\right){L}_{c}$$

(2)

$${H}_{m}=\frac{1}{{\rho }_{a}}\left({\rho }_{a}-{\rho }_{m}\right){L}_{m}$$

(3)

where ${\rho }_{c}$ and ${\rho }_{m}$ are the densities and ${L}_{c}$ and ${L}_{m}$ are the thicknesses of the crust and mantle lithosphere, respectively. ${H}_{0}\approx 2.4 km$ is a reference constant for the buoyant height of sea level at a mid-ocean ridge in Eq. 1. H_c can be estimated according to the density model and crustal geometry; then, the changes in surface elevation are associated with the value of H_m (Eq. 1). However, in subduction zones, the variations in topography are influenced by the subducting slab (coupled or decoupled). The residual topography (observed topography minus the calculated contribution of the H_c), therefore, reflects the contribution from the buoyancy of mantle lithosphere and the forces that are pulled down by descending slabs (coupled) or uplifted by the ascending asthenospheric materials (decoupled). In other words, in subduction zones, the change of the residual topography depends on the status of plate coupling and can be revealed by the variation in H_m.