A Decade of Global Navigation Satellite System/Acoustic Measurements of Back-Arc Spreading in the Southwestern Okinawa Trough

Abstract

Long-term seafloor geodetic measurements are important for constraining submarine crustal deformation near plate boundaries. Here we present an integrated analysis of a decade of GNSS/acoustic data collected at a site 60 km to the east of northeast Taiwan near the axis of the Okinawa Trough back-arc basin. We obtained a time-series of horizontal and vertical positions based on 18 measurements from 2009 to 2019. These data reveal a southeastward movement at a rate of 43 ± 5 mm/yr since 2012 with respect to the Yangtze Plate. The horizontal motion can be explained by the clockwise rotation of the Yonaguni Block and northern Central Range. In addition, the vertical displacement of the transponder array shows rapid subsidence of 22 ± 9 mm/yr from 2012 to 2019. The fast subsidence rate and negative free-air gravity anomaly in this region indicate that crustal thinning is compensated mainly by surface deformation rather than upward migration of the Moho. Taking into account the offset in 2012 owing to the replacement of the transponder array, the horizontal position time series of our site are best explained by two linear lines with a slope change in July 2013. The timing of the velocity change coincides broadly with a change in the nearby seismicity rate and dike intrusion 150 km away from the site. Our results highlight the potential of seafloor geodesy in assessing temporal changes in deformation near the spreading center of the Okinawa Trough, which cannot be one using data from onland GNSS stations.

Introduction

The Okinawa Trough is an active back-arc basin (Sibuet et al., 1987). The axis of the basin has its southwestern end in the Yilan Plain of Taiwan and extends to the Beppu–Shimabara rift zone in Japan. The spreading rate is fastest (60 mm/yr) around the Yaeyama rift in the southwestern Okinawa Trough (Nishimura et al., 2004; Figure 1). Episodic opening of the rift has been reported around the Yaeyama rift. For example, dike intrusion was detected north of Yonaguni Island in 2002 (Nakamura and Kinjo, 2018), beneath the Yilan Plain in 2005 (Lai et al., 2009), and north of Iriomote Island in 2013 (Ando et al., 2015). In contrast to the episodic rifting, onland Global Navigation Satellite System (GNSS) measurements of the Ryukyu Arc (Lallemand and Liu, 1998; Nishimura et al., 2004) have shown an almost constant southward velocity. To investigate episodic opening of the Okinawa Trough, Chen et al. (2018) conducted GNSS/acoustic seafloor geodetic measurements immediately south of the spreading center for four years from 2012 to 2016, and suggested episodic opening of the rift valley occurred between the most active site of rifting and the Yilan Plain.

FIGURE 1

DeMets et al. (2010). Hypocenters of earthquakes larger than M5 in the same period are shown as solid circles and are color-coded according to depth. Earthquake source parameters were determined by the Japan Meteorological Agency (JMA) and Central Weather Bureau (CWB). The black box encloses the study area of Figure 4. The location of the OILN site is shown as a white circle. The velocity at the OILN site was estimated from the data acquired over a 4-years period in Chen et al. (2018). The vectors of the on-land GNSS stations were calculated from the daily time series by linear fitting after removing offsets caused by large earthquakes, antenna maintenances, and offsets over 30 mm even if the cause is not clear.

In this study, we analyzed data obtained using the seafloor benchmark unit installed by Chen et al. (2018) for three additional years (i.e., 2016–2019). In addition, we retrieved data from 2012 and earlier years that were discarded previously because of the replacement of transponders due to battery exhaustion, making it impossible to compare these data with later transponder data. Using these data acquired over the past decade, we determined the vertical and horizontal motions at this site and the associated uncertainties.

Data and Methods

We carried out 18 series of GNSS/acoustic geodetic measurements using a fishing boat or buoys (Figures 2A–C) from 2009 to 2019 at the seafloor site OILN (Figure 1), which includes the 8 series of measurements reported in Chen et al. (2018). We used the drifting measurement strategy proposed by Ikuta et al. (2008), in which the vessel drifts with the sea current and wind for hours with the boat engine shut down during the transmission and reception of acoustic signals. We repeated the drifting measurement with different initial vessel locations until the coverage was sufficient. The observation period for each campaign varied from 11 to 62 h, depending on the weather and sea conditions.

FIGURE 2

History of the Site

The OILN site comprises a triplet or quartet of mirror-response transponders that were installed at a depth of about 1300 m. The site was first built in October 2008 with a triplet of transponders equipped with short-lived batteries (Units 3, 4, and 9 in Table 1 and Figure 2E). In 2010, we deployed an additional transponder (Unit 12). Three more transponders were installed in 2012 (Units 15, 17, and 18) to replace the transponder triplet installed in 2008, as their batteries were exhausted between 2010 and 2012. Only Unit 12 was in operation before July 2012. The number of available seafloor transponders was three and four prior to and after 2012, respectively. The detailed history of the transponder replacement is listed in Table 1.

Data Acquisition

For 15 of the 18 field campaigns, measurements were conducted using the fishing boat Yilonglu-VI, whereas the other 3 campaigns were carried out using two towed buoys (Table 1; Figure 2). Each vessel was equipped with a transducer, a GNSS antenna for determining the position of the transducer, and a GNSS antenna triplet (the boat and buoy I) or a gyrocompass (buoy II) for determining the attitude of the vessel (Figure 2). The positions of the seafloor transponders were measured by acoustic ranging from the onboard transducer, whose position was determined by combining the GNSS antenna position and attitude of the vessel. We solved the kinematic 3D position of the GNSS antenna every 0.2 s using GrafNav Ver.8.60. Two onshore GNSS stations, SCHN and S101, located 110 km SW and 90 km W of the OILN site, respectively, were used as fiducial sites to solve the kinematic GNSS antenna position. Given the signal condition of the GNSS antenna on the top of the transducer pole was sometimes not optimal, owing to obstruction by onboard structures above the antenna, on several campaigns Chen et al. (2018) used a GNSS antenna at the bow of the vessel as the main antenna instead of the on-pole antenna (Figure 2C). In this study, we improved the kinematic positioning of the on-pole antenna by using the bow antenna as a nearby moving fiducial site. This very short baseline analysis improved the yield of kinematic solutions significantly and allowed the on-pole antenna to always be used as a reference to calculate the onboard transducer position. The length between the antenna and transducer was determined with an accuracy of 1 cm. As a result, the repeatability of the estimated vertical position was improved relative to the result of Chen et al. (2018).

According to Chen et al. (2018), the error propagation from the unbiased observation error to the benchmark location is approximately , where ΔG represents the error on the kinematic GNSS positioning, Δθ is the uncertainty of the vessel attitude in 3D orthogonal coordinates, r is the relative position of the onboard transducer with reference to the onboard GNSS antenna, Δc is the travel-time error measured as a distance, and N is the number of acoustic shots. ΔG is 8 cm in the vertical direction and 2 cm in the horizontal direction (Ikuta et al., 2008). Δθ×r is 2 × 2.7 m for Buoy II, 0.8 × 2.6 m for Buoy I, and 0.4 × 6.0 m for the fishing boat (Chen et al., 2018). Δt is 0.04 ms for all systems (Ikuta et al., 2008). Based on these values, the resulting error varies from 1.3 to 2.0 mm for a number of acoustic measurements between 3000–12,000. This error is much smaller than the actual repeatability of the estimated seafloor benchmark positions. The error on the benchmark positioning is larger than that expected from the unbiased observation error, which may be due to unexpected biased error in the observations, including kinematic GNSS positioning, as well as a discrepancy between the observational model and actual signal propagation.

Analytical Methods

We solved the position of the seafloor transponders using a tomographic technique based on Chen et al. (2018). The position of the seafloor transponders, as well as the sound speed variation and onboard configuration, were solved using the penalized least squares method. We used the following equation:where T_i is the ith acoustic travel-time (ms) obtained at time t_i (s), f (r₀,r_i, a (t_i)) is the travel-time predicted by ray tracing from the position of the seafloor transponders r_0j in the jth campaign and the onboard transducer r_i, a(t) is the temporal varying coefficient (no units) of the speed of sound in seawater represented by superposition of cubic B-spline functions (Barnhill and Riesenfeld, 1974), and μ_j is the hyper parameter that balances the smoothness of a(t) with the travel-time residual. We selected the best hyper parameter for each campaign using the Akaike’s Bayesian information criterion (ABIC), based on Chen et al. (2018). The resulting hyper parameters range from 10^4.00 to 10^4.16, although the values of Chen et al. (2018) ranged from 10^3.67 to 10^4.00. Hyper parameter dependence of the seafloor benchmark position solutions is evaluated in Appendix. The total number of campaigns in this study was 17, because two measurements conducted in 2012 were merged into one campaign because of the limited amount of data.

To solve Eq. 1, we used two realistic assumptions. Firstly, in estimating r₀, we assumed that all the transponders were fixed on a rigid body that does not rotate, and moved from campaign to campaign in the same geometric configuration. This rigid body assumption allows us to estimate the motion of the rigid body between campaigns, even if some of the transponders are different, provided there was at least one common transponder in each campaign. Secondly, the temporal variation of the speed of sound was assumed to be controlled by a horizontal layered structure:where z is the depth and V₀(z) is the speed of sound obtained by the CTD measurements.

The configuration between the onboard GNSS antenna and transducer was not measured accurately, but was solved from the acoustic ranging because it was difficult to keep the same instrument configuration on the fiberglass fishing boat between the different campaigns. Although the acoustic ranging has an approximately linear trade-off between the vertical component of the seafloor transponder and the vertical position of the onboard transducer relative to the GNSS antenna, the horizontal onboard configuration of the transducer can be estimated robustly by the acoustic ranging. Given that the onboard GNSS antenna is on top of the transducer pole, we measured the length of the pole L with an accuracy of <1 cm. When the horizontal distance of the transducer from the GNSS antenna is solved as h₀, the vertical distance is . We first solved Eq. 1 by assuming the pole length L, and then resolved it iteratively. This improvement in the on-pole antenna solution allowed us to obtain more accurate vertical seafloor position than Chen et al. (2018).

Results

Figure 3 shows a time-series of the 3D position of the OILN site relative to the Yangtze Plate from 2009 to 2019. The horizontal components show a southeastward motion, but there is a discontinuity between the initial four and the latter campaigns. The campaigns prior to and post-2012 have only one common transponder (Unit 12). Furthermore, the amount of data from the third and fourth campaigns (June 2010 and August 2011), which were measured using the common transponder (Unit 12) before 2012, is limited (Table 1). For this reason, we first focused on the post-2012 results without considering the pre-2012 data. In addition, we discarded the positions from the two campaigns conducted in 2010 and April 2014, due to the limited acoustic data (<3,000), as compared with the other campaigns (average = 7,597).

FIGURE 3

Chen et al. (2018). Thin black and gray lines show the positions of the GNSS sites on Gueishang (GS20) and Yonaguni (0499, 1211) islands, respectively. Station 1211 on Yonaguni Island was built ∼800 m from site 0499 in March 2014, and its absolute position on the diagram is not meaningful. Open and solid purple circles show the campaigns before and after transponder replacement, respectively. Circles with brackets show the positions estimated from acoustic data with <3000 shots. (B) E–W component. (C) Up–down component.

The average rate of movement since 2012 is 41.2 ± 4.8 and 13.5 ± 3.7 mm/yr in the southward and eastward direction, respectively (Figures 3A,B). This result is considerably different from the average velocities of 50.3 ± 11.2 and 34.8 ± 7.2 mm/yr in the southward and eastward directions, respectively, obtained by Chen et al. (2018) from measurements over 4 years (2012–2016). This discrepancy is significant, even taking into account the 1σ error, and may reflect a temporal change in the motion (Temporal Variation in Horizontal Motion). The position of the benchmark for each campaign also differs in the present and previous studies. The vertical positions reveal there has been 22.0 ± 8.9 mm/yr of subsidence over 7 years (Figure 3C). The scatter of data from the individual campaigns from a linear trend is 33 and 59 mm, which is expressed as the root-mean-square (rms) deviation of the horizontal and vertical components, respectively.

Discussion

Observation Errors

The horizontal velocity of the OILN site from 2012 to 2019 has an accuracy of <5 mm/yr, and the positions of the individual campaigns have a scatter of 35.6 and 27.6 mm (rms deviations from the linear trends of the N–S and E–W components, respectively). This error is comparable with the value of 20–30 mm obtained by the Japan Coast Guard (JCG) group (e.g., Yokota et al., 2016). Although the errors result from a combination of numerous factors, our observations appear to have some limitations compared with those of the JCG group. One involves the uncertainty on the position of the acoustic transducer relative to the GNSS antenna, which changed slightly from campaign to campaign. The fishing boat used in this study is made of fiberglass and, as such, rigging is not possible. This caused difficulties in setting the sensor jigs at the same positions with centimeter accuracy. As a result, the configuration between the onboard GNSS antenna and transducer also needed to be solved from the acoustic ranging. This led to an increase in the number of unknown parameters, which resulted in a reduction of the accuracy of the transponder array positioning. A further disadvantage of our data is the uneven ship track coverage over the transponder array. Due to the noise from the engine and crew, we had to cease boat operations and allow the boat to drift during the acoustic measurements. Sato et al. (2013) developed a hull-mounted onboard system, in which the transducer is installed on the bottom of the vessel, and this enables the vessel to sail along predetermined track lines. This improved the positioning repeatability and decreased the observation time for a single site from 2 to 4 days to 16–24 h (Ishikawa et al., 2020). Although our method had some disadvantages compared with that of the JCG group, we compensated for these disadvantages with long observation periods.

Average Rate of Back-Arc Rifting

Although the transponder array motion appears to have some temporal variations in both the horizontal and vertical components (Figure 3), we first consider the average motion (i.e., linear fitting without offsets) during 2012–2019. As proposed by Chen et al. (2018), if we assume the rigid Yonaguni Block is the western end of the Ryukyu arc, the Euler pole calculated from the 0499 (Yonaguni Island) and OILN sites is located off the eastern coast of Taiwan (Figure 4). The velocity predicted from the rotation of the Yonaguni Block at onland GNSS sites in the northern part of the Central Range is consistent with the observed GNSS velocities within ±15 mm/yr (Figure 4B). This suggests that clockwise rotation of the Ryukyu arc is the main driving mechanism for deformation of the Yilan Plain, which is consistent with evidence from earthquake focal mechanisms (e.g., Huang et al., 2012). Most of the earthquakes shallower than 25 km beneath the Yilan Plain and the Okinawa Trough show normal faulting mechanisms (Figure 4B). In addition, the GNSS velocities in northern Taiwan (Figure 4A), with respect to the stable Yangtze Plate, also exhibit a small eastward velocity of up to 10 mm/yr, associated with the clockwise rotation of the Yonaguni Block (Figure 4C). The velocity generally increases toward the trough axis. This E–W extension suggests that this area is also affected by back arc spreading. Schellart and Moresi (2013) used a 3D lithospheric deformation and mantle flow simulation to show that slab-rollback-induced toroidal and poloidal mantle flows are responsible for back-arc spreading. The overriding plate is extended by the basal shear of the asthenosphere, which flows into the space vacated by the rollback of the sinking slab. The E–W extension off the northern coast of Taiwan behind the spreading center may reflect the mantle flow beneath this region (Figure 4D), which has been predicted from a numerical mantle flow simulation based on the specific slab geometry (Lin and Kuo, 2016).

FIGURE 4

Chen et al. (2018), but also showing toroidal mantle flow.

Large Subsidence Rate

Chen et al. (2018) did not discuss the vertical motion of the OILN site because of its large uncertainty. However, we estimated the vertical motion more accurately, with the improved onboard configuration of the main GNSS antenna on top of the transducer pole. The total subsidence during 2012–2019 was ∼150 mm (Figure 3C). The average subsidence rate was 22.0 ± 8.9 mm/yr.

One possible interpretation of the high subsidence rate of the OILN site is sinking of the seafloor transponders into thick sediments. Taking into account that the transponder is ∼70 cm high and comprises a glass sphere that is 43 cm in diameter (Figure 2D), the transponders should have sunk fastest just after installation. However, after the installation of the new transducers in 2012, there was no acceleration in the subsidence. Another possibility for the high subsidence rate is a submarine landslide. The slope of the seafloor around this site is 1.5° to the north. When a landslide occurs, we would expect a large horizontal northward displacement compared with the subsidence. However, no such northward displacement was observed in the study period.

The other possibility is that subsidence is occurring due to crustal thinning at the margin of the back-arc basin, which is supported by the similar subsidence rate observed at the onland site on Gueishang Island (GS20; Figures 3C, 5A). In addition, onland sites in the Yilan Plain also show subsidence with rates of 5–10 mm/yr (Figure 5A). The subsidence rate of the OILN site is much faster than in active rift basins. For example, Steckler et al. (1988) showed from borehole profile data that the fastest tectonic subsidence rate in the central–southern Gulf of Suez was a few mm/yr in the late Burdigalian. This discrepancy may be explained by the unusual geological setting around the OILN site and Yilan Plain. The crust beneath the Okinawa Trough around the Yaeyama rift is continental, and the trough is in a transitional stage from continental rifting to seafloor spreading (Arai et al., 2017). If crustal thinning had a major role in the back-arc rifting around the OILN site, the thinning rate would be controlled by crustal thickness, spreading rate, and width of the trough. Assuming there is no change in the volume of the crust due to extension, the estimated thinning rate for a crustal thickness of 15–25 km (Arai et al., 2017) is 25–42 mm/yr, based on an average spreading rate of 50 mm/yr over a 30-km-wide trough. If isostasy is being attained, then thinning of the continental crust should be compensated by upward migration of the Moho discontinuity and seafloor subsidence. The ratio of the subsidence to thinning rate is governed by the density ratio between seawater, crust, and mantle according to the following equation:where Δh_c is the crustal thinning rate and Δd is the seafloor subsidence rate. ρ_w, ρ_c, and ρ_m are the densities of seawater, continental crust, and mantle, respectively. In this case, when typical density values of 1.0, 2.7, and 3.3 g/cm³ are adopted for ρ_w, ρ_c, and ρ_m, the expected subsidence rate Δd is 6.5–11.0 mm/yr for an estimated thinning rate Δh_c of 25–42 mm/yr.

FIGURE 5

The observed subsidence rate of 22 mm/yr at the OILN site is higher than expected from isostasy. In fact, there is a narrow, negative, free-air gravity anomaly region that extends from the Yilan Plain to the OILN site along the trough axis (WGM2012 Earth Gravity Model; Bonvalot et al., 2012), which does not correlate with the bathymetry. This suggests that gravitational mantle flow cannot compensate for the rapid deformation within the narrow region around the western end of the rift axis. Given the small negative free-air gravity anomaly to the east of the OILN site in the Okinawa Trough, this rapid subsidence should be limited to this westernmost area in the initial stages of back-arc rifting.

Temporal Variation in Horizontal Motion

We estimated an average velocity in Average Rate of Back-Arc Rifting, but Chen et al. (2018) reported that the transponder array position of the OILN site from 2012 to 2016 can be explained by intermittent motion, based on Akaike Information Criteria (AIC). We examined if intermittent rifting of the Okinawa Trough axis as occurred using a longer observation period from 2009 to 2019.

Chen et al. (2018) modeled the horizontal transponder array position from 2012 to 2016 with a linear fit and an offset in 2013. They suggested that the offset was associated with an acceleration of seismic activity around the OILN site in 2013. Figure 6A shows an enlargement of the horizontal transponder array motion from Figure 3, excluding the campaigns with limited acoustic data (<3,000 shots). When including the new campaigns since 2016, no offset is recognized between the sixth (July 2013) and seventh (September 2013) campaigns.

FIGURE 6

Chen et al. (2018), but were treated with our new procedure. (B) Cumulative seismicity (>M2.0) detected within the corresponding colored rectangles shown in (C,D). The red dashed line is a fit to the seismicity from January 2009 to July 2013. (C) Horizontal velocity vectors and seismicity (>M3.0) from Jan 1, 2009 to July 2, 2013. D1 and D2 show the location of dike intrusion modeled by Nakamura and Kinjo (2018), which occurred at the end of 2002 and in April 2013, respectively. (D) Horizontal velocity vectors and seismicity from July 3, 2013 to September 30, 2019.

However, the site appears to have moved at different rates before and after 2013. We estimated the time of the rate change and the offset due to transponder array replacement by least squares regression as the follows:where x is the horizontal position vector at time t, V₁ and V₂ are the horizontal velocities before and after the bending time t₀, respectively; and x₀ and x₁ are the intercepts before and after the transponder array replacement, respectively. The appropriate t₀ value was determined based on AIC. Note that we did not include the vertical component in this regression because uncertainty of the vertical component is much larger compared to the horizontal components.

The best-fit result is shown in Figure 6A, and yielded a t₀ of July 2, 2013. The minimum AIC value is −110.5. The offsets caused by the replacement of the transponder array are estimated as 5.6 and 12.7 cm for the south and west components, respectively. When fitting by a line instead of bending lines, the AIC value is −103.9. The relative likelihood of the linear fitting with respect to the bending-line fitting is calculated as exp{(103.9–110.5)/2} = 0.037. Consequently, it is plausible that the velocity of OILN changed in July 2013. The estimated t₀ value based on the horizontal motion of the OILN site is consistent with the timing of the change in seismicity. The seismicity rate around the OILN site increased abruptly at the end of 2013 (Figure 6B), which was at the time of dike intrusion in the Yaeyama rift (D2 in Figure 6C) reported by Nakamura and Kinjo (2018). Nakamura and Kinjo (2018) analyzed onland GNSS coordinate data and seismicity around Iriomote Island and showed that the shear strain rate changed around the time of dike intrusion, which was also associated with a seismicity change in the forearc region of Iriomote Island. Our results show that these changes also coincided with the seismicity and crustal deformation changes 150 km away from the Yaeyama rift. However, the motion of the OILN site does not show an acceleration in rifting, but instead exhibits a deceleration of rotation associated with an increase of seismicity, which appears to be counterintuitive. One possible explanation for this is activation of a normal fault to the south of the OILN site, which separates the Okinawa Trough from the Ryukyu arc. Arai et al. (2017) described a high-angle normal fault that separates the trough from the arc ∼100 km east of the OILN site, which appears to continue to the south of the OILN site, based on the bathymetry. The increased seismicity may reflect not only the opening of the Okinawa Trough, but also normal faulting south of the OILN site.

Chadwell et al. (1999) conducted a 2-years-long acoustic ranging experiment over a distance of 1 km across the Juan de Fuca ridge axis, and found that the ridge axis was not rifting. Chadwell and Spiess (2008) subsequently studied the plate motion using GNSS/acoustic methods at a location 25 km away from the same ridge axis, and found a motion indicative of constant seafloor spreading. They suggested that the crust between 0.5–25 km from the ridge axis accommodates the transient deformation by normal faults on the ridge slope and elastic strain accumulation. Our site is located 10 km away from the trough axis, and also shows a transient velocity change. Although we did not observe intermittent rifting, the velocity changes almost perpendicular to the rifting direction. The timing of the near-axis velocity change coincides with a change in the strain rate (Nakamura and Kinjo, 2018) and the inter-plate slow slip event (SSE) occurrence rate (Tu and Heki, 2017). The unique location of the OILN site close to the rifting center allows changes in the deformation field to be detected over time, which is not possible with land-based stations.

Differences From the Study of Chen et al. (2018)

Figure 7 shows the horizontal positions determined by Chen et al. (2018) superimposed on the result of the present study. There are significant differences of up to 10 cm in the N–S positions of campaigns conducted in 2012 and 2014. In 2012, the position accuracy of Chen et al. (2018) is degraded because of the limited acoustic data (<3,000). In this study, we merged the data from the campaigns in April 2012 and July 2012 that were previously discarded by Chen et al. (2018) due to lack of data, in order to increase the available data to 5,828. Despite the absence of such differences in the 2014 acoustic data, there are significant differences in the estimated benchmark position. The difference between the results of Chen et al. (2018) and this study for the 2014 campaign reflect the KGPS analysis. Chen et al. (2018) solved for the GNSS antenna position using Interferometric Translocation software (IT; Colombo, 1998), whereas we solved for the same antenna using the bow antenna as the fiducial site with GrafNav software (NovAtel Inc, 2014). As a result, the positions of the onboard transducer estimated in this study differ from those in Chen et al. (2018). On average, the positions determined in this study differ from Chen et al. (2018) by 0.7 ± 5.0 cm in a westward direction, 3.8 ± 3.3 cm in a northward direction (Figure 8A), and 15.6 ± 15.7 cm in a downward direction (Figure 8B). We compared the elevation with the theoretical sea level in order to evaluate which result is more reliable (Figures 8C,D). The results show that the elevations obtained in this study are more consistent with sea level than those of Chen et al. (2018). This trend is also evident in some other campaigns, but is especially large in April 2014. It is possible that the coordinates of the seafloor position might have been influenced by this unexpected bias in the kinematic GNSS positioning of Chen et al. (2018).

FIGURE 7

Chen et al. (2018). (B) E–W component.

FIGURE 8

Chen et al. (2018). (A) Relative horizontal positions from Chen et al. (2018). The red cross shows the average position. (B) E–W cross-section. (C) Time-series of the vertical transducer position. Positions from Chen et al. (2018) and this study are shown as purple and light green dots. One-hour running averages are shown as green and red dots. The predicted astronomical tide is shown as black dots. (D) Height difference between the transducer position and the astronomical tide level.

In the time-series obtained by low-frequency campaigns, the temporal changes are heavily biased by the results of only a few campaigns (e.g., 2012 and 2014 in Chen et al. 2018). This is also the case for our study. It should be noted that the estimated motion before 2012 obtained in this study was determined mainly by the location of the campaign in 2011. However, irrespective of the pre-2011 coordinate values, the eastward velocity appears to have decreased over time, at least when comparing 2012–2016 with 2016–2019 (Figure 6).

Conclusion

We measured seafloor movement near the rift axis of the Okinawa Trough back-arc basin over a 10-years-long period using a GNSS/acoustic technique. We deployed a seafloor site that is 60 km east of northeastern Taiwan in 2009, and conducted 18 measurement campaigns until 2020. Chen et al. (2018) reported the results for campaigns from 2012 to 2016 and found possible temporal variations in the rifting velocity of the Okinawa Trough. In this study, we extended the measurement period to 2019, in order to obtain more reliable estimates of velocity. The seafloor transponder array has undergone southeastward movement at an average velocity of 43 mm/yr from 2012 to 2019, which is 18 mm smaller than the results of Chen et al. (2018) for 2012–2016. This horizontal movement can be explained by clockwise block rotation of the Yonaguni Block and Central Range, which form the Yilan Plain, by back-arc rifting. In addition to the horizontal motion, the transponder array revealed fast subsidence in the 7 years from 2012 to 2019, with an average velocity of 22 mm/yr. The rapid subsidence and negative gravity anomaly suggest that rapid crustal thinning is not being fully compensated for by upward migration of the Moho, but is being accommodated by surface deformation. We also attempted to retrieve the acoustic ranging data obtained from 2009 to 2012, which were initially discarded because of replacement of the seafloor transponders. By fitting the decadal horizontal motion with two linear lines with a break and an offset due to the transponder replacement, the break is estimated to have occurred in July 2013. This timing coincides with a change in the seismicity rate around the site, as well as dike intrusion 150 km away from the site.

References

• 1

  AndoM.IkutaR.TuY.ChenH. Y.LinC. H. (2015). “The Apr 2013 earthquake swarm and dyke intrusion in the Okinawa trough, abstract for, 2015,” in Japan Geoscience Union Meeting, Makuhari, chiba, May 24–28, 2015, SCG62–SCG07.

  • Google Scholar
  • View reference in article

• 2

  AraiR.KodairaS.KaihoY.TakahashiT.MiuraS.KanedaY. (2017). Crustal structure of the southern Okinawa Trough: symmetrical rifting, submarine volcano, and potential mantle accretion in the continental back-arc basin. J. Geophys. Res. Solid Earth.122 (1), 622–641. 10.1002/2016jb013448

  • CrossRef
  • Google Scholar
  • View reference in article

• 3

  BarnhillR. R.RiesenfeldR. F. (Editors) (1974). Computer aided geometric design. New York: Academic, 326.

  • Google Scholar
  • View reference in article

• 4

  BonvalotS.BalminoG.BriaisA.KuhnM.PeyrefitteA.ValesN.et al (2012). World gravity map. Commission for the geological map of the world. Paris: BGI-CGMW-CNES-IRD.

  • Google Scholar
  • View reference in article

• 5

  ChadwellC. D.HildebrandJ. A.SpiessF. N.MortonJ. L.NormarkW. R.ReissC. A. (1999). No spreading across the southern Juan de Fuca ridge axial cleft during 1994-1996. Geophys. Res. Lett.26 (16), 2525–2528. 10.1029/1999gl900570

  • CrossRef
  • Google Scholar
  • View reference in article

• 6

  ChadwellC. D.SpiessF. N. (2008). Plate motion at the ridge-transform boundary of the south Cleft segment of the Juan de Fuca Ridge from GPS-Acoustic data. J. Geophys. Res.113 (B4), B04415. 10.1029/2007JB004936

  • CrossRef
  • Google Scholar
  • View reference in article

• 7

  ChenH. Y.IkutaR.LinC.-H.HsuY.-J.KohmiT.WangC.-C.et al (2018). Back-arc opening in the western end of the Okinawa Trough revealed from GNSS/Acoustic Measurements. Geophys. Res. Lett.45 (1), 137–145. 10.1002/2017GL075724

  • CrossRef
  • Google Scholar
  • View reference in article

• 8

  ColomboO. L. (1998). “Earth Sciences and Geography,” in GPS for geodesy. 2nd Edn, Editors KleusbergA.TeunissenP.PeterJ. G. (Berlin: Springer), 537–567.

  • Google Scholar
  • View reference in article

• 9

  DeMetsC.GordonR.ArgusD. (2010). Geologically current plate motions. Geophys. J. Int.181 (1), 1–80. 10.1111/j.1365-246X.2009.04491.x

  • CrossRef
  • Google Scholar
  • View reference in article

• 10

  HuangH. H.ShyuJ. B. H.WuY. M.ChangC. H.ChenY. G. (2012). Seismotectonics of northeastern Taiwan: kinematics of the transition from waning collision to subduction and postcollisional extension. J. Geophys. Res.117 (B1), B01313. 10.1029/2011JB008852

  • CrossRef
  • Google Scholar
  • View reference in article

• 11

  IkutaR.TadokoroK.AndoM.OkudaT.SugimotoS.TakataniK.et al (2008). A new GPS-acoustic method for measuring ocean floor crustal deformation: application to the Nankai Trough. J. Geophys. Res.113 (B2), B02401. 10.1029/2006JB004875

  • CrossRef
  • Google Scholar
  • View reference in article

• 12

  IshikawaT.YokotaY.WatanabeS.NakamuraY. (2020). History of on-board equipment improvement for GNSS-A observation with focus on observation frequency. Front. Earth Sci.8, 150. 10.3389/feart.2020.00150

  • CrossRef
  • Google Scholar
  • View reference in article

• 13

  LaiK. Y.ChenY. G.WuY. M.AvouacJ. P.KuoY.-T.WangY.et al (2009). The 2005 Ilan earthquake doublet and seismic crisis in northeastern Taiwan: evidence for dyke intrusion associated with on-land propagation of the Okinawa Trough. Geophys. J. Int.179 (2), 678–686. 10.1111/j.1365-246X.2009.04307.x

  • CrossRef
  • Google Scholar
  • View reference in article

• 14

  LallemandS.LiuC. S. (1998). Geodynamic implications of present-day kinematics in the southern Ryukyus. J. Geol. Soc. China.41 (4), 551–564

  • Google Scholar
  • View reference in article

• 15

  LinS. C.KuoB. Y. (2016). Dynamics of the opposite-verging subduction zones in the Taiwan region: insights from numerical models. J. Geophys. Res. Solid Earth.121 (3), 2174–2192. 10.1002/2015JB012784

  • CrossRef
  • Google Scholar
  • View reference in article

• 16

  NakamuraM.KinjoA. (2018). Activated seismicity by strain rate change in the Yaeyama region, south Ryukyu. Earth Planets Space.70 (1), 154. 10.1186/s40623-018-0929-y

  • CrossRef
  • Google Scholar
  • View reference in article

• 17

  NishimuraS.HashimotoM.AndoM. (2004). A rigid block rotation model for the GPS derived velocity field along the Ryukyu arc. Phys. Earth Planet. In.142 (3-4), 185–203. 10.1016/j.pepi.2003.12.014

  • CrossRef
  • Google Scholar
  • View reference in article

• 18

  NovAtel Inc (2014). GrafNav 8.60: user guide OM-20000147 rev 3. Available at: https://hexagondownloads.blob.core.windows.net/public/Novatel/assets/Documents/Waypoint/Downloads/NavNet860_Manual/NavNet860_Manual.pdf (Accessed January 26, 2021).

  • Google Scholar
  • View reference in article

• 19

  SatoM.FujitaM.MatsumotoY.IshikawaT.SaitoH.MochizukiM.et al (2013). Interplate coupling off northeastern Japan before the 2011 Tohoku-oki earthquake, inferred from seaoor geodetic data. J. Geophys. Res. Solid Earth.118 (7), 3860–3869. 10.1002/jgrb.50275

  • CrossRef
  • Google Scholar
  • View reference in article

• 20

  SchellartW. P.MoresiL. (2013). A new driving mechanism for backarc extension and backarc shortening through slab sinking induced toroidal and poloidal mantle flow: results from dynamic subduction models with an overriding plate. J. Geophys. Res. Solid Earth.118 (6), 3221–3248. 10.1002/jgrb.50173

  • CrossRef
  • Google Scholar
  • View reference in article

• 21

  SibuetJ. C.LetouzeyJ.BarbierF.CharvetJ.FoucherJ. P.HildeT. W. C.et al (1987). Back arc extension in the Okinawa trough. J. Geophys. Res.92 (B13), 14041–14063. 10.1029/JB092iB13p14041

  • CrossRef
  • Google Scholar
  • View reference in article

• 22

  StecklerS. M.BerthelotF.LyberisN.Le PichonX. (1988). Subsidence in the gulf of suez: implications for rifting and plate kinematics. Tectonophysics.153 (1–4), 249–270. 10.1016/0040-1951(88)90019-4

  • CrossRef
  • Google Scholar

• 23

  TuY.HekiK. (2017). Decadal modulation of repeating slow slip event activity in the southwestern Ryukyu Arc possibly driven by rifting episodes at the Okinawa Trough. Geophys. Res. Lett.44 (18), 9308–9313. 10.1002/2017GL074455

  • CrossRef
  • Google Scholar
  • View reference in article

• 24

  YokotaY.IshikawaT.WatanabeS.TashiroT.AsadaA. (2016). Seafloor geodetic constraints on interplate coupling of the Nankai Trough megathrust zone. Nature.534 (7607), 374–377. 10.1038/nature17632

  • CrossRef
  • Google Scholar
  • View reference in article