Cartographic Interpretation of the Seafloor Geomorphology Using GMT: a Case Study of the Manila Trench, South China Sea

The study is geographically focused on the Manila Trench, located in the west Pacific Ocean, South China Sea, west Philippines. The research aims at the geological mapping, analysis and visualizing variations in the submarine geomorphology of the Manila Trench. Technically, the work was done using Generic Mapping Tools scripting toolset (GMT). A combination of various GMT modules was applied for geospatial modelling. Methodology includes cartographic data integration and interpretation through approaches of data analysis: topographic plotting, geophysical modelling, geological mapping and statistical analysis. The data included SRTM, ETOPO1, geoid and gravity grids (CryoSat-2, Jason-1). Two sets of the cross-section profiles of the trench were automatically digitized. The profile transects were compared and differences in the geomorphic shape in southern and northern parts revealed. Southern part has steeper slope on the western part. Northern part is steeper on the continental slope part. The submarine terraces are located on the northern segment at depths -2,000 m. The depth and geomorphology of the slope vary for the range -3,500 to -4,500 m: minimals for the northern part with 526 samples (18.2%) for the depths -4,000 to -4,200 m. The histogram for the northern part has bimodal distribution with two peaks. The southern part shows 142 values for the minimals -3,500 to -3400 m. The statistical analyses revealed that northern part of the trench is deeper. The GMT functionality shown in this paper enabled integration and interpretation of the multi-source data: automatically digitized profiles, geological mapping, 2D and 3D bathymetric modelling, statistical analysis, mathematical approximation of the trend modelling. The GMT proved to be capable of visualizing geodata that can significantly improve Earth studies and interpretation of submarine geomorphology of the oceanic trenches through the advanced data analysis.


INTRODUCTION
Understanding a variety of the submarine landforms hidden from the human eyes due to their remote location, and factors affecting seafloor geomorphology, is a very complex task that requires a multidisciplinary approach: geological data analysis, data processing and modelling by advanced algorithms, geostatistical analysis and visualization. The asymmetry of the deep-sea trenches reflects a phenomenon of tectonic plates subduction. Thus, as one plate bends down to the Earth's mantle, another plate is being deformed filling the growing empty space. The depths of trenches are influenced by many processes and factors controlling their actual shape and structure. Many attempts were undertaken to map and visualize general form of the submarine geomorphology and model oceanic deep-sea trenches. Rapid development of the IT technologies in XXI century facilitated geodata processing which contributes to our better understanding of the marine geomorphology.
The significance of the research presented in this paper is associated with the presented approach of GMT as an advanced toolset among geoinformation GIS technologies applied for geological mapping and geomorphic analysis. Having access to the machine learning technologies associated with rapid recent IT progress (scripting coding, geospatial datasets and DEMs), data analysis offered by GMT can be employed for marine geological research to provide new, more detailed insight into the seafloor bathymetry. Presented GMT codes enable to perform a speed yet quality mapmaking and data processing. In turn, this allows to study in details submarine geomorphic landforms of the oceanic seafloor. In turn, this allows an in-depth study of the Manila Trench.
The contribution of this research towards development of the methodologies of seafloor mapping and modelling geodata by GMT consists in the presented and explained GMT scripting codes and snippet.
Rather than traditional GIS interface, a GMT is fully based on the scripting approach where cartographic methodology consists in executing the programming code which results in mapping. Geoinformation is crucial for understanding ocean seafloor and precise bathymetric mapping. Combining datasets of the raster maps with geomorphic modelling increases our knowledge of the factors affecting seafloor landforms.
The research aim is comparative geomorphic analysis of deep-sea oceanic trench: Manila Trench located in the geologically complex region of the western part of the Pacific Ocean, a part of the 'Ring of Fire' where active movements of the tectonic plates and volcanism take place. Complex geophysical settings affect the formation of the trench, high seismicity and geodynamic instability visualized on the thematic maps of geological and tectonic settings as the most important causes of the ocean trench formation.
Cartographic objective of the research was to visualize bathymetry, geomorphology, tectonic and geological settings of the trenches through 2D and 3D modelling and mapping. The statistical analysis aimed to compare datasets of the trench and show differences in geomorphic and bathymetric structure.

STUDY AREA AND MATERIALS
The study geographically focuses on the Manila Trench, an oceanic trench located in the South China Sea, west Pacific Ocean, west of the islands of Luzon and Mindoro in the Philippines (Fig. 1).
The study presents geomorphic mapping and the use of data models of various spatial resolution to visualize certain geophysical behaviors of the crust in the region. The area of the Manila Trench is located in the Manila subduction zone at the Philippine Sea plate boundary where it moves in the northwest direction toward the Eurasian plate with a high convergence rate [1].  The hypocenters of the tsunami of the Manila Trench are located at the depths <100 km [4].
The tsunami hazard from the Manila Trench source has been assessed in more details in several research papers [5][6][7]. The zone of the Eurasian Plate subduction explains the belt of volcanoes in the Manila Trench area (Fig. 2), on the west side of the Philippine island of Luzon. The area between the northernmost Manila subduction zone and southern Taiwan is considered a transition from the subduction to initial collision, and a weakly coupled nature for the northern part of the Manila Trench [8]. Because the structure and submarine geomorphology of the ocean trench is affected by Earth's interior through tectonic plates subduction and movements, the visual shape of the geoid is distinctly repeating the isolines, as shown on image above. Image source: author, made using GMT.
Submarine regions, in contrast with the terrestrial ones, still remain the least accessible areas on the Earth due to their remote location. Therefore, modelling oceanic trenches is one of the most complicated issues in the marine geology explained by their inaccessibility. Oceanic trenches can only be studied using computer based modelling, advanced mapping and algorithms of the data analysis.
Therefore, using special software and tools designed for geospatial data modelling enables to get an insight into the deepest regions of the ocean and to visualize areas with the most difficult access [9]. The morphology of the seafloor is caused by the uneven distribution of the elevations on the Earth with distinctly uneven hypsography. Thus, the majority of the depths on the Earth is occupied by the deep basins (ca. 4-6.5 km), while relatively few areas are covered by the shallow zones [10].
The tectonic plate boundaries are the hotspot areas where the most of the largest earthquakes take place and oceanic slabs descend beneath the continental lithosphere causing trench migration [11]. The largest earthquakes mostly occur in the shallow part of the subduction zones [12]. The earthquakes are sometimes accompanied by strong tsunami waves. Slab dynamics is one of the most important drivers for the trench formation, dynamics and migration [13]. Many factors affect speed and direction of the trench migration. These include plate-mantle coupling, slab interactions with the mantle transition, plate geometry, kinematics, strain rate, temperature, fluid pressure, deformation mechanisms, as well as mantle rheology [14]. The geological complexity of the Manila Trench is also expressed by the connection with the Philippine Trench where a horizontal mantle flow exists between the Manila Trench and the Philippine Trench. It is caused by the collision between the Palawan block and the Philippine Mobile Belt, and movement of the South China Sea slab [15]. Various factors cause the formation of the Manila Trench. Being an oceanic trench, it presents a special area of the ocean seafloor with distinct geomorphological structures characterized by the notable depths and steep gradient angles, located in the zones of the continental margin tectonic plates bending [16]. The convergence between the two plates forming the Manila Trench is roughly northwestward. There is a high-velocity zone present in the crust and upper mantle beneath the Luzon arc where the Manila Trench is stretching [17]. All these factors briefly mentioned to illustrate the tectonic settings of the study area indicate high frequency of earthquakes and special tectonic zone where the Manila Trench has been formed. The seafloor of the deep-sea trench presents a complex structure combined of various landforms: mid-ocean ridges, transform faults, ocean plains complicated by chains of seamounts, minor ridges, trenches and plateaus [18]. Factors influencing geomorphic development, structure and bathymetric patterns of the trench are diverse. To mention some of them: geological, hydro-chemical, geothermal, climatic, tectonic and bathymetric determinants [19].  The structure of the northern Manila Trench has been studied in various papers focused on the problems of crustal structure and deformation in the north of the Manila Trench [22][23][24]. Among other findings, the increasing dip angle of the Manila Trench from north to south has been discovered. The Manila trench loses its surface roughness in the early collision zone and gradually becomes a less well-defined deformation front [25]. However, the comparative analysis of its southern and northern parts is still missing. Therefore, current paper contributes to better understanding of the Manila Trench by comparative analysis of its geomorphology in the northern and southern segments. The data used in the current study included SRTM, ETOPO1, geoid and gravity grids (CryoSat-2, Jason-1). Topographic data used in this research include ETOPO1 Global Relief Model grid containing elevation data (topography/bathymetry) for the Earth.

METHODOLOGY
Data analysis applied for marine geology may be supported by programming languages, e.g. R [26], Octave, AWK [27], statistical libraries, e.g. Gretl [28] and other advanced tools, such as statistical software, e.g. SPSS [29]. For this study, a Generic Mapping Tool (GMT) scripting toolset [30] was selected as a main tool. It enabled to perform all the steps of the research: data analysis, modelling, mapping and statistical analysis and visualization. The functionality of the GMT and its powerful cartographic possibilities explain the choice for the GMT. The use of the GMT specifically for the Manila Trench was presented as mapping in various works, e.g. [31]. Examples of other advanced tools for geological data analysis include PHASER diffractometer, direct observations received from the Research Vessel (R/V) cruises, observations from Digital Elevation Models (DEM) [32]. The mapping was based on the available data sets embedded in the GMT [33] as well as Shuttle Radar Topographic Mission (SRTM) data [34][35]. The geologic characteristics was performed by modelling seafloor geomorphology, mapping tectonic slabs and submarine volcanoes.
The geoid gravitational regional modeling was done through the GMT modules 'psbasemap', 'grdcontour' and 'grdimage'. Plotting contours of the terrestrial and water areas, geoid and bathymetry, net grids and basic cartographic elements such as titles, scales, annotations, was performed for the mapping (Fig. 3). Each GMT module consists of the set of code lines, which, taken together as a script, produce and visualize maps. The gravity mapping shows free-air gravity anomalies covering study area including terrestrial and marine regions. Modelling was based on the data sets based on the CryoSat-2 and Jason-1 satellite-derived data [36], for the Manila Trench and surrounding regions (South China Sea, the Philippines). The gravity modelling revealed ground information for the geologic and geophysical settings of the study area, as free-air gravity anomalies reflect subsurface density variations [37]. The gravity dataset covering research area is based on the map gridded with a resolution of one arc-minute (Fig. 4). The gravity map was plotted using following GMT code (Code 2): gmt grdimage gravRT.grd -I+a45+nt1 -R105/123/8/24 -JY114/12/6.5i -CgravRT.cpt -P -K > $ps The GMT modules used for plotting 3D model (Fig. 5)  Using this code, the 2D initial meshing was divided into small simple polygons with curved contour lines and color palette highlighting the topography.
The geomorphological modelling consisted in plotting two sets of the cross-section profiles and then comparing their parameters (Fig. 6). The raster image as a background (Fig. 6C)  The graphs (Fig. 6A and Fig. 6B) were plotted using following GMT code (Code 7):

RESULTS AND DISCUSSION
Based on the patterns of the free-air gravity map (Fig. 4), the high density materials are distributed continuously from the northern the southern Taiwan vertically southwards to the Luzon Island, the Philippines from -50 mGal to -200 mGal. The majority of the values in the South China Sea areas lie between the 10 to 30 mGal, while values reaching >60 mGal can be noted in the island areas and southern Philippines. The significance on the free-air gravity differences for the two profiles consists in the principle of the free-air gravity map, which contains signals from the seafloor topography, sediments, and crust and mantle density anomalies. More specifically, the free-air anomaly is dominated by short wavelength variations which reflect the density contrast at the seafloor. Therefore, the differences highlight the unevenness in the geophysical properties in the study areas. Thus, the area crossed by the southern profile (Fig. 6, green transect) has dominant values of -50 to -220 mGal along the profile while northern profile (Fig. 6, yellow transect) has more values of -20 to 20 mGal in the surrounding coastal area. Asymmetric profile shapes and more gentle slope in the northern profile ( Fig. 6) can point at local geophysical variations, e.g. pulsing, dehydrating, and radially flowing mantle plume. Furthermore, bathymetry and gravity data demonstrate that V-shaped profiles in the deepest parts of the cross-section (<-4,000 m) may be the result of sublithospheric flow. It is formed because of the displacement of large volumes of magma at the trench axis which is reflected in its geomorphology as east-west asymmetry (Fig. 6).
The statistical results (Fig. 7)  However, a gradually slanting terrace on the southern part of the trench can be seen on the fragment 50-100 of the cross-section with depths from -1,100 to -1,300 m. Hence, the comparison between two segments of the trench located on the northern and south-eastern part of the trench show that northern part of the trench has shallower values of depths comparing to the southern yet more steep slope on the Luzon Island side.
The approaches to the data analysis, visualization and geological modelling are diverse. To mention some of them: tomographic, seismic and bathymetric 2D and 3D modelling, geologic cross-sectioning, seismic cross-section profiling, R programming language for analysis of correlation between variables [38,39], ArcGIS based assessment and calculation of measured data [40][41][42], ILWIS GIS [43], Python statistical libraries, such as Matplotlib, NumPy, SciPy, Seaborn, Pandas, StatsModels [44,45]. Using special software for automatization of the geodata can be illustrated by Autotrace [46], R packages [47], [48] or geological software [49,50]. In view of a variety of the approaches and methods illustrated above, the advantage of the GMT consists in its embedded statistical module that enable to visualize data and to perform descriptive statistical visualization and modelling as histograms showing data distribution.

CONCLUSION
The presented research sets out to demonstrate the utility of the Generic Mapping Tool (GMT) for cartographic visualization in support of interpreting geophysical data of the seafloor geomorphology of the Manila Trench. This investigation of the tectonic dynamics of the Manila Trench is an important problem in the field of geodynamics and relevant for understanding the potential origins of hazards related to earthquakes and tsunamis in the region. The study produces useful visual cues to interpret the tectonics of the crustal structure based on the derived data models such as the geoid, free-air gravity anomalies, bathymetry, and geology aimed at modelling and visualizing tectonic lithospheric settings and geological situation in the study area.
The GMT cross-section stacking methodology [51][52][53] proved to be a successful means of visualizing and plotting geomorphological models applied for the submarine bathymetry by effectively minimizing the hand-made cartographic routine. Using GMT enables to automatically digitize profiles, perform statistical analysis by plotting histograms and model trends for the general shape of the profiles using mathematical models of the lines approximation. Several GMT modules were tested and applied for the cartographic visualization of the marine free-air gravity, geoid, bathymetric mapping, geomorphic modelling and statistical data analysis. Two different trench segments were analyzed aimed at detecting their spatial changes for an area covering Manila Trench and adjusting Luzon Island (northern Philippines). The results provide useful recipes on GMT scripts that produce visualization outputs for interpretation and explanation of the tectonic dynamics of the Manila Trench. The application of GMT scripts, produced in this paper, can easily be adapted to other investigations using similar datasets. This paper also aims to close the knowledge gap on the geomorphological differences between the northern and southern part of the Manila Trench. Therefore, the goal of the GMT-based mapping was to relate geomorphological mapping to underlying tectonics conditions. These fine-resolution data such as SRTM grid was used to map bathymetry and topography of the study area and 1-min resolution from the gravity data (Jason-1, CryoSat).
Future work should consider applying other GMT modules and their combination. Developing a complex approaches of the cartographic-geophysical investigations include earthquake intensity modelling, mapping detailed volcanic arcs. The combination and overlapping of various geospatial data as layers aimed to examining and quantifying geomorphological submarine landforms in connection with tectonic settings would enable to get more information of these data using advanced geospatial data analysis. Other recommendation for the future research could be integration of the GMT with programming and advanced statistical packages that absent in GMT: clustering, dendrogram plotting, factor analysis, Principal Component Analysis, etc.