Using GMT for 2D and 3D Modeling of the Ryukyu Trench Topography, Pacific Ocean

The Ryukyu Trench stretches along the southeastern edge of the Ryukyu Islands Arc between the Philippine Sea, Taiwan, the Japanese island of Kyushu and the Pacific Ocean with a length of 1,400 km (Ando et al. 2018) (Figure 1). The system of the Ryukyu Trench and Ryukyu Island Arc serves as a natural border between two seas with highly contrasting bathymetric settings: the East China Sea (ECS) and the Philippine Sea. The ECS differs from other marginal seas by its highly developed shelf, which occupies most of its seafloor and extends to the Yellow Sea. The shallow depths of the ECS clearly contrast with the depths of the Philippine Sea basin (Figure 1). The Ryukyu Island Arc is located on a submarine ridge with depth ranges of between 500 and 1000 m. The slopes of the ridge have a convex geomorphic profile and a stepped structure, formed by ancient folded areas bounded by Cenozoic and modern geosynclinal strata.

Figure 1

The Ryukyu Trench, located between the two seas with such contrasting bathymetry, extends along the foot of the outer slope of the Ryukyu Island Arc and reaches a maximal depth of 7,790 m (Litvin 1987).

The bathymetry of the Ryukyu Trench varies in its N and S segments, reflecting the correlation between its geomorphic shape, gradient steepness and local tectonic settings (Lemenkova 2019g). In its N part, the Philippine Sea Plate subduction is notable for shallow depths and low levels of steepness, reaching only 11° approximately. The Ryukyu Trench reaches depths of 7,000 m below sea level (BSL) (Figure 2). The N part of the Ryukyu Island Arc is divided into three parallel ridges. The forearc slope between the forearc island chain and the trench is between 130 and 160 km wide. The N group has Yaku-shima and Tanega-shima and the 200-km-long Tokara Channel gap between the island groups (Okamura et al., 2017). The depression of the Okinawa Trough is located north of the Ryukyu Island Arc (Figure 2). The continental crust decreases under the Okinawa Trough and increases again under the Ryukyu Island Arc. The basement surface of the Paleo-Mesozoic and Cenozoic geologic structures in the ECS is complicated by several sub-meridional depressions (with depths of >2 km) separated by uplifts (with depths of <1 km).

Figure 2

The Ryukyu Trench was formed as a result of the subduction of the oceanic crust of the Philippine Sea Plate beneath the continental crust of the Eurasian Plate (Doo et al., 2018). The tectonic settings of the Ryukyu Trench reflect the evolutional processes of the building seamounts and structure of the Earth's crust. Therefore, they caused variations in the lithological patterns, the geological conditions of the rock formation, the stratigraphy and sedimentary settings of the ocean seafloor during the Jurassic, Mesozoic and Cretaceous (Lemenkova 2019h). The double subduction system of the Ryukyu Trench reflects complex trench motions as the expression of the interaction between the subducting plates and the underlying mantle of the Ryukyu Trench. The Pacific Plate subducts beneath the Philippine Sea Plate, which in turn subducts beneath the Eurasian plate along the Ryukyu Trench (Faccenna et al. 2018). The subduction of the Philippine Sea Plate causes the eruption of approximately 34 volcanoes (Figure 2) and earthquakes that produce tsunami.

Therefore, the Ryukyu arc-trench system is marked by active submarine margins of the continents located near northern Taiwan in the southern Ryukyu forearc (Figure 2). Southwestern Ryukyu is characterized by the edge regime influenced by the juxtaposing continental lithosphere (Kuo et al. 2012). Variations in the geomorphology of the trench follow the geometry of the slab dip which differs regionally: gentler in central (C) and southern (S) segments (40–50°) at a 7 km depth and steeper (70°) at a >7 km depth in the northern ones (Kodaira et al. 1996).

The sediment types of the Ryukyu Trench include mud diapirs, hemipelagic sediments, mud-supported breccia, fossil planktonic foraminifera, and calcareous nannofossils (Ujiie 2000). The sediment composition of the Ryukyu Trench changes spatially along the Ryukyu arc from smaller in the N segment arc to heavier along the C and S parts. This is explained by the dominance of the terrigenous and pelagic sedimentation along the N and C Ryukyu arc (Shu et al. 2017). The deep-sea trench presents a kind of structural trap in the continental margins of the Pacific, as sediments are carried by the ocean currents passing through the system of trenches of the Pacific Ocean (Lemenkova 2018). At the same time, the sediment thickness and composition of the trench seafloor strongly depend on its location and curvature (slope steepness). According to the data obtained from deep-sea drilling (Litvin 1987), the sedimentation structure of the northwestern section of the Philippine Sea is composed of Eocene–Pleistocene low carbonate silts and clays, with volcanic tuffs, conglomerates, and basalts beneath them. The carbonate content of precipitation increases markedly on the geomorphic positive relief forms (small seamounts and hills).

Although the seafloor of the Ryukyu Trench is separated vertically from the euphotic zone, it serves as a depositional center for organic carbon accumulation. This is closely connected with geomorphic factors: steep slopes enable lateral transport from the continental margins (Itoh et al. 2011). As a consequence, the spatial biodiversity of the Ryukyu Trench has lower values at shallow depths, but these increase with water depth reaching maximum levels between the trench floor and the abyssal plain (Kitahashi et al. 2014). This proves that the distribution of the deep-sea biomass of the Ryukyu Trench spatially differs in various regions of the trench following geomorphic variations. Such discrepancies are a result of differences in environmental factors: primary productivity level. Moreover, migration patterns that have been discovered in the Ryukyu Trench area show that reefs and geomorphic landforms serve as land bridges for the migration of terrestrial organisms from Okinawa-jima to Miyako-jima (Arai et al. 2016).

The dataset used in the present study includes the GEBCO global terrain model downloaded from the web (GEBCO Committee 2016). The GEBCO project was originally initialized by HSH Prince Albert I of Monaco in 1903 to enable studies of the physical shape, geomorphology and bathymetric characteristics of the seafloor (Weatherall et al. 2015). Since the initiation of the GEBCO project at the onset of the 20th century, it has been significantly improved due to the technical and cartographical tools that have been developed, while the goal of mapping the bathymetry of the world's oceans has been maintained until recent times (Mayer 2006). The GEBCO project is focused on the high-precision mapping of the ocean seafloor and is developed and maintained by several oceanographic organizations (Weatherall et al. 2015). The GEBCO_2019 grid used in the current research is methodologically based on the development of the previous versions of GEBCO (Schmitt & Weatherall 2014) and the SRTM15_PLUS dataset (Olson et al. 2014). Besides, GEBCO is collaborating with regional groups to continually improve its global model: IBCAO and IBCSO gridded datasets (Jones 1994; Schmitt & Weatherall 2013; IOC, IHO & BODC 2003; GEBCO_2010). An ETOPO1 15 arc-sec grid was downloaded from the publicly available NOAA resources (Amante & Eakins 2009).

Geoid and marine free-air gravity field grids were based on Geosat and ERS-1 altimetry raster grids (Smith & Sandwell 1995; Sandwell et al. 2014). The geographic names were verified in the GEBCO Gazetteer of submarine features (IHO-IOC 2012). The vector lines on the maps (tectonic slabs, trench lines, coastlines, tectonic faults) were taken from the embedded layers of GMT (Wessel & Smith 1996). The data were converted to several projections: Mercator (Figure 1, Figure 6, C), Lambert Conic Conformal (Figure 2), Equidistant Conic (Figure 3), Cylindrical Equal-Area Gall-Peters (Figure 4), and a 3D model (Figure 5). A subset covering the area with geographic coordinates of 120°00′E 134°00′E 20°00′N 33°00′N was used for the study.

Figure 3

Figure 4

Figure 5

GMT is an open-source set of over 50 command-line cartographic commands suitable for manipulating big data both in 2D and 3D formats (arrays of XY and XYZ datasets in ASCII format). Shell scripts are executable on UNIX, Linux or Mac OS X. The GMT-integrated scripts provide a robust solution for generating high-accuracy maps based on GEBCO and ETOPO1 grids, extracting isolines, coastlines, relief, and plotting publication-quality maps. Eight GMT commands - integrated shell scripts - were written for plotting maps and graphics. Several GMT modules were used for cartographic visualizations.

For instance, ‘grdcontour’ was used to extract isolines from a GEBCO grid; ‘psconvert’, ‘psscale’ to visualize scale bars; ‘pstext’ to add annotations; ‘psxy’ to plot geologic lineaments and sample points; ‘psbasemap’ to map cartographic elements (isolines, geoid, bathymetry, net grids and titles), and the module ‘grdtrack’ to automatically extract data for plotting cross-section profiles from a GEBCO grid.

Each GMT module consists of a small line of code. Combined as a script, they were used to produce maps. The ‘grdimage’ module was used to map Figures 1, 3 and 4 (bathymetric grid, geoid, and gravity, respectively) using a sequence of the GMT codes. The 3D modelling was done using the ‘grdview’ module with the command: grdview rt_relief.nc -J -R -JZ2i -CrainbowRT. cpt -p115/30 -Qsm -N-7500+glightgray -Wm0.07p -Wf0.1p,red -Wc0.1p,magenta -B4/4/2000:”Bathymetry and topography (m)”:ESwZ -S5 -Y6.0c -O -K >> $ps. The 3D map was overlaid with the 2D geoid model (Figure 5) presenting a combined 2D–3D modelling overlay approach. The trend curves were modelled using the ‘trend1d’ module. Statistical modelling and digitizing cross-section profiles were done using the ‘pshistogram’ and ‘grdtrack’ modules.

The design of GMT follows the UNIX approach in which multi-tasks are broken down into small and independent ‘self-contained’ commands and small sub-tasks. Consequently, taken together as a chain of such sub-tasks, individual GMT module commands were combined into a shell script for batch processing. A workflow of the GMT shell scripts contains a sequence of modules; each of them was used for a specific task. Visualizing raster grids as maps was performed using the ‘grdimage’ GMT module with local arguments set up for plotting grids in specific map projections and coordinate systems. For example, a command -R120/134/20/33 -JM6i means selecting a square region with 120°E 134°E 20°N 33°N coordinates and plotting a map in Mercator projection with a polygraphic (printing) dimension of the map layout of 6 inches. Applying color palettes was done using GMT module ‘makecpt’, cleaning data was performed using ‘rm’ from the UNIX utility. A map legend was plotted using the ‘pslegend’ module with annotating texts using the ‘pstext’ module, generating a study area through ‘cutting’ selected squares (‘grdcut’). Delineating isolines and adding bathymetric contour lines used to represent elevation from DEM was performed using the ‘grdcontour’ module.

3D mapping allowed the profiling of the topographic objects (cross-sections) in three dimensions. The map was technically based on a ‘grdview’ module which uses an advanced cartographic method for surface visualization and gathering information from a 3D scene. Extracting cross-section profiles was done by ‘grdtrack’ which enables linear graphs of the topographic curve from the start and end points of the cross-sectional segment to be plotted. The curve of the selected segment on a profile shows elevations perpendicular to the surface. Cross-section profiles show the gradient, slope curvature and bathymetry represented on a 2D model. The 3D modelling of the plot with a mesh surface was done using a ‘grdview’ module. For 3D modelling, bathymetric data were exported as a NetCDF (Network Common Data Form) file with an .nc extension (@rt_relief.nc) and imported into GMT for raster processing, 3D visualization and evaluation. Bathymetric data were gridded at a 1 arc-minute resolution from the original ETOPO1 file (earth_relief_01m.grd).

Error bars represent uncertainty in bathymetric data and the statistical range of data distribution. Bathymetric uncertainty can be also analyzed using QGIS plugins, a combination of R and Python, to implement methods and tools aimed at generating an estimate of the uncertainty (Schmitt et al., 2015). The estimation of uncertainty depends on the methodology (in particular, on its accuracy) of the interpolation algorithms adapted for bathymetric use. For instance, interpolation algorithms of topographic isolines include inverse distance weighted interpolation (IDW), triangulated irregular network (TIN), B-Spline and Spline under tension (Charpy et al., 2015). Error bars (thin gray vertical lines on Figure 6, A and B) accompanying median values of the cross-section profile (red line on Figure 6, A and B) graphically present the variability of the bathymetric data used to indicate error or uncertainty in the input raster grid. A collection of GMT modules was used to perform mapping, projecting, rasterization, plotting isolines using techniques of interpolation and smoothing of raster grids (Wessel & Smith 2015).

Figure 6

The GMT ‘psconvert’ command was applied to convert PS files to JPEG (and can also be used for other desired extensions, such as EPS, PDF, and TIFF) through the following command: gmt psconvert RTtrends2.ps -A20.0c -E720 -Tj -P -Z, where ‘RTtrends2.ps’ is the name of the file, and -E720 is a 720 dpi resolution file in a JPG format (-Tj). Adding essential elements on the base map (title, legend, directional rose, grid, scale, time stamp) was technically possible using ‘psbasemap’; calculating statistical mean and median, using ‘pshistogram’. Error bars showing the standard deviation of the elevations were plotted using the ‘gmt grdtrack trenchRTn.txt -Grtrelief.nc -C400k/2k/20k -Sm+sstackRTn.txt > tableRTn.txt’ command.

Each of the eight maps presented was created using eight GMT scripts, derived from Lemenkova (2019i) and modified to carry out this research. Within an individual GMT shell script, the attributes of each of the maps (coordinates of the map square, color palette, annotations, border ticks, legend and other elements) were strung together as a series of commands. Codes for text annotations were used as a template for iterations to repeat texts with modifications (coordinates, font size, color, rotation), as can be seen in Figure 2. The serious polygraphic advantage of GMT is its operation with output files. It creates print-quality outputs and graphics in a PostScript format, an industry standard for high-quality printing, which can be converted to other graphical formats (JPG, TIFF, BMP).

The geomorphology of the Ryukyu trench-arc seafloor results from the complex geodynamic settings and tectonic processes affecting its structure. These demonstrate traces of the continuous evolution and shape of the seabed surface as a system of mid-ocean ridges, trenches, troughs, transform faults, vast ocean plains complicated by ridges, plateaus, and such minor geomorphic forms as terraces, seamounts and hills (Figure 5). Perspective 3D scene visualization focused on the evaluation of such parameters as submarine topography and the comparison of the bathymetric ranges in different parts of the trench. Various factors affect variations in the submarine topography of the trench. Along with regional patterns, local variations are affected by the geological settings which include the character of the sedimentation processes. For example, the seafloor of the ECS is mainly covered by a thick continental crust with a thickness of 30 km approximately.

A raster layer was generated as a color-coded surface at the same resolution as for the visualized 3D model of the topographic slope, overlaying the geoid contour plot (Figure 5). Azimuthal rotation from a 115°/30° direction was added by rotating the initial 3D grid surface in an azimuth using the ‘grdview’ module. Although GEBCO presents the highest quality for the global terrain model for ocean and land (15 arc-second intervals), the ETOPO1 raster grid was selected as a base map for 3D modelling, since it has enough density to grid at an acceptable resolution (1 min). The ETOPO1 reduced the modelling and computation time for the 3D mesh and the limitations for data processing and visualization. The seafloor relief correlates with the core basement: the surface of the Ryukyu island arcs with steep ledges descends to the bottom of the Philippine Sea basins and Ryukyu Trench. Depths at the bottom of the seafloor generally exceed 5 km, in some places up to 6 km, and reach 7,790 m in the trench. Spatial differences of the Ryukyu Trench topography are mirrored in its sediment composition that changes systematically from lighter in the northern Ryukyu arc to heavier along the central part of the arc.

The comparison between two spatially different segments of the trench show that the N segment has higher depth values compared with S. This can also be noted in Figure 1 where the deepest values of the central part are clearly seen in the area extending seaward from the Miyako Fault (Figure 2). As can be seen when comparing Figures 6 A and B, the median stacked profile on the S part (Figure 6 A) demonstrates a slightly gentler slope shape with submarine terraces, which is not clear on the N part of the forearc slope (Figure 6 B). At the same time, the depth of the forearc slope increases in the N segment of the trench (Figure 6 B).

The geometric approximation of the trend curves (Figure 7) shows similarity between the segments in the applied coarse linear model and differences in geometric variations at applied polynomial functions.

Figure 7

Statistical data processing (Figure 8) shows differences in the depth ranges of the Ryukyu Trench along the island arc. Depth and geometry of the slope vary considerably: maximal values of −5,500 to −5,700 m are notable for the N segment where 334 sample points were recorded (Figure 8 A).

Figure 8

The forearc basin at a depth of >2,000 m and <1,000 m BSL shows more samples in the N segment than in S (Figure 7). In contrast, a wide terrace with depths that are shallower than 25 m is generally to be found more frequently (161 samples) in the S part of the trench.

Variations of the slope geomorphology are highlighted by the contour of the increased marine free-air gravity anomaly, as a triangle shaped ‘island’ located at 126°E 24°N (Figure 4). Visualized gravity anomalies (Figure 4) show local variations of the gravity field around the Ryukyu arc-trench system and region of the East China and Philippine seas. Positive anomalies in the range of 40 to −240 mGal, (orange to red areas on Figure 4) demonstrate more gravity than predicted by the model which is explained by the presence of a sub-surface positive mass anomaly in the Ryukyu Island Arc and Taiwan Island. On the contrary, negative anomalies in the range of −20 to −240 mGal (colors blue to purple on Figure 4), clearly depicting the stripe of the Ryukyu Trench area showing lower values than predicted, which indicates a sub-surface mass deficit in the area.

In general, transition zones in the region from the ECS eastwards along the Philippine Sea towards the Mariana Islands are characterized by increased free-air gravity value (Faye's) anomalies of 0–25 mGal (green and light blue-greenish colors in Figure 4 which cover the majority of the area). A clear exception can be seen in the Ryukyu Trench: a purple-colored stripe with values below −42 mGal (Figure 4). Furthermore, Faye's anomalies on the shelf of the ECS are in the range of 20–30 mGal; in the Okinawa basin, they increase to 80–100; over the island arc of Ryukyu, they decrease again to 20–40 mGal; and in the Philippine Sea basin, they reach 40–80 mGal (Figure 4). These results can be compared with previously measured Bouguer anomalies which are slightly positive on the shelf of the ECS, reach 160 mGal in the Okinawa Trough and 60–80 mGal over the Ryukyu Island Arc (Gaynanov 1980).

Being a joint structure, the Ryukyu Trench is located on the slopes of the depression of the surface of the upper mantle, at the junction of the subcontinental and oceanic crust types. The thickness of the Earth's crust in the Ryukyu Trench is 8.2 km (in comparison, it is 12.7 km in the Kuril-Kamchatka Trench). The thickness ratio of the sedimentary and basalt layers varies from 1:1 to 1:3 (Milashin 1971). The thickness of the basalt layer of the trench's seafloor exceeds the thickness of the same layer in the oceanic crust and is less than that of the subcontinental crust of anticline ridges. The mosaic distribution of crust types in such geosynclinal systems as oceanic trenches corresponds to the extreme values in gravity anomalies. Positive gravity anomalies (Bouguer reduction) are generally characterized by geosynclinal basins, which corresponds with the results presented. The belt of clearly expressed negative gravity anomalies along the Ryukyu Trench demonstrated corresponds with the general pattern of free-air gravity in deep-sea trenches (a rare exception is the Cayman Trench). The geoid model shows geoidal undulation in the SW part of the area east of Taiwan (Figure 3). Generally, a major part of the area has geoid values ranging from 25 to 35 mGal. Increased values (>35 to 85 mGal) can be seen on the SE part in the Philippine Sea, while lower values (15–25 mGal) bordering the continental shelf stretch along the ECS, and the lowest (<15 m Gal) mainly occupy land areas.