**Geologic development** Geologic development of the MT took place in several phases shaped by the influence of various factors: geologic conditions, volcanic submarine areas, hydrological settings of the currents and sedimentation. Seafloor bathymetry of the trench varies significantly in seaward axis. The Caroline Ridge is located off the trench axis at a distance of ca. 250 km. A prominent trench-parallel belt of seamounts with trench-perpendicular width of ca. 250 km, is located <1350km. Another prominent group of wide seamounts intersects the trench at distance of 1600-2000 km [774]. The effects of the volcanism on the biological communities [321] proved correlation between earthquakes and the assemblage structures. Because rapid sedimentation is induced by turbidity currents, the organic matter content increased in the topmost layers which indicate that disturbances to the deep-sea sediments affects vertical distribution structure of meiofauna in the trench.

The seabed structure of the MT is composed in the following way from top to bottom: a sedimentary cover, a basalt of the leitic composition, a complex of parallel dykes of diabases, an isotropic gabbro, a banded gabbro-ultrabasic complex lying on the mantle ultrabasites [77]. The trench advances toward the upper plate corresponding to the subduction of very old, Mesozoic, oceanic lithosphere. Mariana Arc lavas are enriched in Molybdenum (Mo) and have 89/95 Mo >Mid-Ocean Ridge Basalts (MORB). This implicates the addition of a Mo-rich fluid with 89/95 Mo ~ +0.05% to +0.05% to the mantle wedge beneath most islands [175]. Thin Quaternary clayey sediments overlie a 400 m thick alternating sequence of the Early Cretaceous quartz-trachyte pillow lava and Early to Late Cretaceous radiolarian cherts in the oceanward slope of the MT [244]. MT belongs to the one of the deepest 45% of the oceans depth range, and its unique topography disrupts the continental shelf-slope-rise to abyssal plain continuum, resulting in an array of deep isolated habitats. The often abrupt and distinct topography is likely to further promote speciation through high hydrostatic pressure, remoteness from surface derived food sources and geographic isolation [159].

**Lithology** Lithological structure and sediment thickness are among other factors affecting MT formation. Relatively few other trenches, comparing to MT, also formed on the boundaries of the tectonic plates, but usually at a greater distance from the continents in the rift zones. Therefore, sedimentation processes within MT are deeply impacted by the location near the rift zones that are mainly associated with the formation of the underwater mountain ranges and spreading. Their expansion in the sides of the neighboring lithospheric plates is a result of the rise towards MT. The outflow to the surface of the ocean floor of the substance of the upper layer of the Earth's mantle. As a result of the unevenness of this process, the transversed faults are formed that have excellent conditions for sedimentation accumulation. According to available geophysical data [145], the oceanic crust of the tectonic plates surrounding MT is composed by a number of layers. Therefore, the impact of the sediment subduction on the trench dynamics is relatively high [253].

MT is the largest structural trap located in the continental margins of the Pacific Ocean: the sediments are being carried by the ocean waves in a clockwise direction, passing through the trenches on the west of the Pacific (i.e. Kermadec Trench and the TT through the Samoan Passage) flowing northwest across the equator to the MT. Its geomorphic conditions create ideal conditions for the sedimentation intercepting almost all clastic material originating from the adjacent land. Accumulation of the terrigenous sedimentation on the MT seafloor depends on various factors: climatic and hydrodynamic local settings, tectonics, seafloor bathymetry, geomorphology, type and size of particles. Large areas of active sedimentation can be distinguished where speed and volumes are very high, as well as local areas of active sedimentation, or on the contrary, dominating erosion areas. Sedimentation occur in various sediment traps, with some areas noted by areal sedimentation. The processes of the sedimentation within the MT are highly uneven.
**Geomorphic structure** MT has a complicated geomorphic structure with various types of the continental tectonic plate boundaries producing a unique type of its seafloor bathymetry. The seafloor of the trench is a background, on which all the processes are reflected [131]. Main morphostructure elements of the seafloor are the outskirts of the continents and of the ocean bedrock. The structure of the trench and the nature of its relief are greatly complicated by the multiple secondary tectonic disturbances, i.e. by the occurrence of faults and displacements on of graben, horsts and lateral geologic shifts [507].

Among other trenches, MT is distinct for its edge type associated with the marginal tectonic plate subduction processes [69]. It is formed in the process of plates subduction. Its submarine terrain has a complex structure, the scale of which sharply exceeds the dimensions forms of the Earth's surface on the margin continents. These include the system of the mid-ocean ridges, crossing transform faults, deep-sea ridges, vast ocean plains complicated by chains and mountain groups and plateaus [293]. The width of the continental shelf of the margin plates connecting to the MT occupies varies from a few tens to thousands km with average bottom slope 17-25° [606].

Areas with seamounts at the trench axis are marked by gentle trench geomorphology, smaller axial vertical loading, and bulgy topography at the outer-rise. Mariana Arc and West Mariana Ridge border the trench on the west side. The Caroline Ridge and a group of Caroline Islands are located southwards off the trench axis at distance of about 0-250 km. A prominent trench-parallel belt of seamounts with trench-perpendicular width of about 250 km, is located at around 1350km BSL. Another notable group of wide seamounts, intersects the trench at a distance of 1600-2000 km BSL [774]. Morphological structure of the MT has following pattern: a majority of the seafloor is occupied by the abyssal depths (3-6 km) followed by the extreme depths (> 6 km). Depths of > 6 km are restricted mostly to the deepest part of the trench located in the south-western part of the trench (area of Challenger Deep), although depressions of 6-7.5 km BSL, occasionally occur in the central part of the trench crescent [557]. Typical depths of the MT seafloor are 4-5 km [422], [424].

Seafloor of the MT is composed by the oceanic crust forming rift zones of the mid-ocean ridges with a capacity of 5-10 km [557]. MT is complicated and constituent of the interrelated factors forming its tectonic structure. Deformations of the trench respond to the coupling between the upper and lower plates: coupling relating to the continental slab age-buoyancy [642]. Back-arc deformation roughly correlate with upper continental tectonic plate velocity [240]. Trench migration rates are chiefly controlled by lower continental tectonic plate velocity [350], which depends on the tectonic slab age buoyancy [92]. Main ridge axes of the MT represent a place where the upwelling magma encounters cold seawater in the location of tectonic plates. Hydrothermal cooling of magma at the ridge axes causes recycling of the entire volume of the water in the trench column. In this way it provides nutrients for biological communities in the abyssal zone of trench [353], [522], [269]. Examples of interconnection of the system are illustrated by biological communities [186], [348].

Geomorphic structure of the MT is complicated by the longitudinal ridges, steps, large landslide, bodies and ledges. Transverse profile of the MT is strongly asymmetric: slopes are higher on the side of the island arc, dissected by deep underwater canyons with narrow steps. MT has complicated steps of various shapes and sizes, caused by active tectonic and sediment processes. Continental slopes are high, reaching up to several thousand m and inclining to 3-6° (in the SW part up to 30-40°), upper boundary of which coincides with the edges of the shelf (depths 150-200 m). Slopes of the passive margins are strongly complicated by the terraces, ledges, marginal plateaus and canyons. Slopes of the margins are steeper, reaching can 5-7 km in height. Seafloor geomorphology with very steep escarpment of the narrow MT does not play a crucial role for ocean water mixing [223]. Interplay between topography and turbulence of the ocean waves includes impacts of the inertial and meridionally directed tidal inertia-gravity waves, responsible for turbulence and variation of water masses mix according to the depth stratification. At the base of the continental slope of passive margins formed by the adjacent tectonic plates, continental foot forms accumulative body. In turn, it is formed by merged cones of the removal
and plumes of the suspension flows and submarine landslides with abyssal sedimentation.

<figure><img src="image_56.png" /><figcaption>**Map 2.16.3:** Earthquakes event map at PSB area: prime hypocentres and magnitude values</figcaption></figure>

Submarine margins of continents, or transition zones, depend on the features of the relief and geology. Accordingly, they are divided into 2 types: the passive and active ones. The 1st ones include the shelf, that is mainland bank, the slope and the mainland foot. Having complicated relief structure with marginal seas, island arcs and deep-sea gutters, MT belongs to the active type of trench. Due to the geodynamic reasons, MT motion, migration and upper plate deformation can be described as the result of the response to the mantle flow and imbalance between the forces exerted by lower and upper continental plates at the plate interface.

**Cracks** MT is notable for well-developed systems of cracks in the structure. Cracks of the MT are abundantly developed at the edge, in the middle and upper parts of the terrace [240]. The formation of these cracks is caused by gravitational instability: horizontal extension is due to the trench-ward increase of the dip angle of the seafloor surface. Besides, it cause slope instability due to gravity pull. The cracks on the surface of the seafloor on the oceanward slopes of the MT are attributed to tensile rock failure. Cracks are elongated nearly parallel to the strike of the cliffs, some are aligned, branching or merging have no vertical nor lateral displacement. This can be explained by open tension fractures. Surface edges of the MT cracks are sharp, indicating young origin with a noticeable pressure ridge caused by mud
overflow during closure of cracks. Cracks are formed at topographically specific areas, at the horizontally
stretched surface of the down-going plates under tensional stress, on the horizontal or very gentle slopes
above steep cliffs. Cracks are caused by a combination of gravitational slope instability at the edges of
the slope or top of the ridge [280] and additional inertia during earthquake shaking [129], [658].

**Submarine earthquakes** Deep earthquakes are important factors causing system of cracks of the MT. Transformation faulting is the mechanism for deep earthquakes [319]. Experimental observations [209], [75], [318] indicate a shear instability that develops during the incipient transformation of metastable olivine to wadsleyite or ringwoodite. Horizontal tensional stress during the earthquake shaking triggers development of the system of cracks in the MT, as stress difference under static conditions is caused by the slope inclination. The majority of the cracks within the MT are open tensile failures and are situated in areas of very gentle slope or horizontal parts of the seabed, in most cases near the edges of steep cliffs which delimit the terraces. There is no lateral nor vertical dislocations on both sides of the cracks, so they are open fractures.

**Igneous volcanic areas** MT is important integral feature of the active continental margins that naturally includes deep plate subduction. Subduction of slabs of oceanic lithosphere into the deep mantle involves a wide range of the geophysical and geochemical processes and is of major importance for the physical and chemical evolution of the Earth [546]. Thus, it is noticed [277] that the continuity of the slab across gaps in the Wadati-Benioff seismic zone between deep earthquake clusters and shallow seismicity and they argue that only the deep earthquakes beneath the north Fiji basin and the deep events beneath New Zealand may occur as detached events with no mechanical connection to the surface [388]. Using the observation and triggering mechanism of high-frequency oceanic T-waves, [28], it has been argued [213], [243] for the mechanical continuity of most slabs across the deep seismicity gap.

Recent studies revealed [165] that MT advance toward the upper plate corresponding to the subduction of very old, Mesozoic oceanic lithosphere. The Mariana Arc lavas are relatively enriched in Molybdenum (Mo) and have 898/95 Mo significantly greater than MORB (apart from the samples from the island of Agrigan). This implicates the addition of a Mo-rich fluid with 898/95 Mo ca. -0.05% to +0.05% to the mantle wedge beneath most islands [175]. Recently updated analyses of global plate motions indicate that significant trench advance is also rare on Earth, being largely restricted to the Marianas-Izu-Bonin arc [101]. The effects of trench migration on the descent of subducted slabs are discussed by [211]. The global plate tectonics movements shaping the seafloor bathymetry can be described as follows.

The seafloor spreading fabric in the MT is obscured by the absence of magnetic anomaly lineations. Elongated topographic structures accompanied by plate bending exist in the oceanward slope of the Challenger Deep [774]. The Japanese Lineation set is a Mesozoic magnetic anomaly lineation, it exists in the East Mariana Basin, east of the MT [566], which means that the subducting plate along the southern part of the MT is a part of the Pacific plate formed in the Mesozoic.

Other remarkable topographic feature of the MT includes a set of many elongated ridges and escarp-
ments accompanied by the plate bending on the outer slopes with two existing distinct strikes for these
structures. The elongated structures have two strikes: N85°E and the same as the trench axis; N70° E.
Most bending-related structures except those near the trench axis, were formed by the reactivation of the
seafloor spreading fabric. The presence of cracks is another distinguished geological feature of the MT
[580]. MT has many open prominent cracks on its surface, up to a few meters deep and a few hundred
meters long on the diatomaceous clayey sediment surfaces of its oceanward slopes. The shape of the
cross-section is Y-shaped. More precisely, the bottom is a very narrow (1-3 km), yet in some places wide
(several tens of kilometers) flattened.
**2.17 Izu-Bonin Trench**

The Izu-Bonin Trench (IBT), also known as Izu-Ogasawara Trench is a hadal trench located in west Pacific Ocean and lies south of Japan (Map 2.17.3), [405].

<figure><img src="image_57.png" /><figcaption>Map 2.17.1: Izu-Bonin Trench: bathymetric map</figcaption></figure>

Longitudinally, it extends from near Tokyo, Japan to Guam [30]. In its structure, the trench continues the Japan Trench. The **IBT** is formed in the area of active tectonic subduction of the Cretaceous (135 Ma) Pacific Plate beneath the Philippine Sea Plate [652], [405]. As a consequence of this subduction the following structured are formed: a volcanic arc of the Izu and Bonin Islands on the Izu-Bonin-Mariana Arc system and the Shikoku Basin as a backarc basin (Map 2.17.2). The **IBT** is boarded on the east by the Ogasawara Plateau [573], [405].

In its historical development, the Izu Bonin arc has two major periods of arc formation that differ in duration, composition, and volcanology. The first lasted ~5 to 6 million years and produced voluminous boninitic-tholeiitic magmas. The second post-42 Ma Izu Bonin arc formation period produced a classical volcanic arc with stratocones [30]. The formation of the Shikoku backarc basin took place in the early Miocene, ca. 30-15 Ma [707]. A new period of along-strike arc rifting began about 2.5 Ma ago [279]. The historical development of the oceanic trench is mirrored in its seabed structure, which is well structured, including valious layers: a sedimentary cover, a basalt of the leitic composition, a complex of parallel dykes of diabases, gabbro and gabbro-ultrabasic complex lying on the mantle ultrabasites.
<figure><img src="image_58.png" /><figcaption>**Map 2.17.2:** Izu-Bonin Trench: geologic map</figcaption></figure>

The trench extends to nearly 10 km depth. The 200 km wide forearc region includes inner trench slope with a basal terrace and seamounts, outer arc connected with the Bonin Islands and a thick forearc sedimentary basin east of the active arc volcanoes. The backarc region encompasses active rift basins and cross-chains of submarine arc volcanoes extending into the Shikoku Basin [127]. Northern end of the Izu-Bonin arc is a juvenile oceanic arc colliding with the Honshu arc, a mature island arc since Middle Miocene Saito and Tani, 2017).

A chain of serpentine seamounts is exposed at the forearc slope 40-50 km west of the IBT axis. Seismically, few large earthquakes (> M7.0) [306] are detected at depths < 100 km (red circles on Map 2.17.2) and many large earthquakes at depths 100-300 km (green circles on Map 2.17.2). The Boso triple junction located on the north of the Izu-Bonin Collision Zone (Map 2.17.2) is the northern end of the arc colliding with the Honshu Island arc as a consequence of the northwestward migration of the Philippine Sea Plate [579], [405].

The IBT is boarded by the Philippine Sea tectonic plate in the west. The Philippine Sea is a geologically complicated region notable for complex bathymetry and plate shape [753] and a complex patchwork of seafloor ages. The bounding subduction zones of the Philippine Sea Plate exhibit significant NS variation in slab morphology [101]. As a result of the interplay of these factors, the trench is consisting of the Izu (north) and the Bonin (south) separated roughly by the Ogasawara Plateau. Thus, seafloor ageing affects advancing trench through the subduction of progressively older and stiffer lithospheric
material [165], which was further discussed by [281] who pointed on asymmetric spreading caused by
ridge migration in the Izu-Bonin region in view of age difference between the oldest forearc basalt and
the crust of the Amami Sankaku Basin.

<figure><img src="image_59.png" /><figcaption>**Map 2.17.3:** Izu-Bonin Trench</figcaption></figure>

More detailed description of the Izu-Bonin tectonics, geophysical, petrological and geochemical set-
tings, volcanism and seismic characteristics, analysis and interpretation of geomorphological, geologi-
cal, and thermochronological data, are reported in relevant publications [270]; [18]; [670]; [587]; [228];
[202]. Hence, the region of the Izu-Bonin-Mariana Arc is now one of the best and intensively studied
oceanic arcs. However, lesser attention is given to the submarine geomorphology of the **IBT** as a deep-sea
geomorphic structure. Referring the lack of such research, current study aims at visualizing the trench
by mapping and modelling its geomorphic form using cartographic approach of the **GMT**.

Morphologic development of the **IBT** took place in several phases shaped by the impact of various fac-
tors: geologic conditions, volcanic submarine areas, and inner ocean currents and sedimentation. Factors
affecting formation, geomorphic development and bathymetric patterns of the trench are diverse. The
most important ones include tectonic development, geological, hydro-chemical, biological, geothermal,
climatic, tectonic, bathymetric and geomorphological determinants, deep ocean currents that bring sedi-
ments to the trench bottom and contributing towards accumulation of the sediment thickness. Impact of
lithosphere is illustrated by a constant exchange of matter and energy between the submarine volcanoes
located nearby [266]. The hydrosphere influence on the **MT** is reflected by deep ocean currents bringing
sediments to the trench bottom and contributing towards accumulation of the sediment thickness [253].
The trench is located in the collision of the tectonic plates. It has a distinctive morphological feature of
the convergent Pacific Plate boundaries, along which lithospheric plates move towards each other. Not
all seismic zones in subducted slabs of the **IBT** are continuous and that some deep events seem to occur
isolated from those at shallower depths.
## 2.18 Japan Trench

The Japan Trench (JT) is located in western part of the Pacific Ocean, eastwards off Japan (Map 2.18.1) [396]. The geologic and tectonic settings in the region of the JT are strongly influenced by the subduction of the Pacific Plate beneath northeastern Japan and North American Plate along the JT (Map 2.18.2), which is well described in the previous works [261], [393]. The region of the JT is a convergent zone associated with many earthquakes [762]. It belongs to the circum-Pacific seismic belt (a.k.a. 'Ring of Fire') where large earthquakes originate within the areas of sedimentation on upper forearc slopes [17], [396], [393].

<figure><img src="image_60.png" /><figcaption>Map 2.18.1: Japan Trench: topographic map</figcaption></figure>

The structure of the trench presents a typical arc-trench system [520]. In terms of tectonic geomorphology, the arc-trench system is a type of the mobile belts with bathymetric and topographic relief formed under the impact of tectonic plates subduction, active volcanism, repetitive earthquakes and orogenesis. The geological complexity of the JT results in the high seismicity of the area, which is well studied in various works [565], [259]; [500]; [67], [393]. The geology of Japan trench-arc system consists of Jurassic to Neogene accretionary prisms and Neogene to Quaternary sedimentary and volcanic rocks caused by volcanism during the opening of the Sea of Japan in the middle Miocene [730], [396].
<figure><img src="image_61.png" /><figcaption>Map 2.18.2: Japan Trench: geological map</figcaption></figure>

The Cretaceous-Paleogene granitic rocks are exposed in several regions, especially on the fore-arc side. High seismicity of the region can be illustrated by several large interplate earthquakes that have affected the Tohoku region along the JT during the last century [752]. Because JT is located on the tectonic plates subduction zone, the earthquakes here belong to the largest and most destructive earthquakes in the world [292]. Many open cracks on the diatomaceous clayey sediment surfaces of oceanward slopes of the northern JT were detected by [580]. The surface of the seabed outside the cracks (i.e. normal seabed) appears to be coated by ferroman-ganese oxides [283], and is overlain by fine pelagic mud, containing dominantly diatoms, silicoflagelates, radiolarians and sponge spicules, together with a considerable amount of detrital quartz and ash fragments.

Northern and southern parts of the JT demonstrate variations in their geologic features with heterogeneously distributed seismicity around the JT. Most of the large earthquakes with magnitudes > 7 occurred in the northern forearc region of the trench. For example, the Tohoku-Oki earthquake of 2011, is the largest earthquake that took place in forearc region of the Northeast Japan subduction zone, west off the JT [356]; [274], through the shallow part of the subduction zone to the seafloor surface. It triggered a devastating tsunami with massive damages along the northeast coast of Honshu, Japan [292].
<figure><img src="image_62.png" /><figcaption>**Map 2.18.3:** Japan Trench: geoid model map</figcaption></figure>

Southern part of the trench is seismically more 'quiet': no large earthquakes were recorded during recent centuries [552], except for the 1938 Shioya-oki earthquake. The subduction zone characterizes northern part of the trench [718]. In a 60 km buffer gap off the trench axis, a decollement structure with soft sediments is located. The zone of this structure between two tectonic plates decouples shear stresses, which in turn, forms an aseismic zone. Comparing two devastating earthquakes: the 1968 Tokachi-Oki and the 1994 Sanriku-Haruka-Oki, two distinct heterogeneities along the plate boundary at 40°10′N, located exactly at the southern boundary of the two earthquakes were detected [230]. This well illustrates the direct affect of the plate tectonics triggering earthquake events that naturally lead to significant deformations in the geomorphic landforms. It further explains the coincidence of the focal regions and the structure with slower velocities in the lower crust.

The subduction of the seamount chain penetrating to the Pacific Plate possibly take place. It can further create irregular geomorphic landforms. When these seamounts subduct deep, they collapse due to the mechanical affects from the free water released by a collapsing seamount. In turn, it may rise into the island-arc mantle and the lower crust. This brief illustration of the deep correlation of the tectonic and geological processes that result in actual land forms and geomorphological patterns. Seismic
investigations of the JT axis shown [66] variations in the frictional resistance along a décollement on the structural development of the frontal wedges near the subduction zones.

<figure><img src="image_63.png" /><figcaption>Map 2.18.4: Japan Trench: gravity model map</figcaption></figure>

The observed slope sediments between the inner and outer wedge may have accumulated during the growth of supercritical wedge. Seismic surveys indicate [334] that shallow areas at the Tohoku earthquake experienced significant deformation. During this event multiple faults slipped above and below the pelagic clay layer [610].

The geophysical features of the JT area were discussed [500]: comparison of the JT Nankai subduction zone for the distribution of the stress state [396], [393]. Based on the correlations between strength and temperature, the deep seismicity is likely to result from plastic instabilities, with the distribution of earthquakes being related to strength distribution and therefore slab mineralogy. A complex rheological structure for subducting slabs has been proposed by Karato, Riedel, and Yuen [308], due to the effects of grain size reduction during phase transformations. Based on their model, rapidly subducting and therefore relatively cold slabs should be weaker than slowly subducting warm slabs. Not all seismic zones in subducted slabs are continuous and that some deep events seem to occur isolated from those at shallower depths [277].
<figure><img src="image_64.png" /><figcaption>**Map 2.18.5:** Japan Trench: 3D mesh model map of geomorphology</figcaption></figure>

Structure, magnetism, and dynamics of subduction zones, the continuation of the low wave speed, high attenuation (low Q), and seismically anisotropic wedge are normally near 400 km depth beneath some back arcs. This may be explained by the persistence of hydrous phases to that depth and that magnetic systems are not limited to near-surface regions of the mantle. The correlation also exist between the dehydration of the slab, rupture nucleation and crustal earthquake triggering.
**2.19 Puysegur Trench**

Puysegur Trench (PsT) has ca. 800 km in length, stretching from the most southern tip of the South Island to the Macquarie Triple Junction [387].

<figure><img src="image_65.png" /><figcaption>Map 2.19.1: Hikurangi, Puysegur and Hjort trenches, New Zealand: GEBCO based topographic map</figcaption></figure>

The Macquarie Fault Zone is associated with the Alpine Fault, which is the right-lateral transform fault boundary separating the PsT and the northern KT [37], [387]. The Macquarie Ridge Complex (MRC) is a dominantly translational structure accommodating right-lateral motion between the Australian and Pacific Plates and PsT, where Australian lithosphere subducts beneath the southwest corner of New Zealand on the Pacific Plate. The plate boundary south of New Zealand is located along the MRC which is crucial for understanding the tectonic structure of the region. Comparing to the Hikurangi, PsT and Hjort are younger (started at 168 Ma) and less developed [700], [387]. The Indo-Australian Plate subducts beneath the Puysegur Bank and the Fiordland Massif forming PsT.

Hjort Trench continues PsT southwards from the Macquarie Island (Fig. 2.19.1) [387]. It presents a linear topographic depression south of Macquarie Island in the southwest Pacific Ocean (Fig. 2.20.1).
Hjort Trench lies in an area of transpression where plate boundary transitions from a transform boundary to a convergent one. Hence, this region shows how a transform boundary with a vertical (near-vertical) transform fault becomes an area of under-thrusting [540]. The deepest point of Hjort Trench is ca. 6.3 km b.s.l. Eastwards, the Hjort Ridge follows the general curve of the Hjort Trench, separating it from the Hjort Plateau, [387].

<figure><img src="image_66.png" /><figcaption>**Map 2.19.2:** Hikurangi, Puysegur and Hjort trenches, New Zealand: geoid model</figcaption></figure>

Unique tectonic setting of the region consist in the connection of two subduction zones in the north (Hikurangi margin) and in the south (Puysegur margin), connected by an oblique continental collision along the Alpine Fault located in the South Island of New Zealand. This cause variations in the sub-marine geomorphic structure of the trenches which is also referred as 'subduction scissor' [608]. The
Puysegur/Hjort subduction is highly oblique (dextral) and directed southwards. The Hikurangi subduction is directed northwestwards. Hence, the South Island is caught in between by the 'subduction scissor'. In the absence of available high-resolution information and detailed mapping on the seabed geomorphology of the Hikurangi, Puysegur and Hjort trenches, a GMT approach to map their seafloor on a local scale were undertaken. Since seafloor of the deep-sea trench represents an integrated response to a complex range of factors, such as oceanographic, geological and tectonic conditions which eventually define the present shape of the seabed landforms, an integrated approach was used to visualize, model and map three trenches on the New Zealand margin [387].

Hjort Trench is characterized by an oblique subduction zone geoid anomaly. Large 1924 earthquake occurred in the northern Hjort Trench which suggest a thrust focal mechanism and active subduction of young oceanic lithosphere beneath older oceanic lithosphere [87]. Frequent seismic events, most less than 20 km deep, characterize the transpression along this plate boundary [515]. The southern segment of the Alpine Fault on the South Island continues into the Tasman Sea Basin and finishes at the PsT (Fig. 2.20.1). The oceanic Indo-Australian Plate subducts eastwards through the PsT under the Pacific Plate [387].

Geologic events during the tectonic historical records of this region, such as Eocene-Oligocene rifting as well as Miocene-Quaternary strike-slip, crustal thinning in the Solander Basin, development of the oblique subduction zone in the PsT, are studied by processing seafloor relief data and interpretation of the seismic and bathymetric cross-section profiles acquired over the PsT area [541]; [351]; [352]. Effects of slab mineralogy and phase chemistry on the subduction dynamics (buoyancy, stress field), kinematics (rate of subduction and plate motion), elasticity (deformation and seismic wave speed), thermometry (effects of latent heat, isobaric superheating) and seismicity (due to adiabatic shear instabilities) are well discussed [56]. The morphology and development of the subducted slabs are more complex than modeled [97], on either whole mantle flow or convective layering at 660 km depth. Apparently, slabs can deflect horizontally in the upper mantle transition zone beneath some convergent margins whereas penetration to lower mantle depths can occur beneath other island arc segments. Different styles of the subduction across the upper mantle transition zone and trench formation.

Various physical and chemical processes interact during tectonic plates subduction, and the complexity of the system necessitates an interdisciplinary approach involving seismology, mineral physics, geochemistry, petrology, structural geology, rock mechanics, and geodynamic modeling. For example, phase transformation kinetics determine buoyancy, rates of subduction and, therefore, thermal structure, with the latter feeding back to affect the kinetics.

Rheology and large scale slab dynamics is affected by mineral transformations and their kinetics. The rheology of mantle minerals, though having considerable importance in controlling the subduction process, is exceedingly difficult to study experimentally, due to the experimental limitations, since rheology at pressures of the transition zone and lower mantle [309]. An increase of intrinsic density or viscosity with depth as well as phase transitions with a negative Clapeyron slope can all inhibit or delay deep subduction, tectonic conditions and relative motion of the tectonic plates at Earths surface play important role. The effects of subduction of crustal rocks during continental collision consists in following subduction to depths >90km. Slices of continental crust exhumed back to the surface very rapidly, thus ensuring the survival of the high-pressure minerals.
**2.20 Hikurangi Trench**

The Hikurangi Trench (HT) (sometimes also referred as trough) is located on the Hikurangi margin, 500 km long 480 km wide [572], northeast from the North Island [387].

<figure><img src="image_67.png" /><figcaption>Map 2.20.1: Hikurangi Trench, New Zealand: tectonic and geological settings</figcaption></figure>

The Pacific Plate subducts beneath the Indo-Australian Plate at the **HT**; Indo-Australian Plate subducts beneath the Pacific Plate at the Puysegur and Hjort trenches, Fig. 2.20.1. Two subduction margins are
joined by the dextral transform Alpine Fault west off New Zealand. The area selected for the dataset lies
in southern Pacific Ocean and encompasses region off New Zealand coast (Fig. 2.19.1) [387].

<figure><img src="image_68.png" /><figcaption>**Map 2.20.2:** Hikurangi, Puysegur and Hjort trenches, New Zealand: Marine free-air gravity anomaly</figcaption></figure>

This area is characterized by a geologically complex history and current situation with oblique sub-
duction zone and subduction scissors. Such complex bathymetry exhibits considerable geomorphological
variability among the Hikurangi, Puysegur and Hjort trenches: there are differences caused by seamounts
and seafloor obstructing natural curvatures of the trenches. The toponymy of the HT vary in the existing
literature due to its specific geomorphology: in some works it is referred as 'trench', e.g. [102]; [619];
[620]; [621], in other works - as 'trough, e.g. [31]; [210]; [295]; [296].

The tectonic location of the Hikurangi is at the southern end of the 1000 km long Tonga-Kermadec-
Hikurangi subduction systems, where Pacific Plate subducts obliquely southwards beneath the Indo-Australian Plate [32] and Indo-Australian Plate has a strike-slip motion at 4 cm/yr, Fig. 2.19.1 [387]. Comparing to the KT reaching up to 9000 m, HT is shallow (ca. 3000 m).

Pacific Plate dips at **HT** at a gentle angle of about 3° for at least 100 km beneath the trench and then further steepens beneath the North Island [29]. The Alpine Fault stretches sub-parallel to the western edge of the **HT** (Fig. 2.19.1). The speed of subduction on the Hikurangi Margin varies with a general trend of increase northwards: 34 mm/yr in the south (westward off Bounty Trough), 38 mm/yr in its central part (North Island) and 47 mm/yr in the north (North Island) according to [620]. The Hikurangi north-west dipping subduction zone is old (started ca. 20 Ma). The subducted slab has a well defined Wadati-Benioff zone and steepens from shallow angles near the surface to near-vertical at depths >100 km [26].

The geophysical and geological settings of the Hikurangi margin vary along the trench from south to north [731]. Thus, interseismic coupling is strong and deep in the south, then gradually shallows northwards and becomes shallow and weak in the north. Southern region is characterized by accretionary wedge and tectonic contraction in the forearc. Northern region is characterized by frontal subduction erosion, numerous seamounts, and a slightly extensional upper plate.

Kinematics of the Hikurangi margin studied by Global Positioning System (GPS) measurements of the horizontal velocity fields [572] shown changes in faulting along the margin and clockwise tectonic block rotations caused by collision of the Hikurangi Plateau and Chatham Rise with **HT** and subduction of oceanic crust northeast of the North Island [387].

The Hikurangi Plateau initially formed a part and now separated from the largest oceanic plateau on Earth, the Ontong Java Plateau, located in the southwestern Pacific Ocean [569]. The Hikurangi Margin and continental shelf of the North Island of the New Zealand is well studied through various geologic investigations, more specifically: bathymetric, hydro-acoustic investigations [31], active tectonics. Palaeogeographic reconstruction of the region of Hikurangi introduced by [296] with respect to the evolution of its current topography/bathymetry was done using advanced methods of the thermochronology.

Slab dynamics is one of the most important drivers for the trench formation, dynamics and migration. 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 [166]. Understanding mechanisms of trench migration (retreat or advance) is crucial to characterizing the driving forces of Earths tectonics plates, the origins of subducting slab morphologies in the deep mantle, and identifying the characteristics of subduction zones systems, which are among the fundamental issues of solid Earth science [761].

Further findings in exploration geology in the **HT** with a special focus on methane seepage and gas hydrates in the underlying sediments are relatively well studied (e.g. [210]; [732]; [323]; [591]; [337]; [338]; [534]); studies on past (Holocene) subduction earthquakes, coseismic coastal deformation and tsunami [102]; [50] which later on largely affected submarine geomorphology of the Hikurangi margin forming submerged terraces [49].
# Chapter 3

## Methodology

Current dissertation contributed towards the development of the advanced approaches of data analysis demonstrating the combined use of GMT, Python, AWK, R and Octave/Matlab. IT approaches in Earth sciences are effective method. Statistical data analysis was implemented by ML algorithms. Seamless integration of the multi-dimensional data by R and Python statistical packages enabled to perform marine geodata analysis and proper selection of the optimal algorithms to process the data features. The fully provided codes written on the programming language make this research easily repeatable and accessible for other scientists. The adoption of the open source programming rapidly replacing the commercial technologies, Python and R will soon become popular among the marine geologists. This research contributes towards the development of the technical implementation of the ML algorithms for effective processing of the geodata.

### 3.1 Data

DATA used in the dissertation for modelling and visualization include raster and vector types. Data are listed in details in Table 3.1.1. High-quality geospatial data sets are crucial for the geospatial data analysis in a big research project. The importance of the quality of the input data is explained by the fact that it affects the final results: high-precision data always results in the reliability of the output. Therefore, the selection, assessment and integration of relevant data was initial step of the cartographic project. Selecting correct data significantly improves the quality of the results, and the interpretation of the trench geomorphology. The more it is actual for the bathymetric and geomorphic mapping where precision is the key factor for the research output [686], [744], [401].

Raster data include bathymetric and topographic grids. Based on an analysis of the ETOPO1 maps for the 20 trenches, suitable segments were selected for plotting sample data on the cross-section profiles. The geoid and gravity grids were visualized based on the CryoSat-2, Jason-1 raster grids. Usage of such data as GEBCO (supported by British Oceanographic Data Centre (BODC)), ETOPO1 (from National Oceanic and Atmospheric Administration (NOAA)), [11]), SRTM (supported by National Aeronautics and Space Administration (NASA)), CryoSat-2 and Jason-1 raster grids enables to perform advanced spatial analysis that were demonstrated further in relevant sub-sections. Data analysis and processing grids were performed using GMT,

Vector data include global layers of coastlines, continents and shorelines embedded in the GMT and downloaded high-resolution American Standard Code for Information Interchange (ASCII) xyz topographic and gravity data sets from open source repositories (USGS). Coastal borders on all maps were based on the Global Self-consistent, Hierarchical, High-resolution Geography Database (GSHHGD) to visualize thematic geologic maps for the area of the trenches. High-resolution geo data collection was performed from the web-site USGS and available embedded GIS layers of GMT [740], [741], [687]. A
subset of regional geological settings was collected from the unified catalog of the Scripps Institution of Oceanography (SIO). Geologic lineaments, ridges, volcanic spots, seafloor fabric, lines of tectonic plates boundaries, rate of movements and fracture zones were mapped through the batch use of several GMT modules and converted to image files by GMT 'psconvert' module. Topographic descriptions of the submarine features were referenced to the existing IHO database of names.

<figure><img src="image_69.png" /><figcaption>**Map 3.1.1:** Comparison of the resolution in datasets: GEBCO, SRTM, ETOPO1 and ETOPO5</figcaption></figure>

Other data used in current dissertation included high-resolution bathymetric dataset GEBCO [191]; [553]. The GEBCO project aimed at high precision mapping of the ocean seafloor originated as early as origins at the beginning of the 20th century and since then up to now have been developed by several oceanographic organizations [537]; [735].

Regular updates of the GEBCO dataset grid, Research Vessel (R/V) based bathymetric observation and hydrographic meetings [170] maintain both the accuracy of GEBCO grids [705] and the actuality of
names of the submarine objects [693]. There are other data sets commonly used in the geographic studies (SRTM, ETOPO1, 5 Arc-Minute Global Relief Model of Earth's topography (ETOPO5)). However, GEBCO resolution of the 15 seconds grids gives a reliable accuracy comparing to the other raster data (Fig. 3.1.1). Using GEBCO grid was done by GMT scripting toolset using existing methodologies and technical documentation available both on GEBCO [343] and GMT [743].

<table><caption>**Table 3.1.1:** Data sources, types and precision</caption><thead><tr><th>No</th><th>Data</th><th>Origine</th><th>Type and Source</th></tr></thead><tbody><tr><td>1</td><td>ETOPO1, ETOPO5</td><td>NOAA</td><td>1 arc-min global relief model raster grid [11]</td></tr><tr><td>2</td><td>Global CMT Catalog</td><td>Harvard CMT</td><td>[154], [149]</td></tr><tr><td>3</td><td>Earth Gravitational Model of 1996 (EGM96)</td><td>National Geospatial Intelligence Agency (NGA)</td><td>described by [651]; [648]</td></tr><tr><td>4</td><td>GSHHGD</td><td>GMT embedded dataset</td><td>[740]</td></tr><tr><td>5</td><td>GEBCO</td><td>BODC, IHO Data Center for Digital Bathymetry (DCDB)</td><td>15-sec DEM raster grid [537], [553], [207]</td></tr><tr><td>6</td><td>SRTM</td><td>NASA</td><td>15-sec DEM raster grid</td></tr><tr><td>7</td><td>DEM</td><td>National Geophysical Data Center (NGDC)</td><td>NOAA</td></tr><tr><td>8</td><td>Toponimic Gazettier</td><td>IHO</td><td>[271], [272]</td></tr><tr><td>9</td><td>Global Volcanism Program (GVP) database</td><td>compiled by Smithsonian Institution</td><td>U. S.</td></tr><tr><td>10</td><td>International Seismological Centre (ISC) database</td><td>catalogue of the seismicity and earthquakes</td><td>[276]; [157]; [745]</td></tr><tr><td>11</td><td>Google Earth Pro</td><td>Keyhole, Inc.</td><td>U. S.</td></tr><tr><td>12</td><td>Geologic datasets</td><td>tectonic plate boundaries, slabs, lineaments, trenches, ophiolites, volcanoes</td><td>SIO</td></tr></tbody></table>

High precision data resources and grids were selected and chosen from the variety of the available datasets, as a ground basis of this research as summarized in the Tab. 3.1.1. General bathymetric maps were based on a cartographic visualizing of the ETOPO1 raster grid [11], enlarged visualization by SRTM [168] and GEBCO, both 15 seconds resolution. High-quality bathymetric ETOPO1 and ETOPO5 data were used and compared to perform mapping and geomorphological modelling (Fig. 3.1.1).
## 3.2 Approaches

GENERAL methodological scheme of the presented research is summarized in Tab. 3.2.1. Traditional methods of the marine geological modelling include using GIS based processing of the RS images, such as for instance aerial photos, Satellite Pour l'Observation de la Terre (SPOT)3, SPOT4 and Environmental Satellite (ENVISAT) data, or producing digital maps based on the data capture in the field [62]. The methodology of the research in submarine geomorphology is traditionally being supported by the GIS software, mainly in geological modeling. Geomorphological studies of the hydrological effects often use GIS modeling and digital terrain analysis for geospatial analysis [48].

<table><caption>**Table 3.2.1:** Software and tools used in dissertation</caption><thead><tr><td>No</td><td>Tools</td><td>Methods</td><td>Approach</td></tr></thead><tbody><tr><td>1</td><td>GMT</td><td>Cartography</td><td>GIS mapping, modelling. Profiles automatic digitizing by modules. Comparative geospatial analysis. Classification</td></tr><tr><td>2</td><td>QGIS</td><td>Cartography</td><td>GIS mapping. Profiles manual digitizing by plugins</td></tr><tr><td>3</td><td>AWK</td><td>Data conversion</td><td>Table reshaping</td></tr><tr><td>4</td><td>R</td><td>Statistical libraries</td><td>Modelling, statistical analysis, graph plotting</td></tr><tr><td>5</td><td>Python</td><td>Statistical libraries</td><td>Modelling, statistical analysis, graph plotting</td></tr><tr><td>6</td><td>Octave</td><td>Visualization</td><td>Modelling, statistical analysis, graph plotting</td></tr><tr><td>7</td><td>BIBLITEX</td><td>Bibliography</td><td>Literature review</td></tr><tr><td>8</td><td>LATEX</td><td>Compiling</td><td>Typesetting</td></tr></tbody></table>

Riverine morphological analysis [51] uses methodology based on experimental fieldwork, setup of network development combined with spatial data analysis by GIS. Other examples include application of aerial imagery of the analysis of the river deltas discharge time series, and land cover and morphological planimetric changes by methods of combined GIS based overlay of the topographic maps, air photographs, and satellite images [19], [395].

Recent methodological example of the Environmental Systems Research Institute (ESRI) ArcGIS application [40] shown riverine geomorphology covering large spatial extents by assessment of their structure, functioning and geomorphology using free datasets, GEarthView plugin for QGIS and algorithms of the RS data processing. ArcGIS based research is presented by [532] focused on a regional scale analysis of changes in channel width and pattern. SAGA GIS are suited for environmental analysis and remote sensing data processing [395], [384], [390]. The assessment of the magnitude and direction of changes caused by flood management in the low land rivers is presented by [258]. By cross-section profiles plotting large datasets of the multi-source data were analyzed: sedimentary cores, hydrographic surveys and historic map to visualize degrees of the involved engineering activities and infrastructure resulting in channel and floodplain modification that affect fluvial processes.

Combined methods of field surveys, seismic profiles echo sounding and classification of the Landsat TM imagery is a core approach [268]. The approach towards submarine canyon mapping is presented as a Remotely Operated underwater Vehicle (ROV)-based habitat mapping program using GEBCO map cartography [267]. The proposed methodology of modeling hydro-morphological degradation [328] is based on various environmental parameters of particular ecological importance (longitudinal and transverse banks, flow diversity, depth and width variance, substrate type and diversity). An application of the Matlab ProcessV3 for hydrological modeling was presented by [357]. Effective approach of the assessment of status and trends in marine fisheries and ecosystems considering marine geomorphology as proxy for marine habitats is proposed by [241].

On the contrary with given above examples, current dissertation makes an accent on using high-level
programming languages and scripting approach for processing large data frames imported from GIS.
## 3.3 Generic Mapping Tools (GMT)

MAIN toolset used in this research is GMT with modules functionality described in relevant documentation and software release announcement [743]. Comparing to the traditional GIS software, the advantages of the GMT consist in its high efficiency for the scripting based cartographic workflow rather than Graphical User Interface (GUI) as well as open source availability. General scheme of the methodology in GMT included following steps:

* Bathymetric mapping based on the extracted ETOPO1 grid by GMT modules 'grdcut' and 'grdim-age';

* mapping, visualizing and analysis of the thematic data sets: marine gravity and geoid,

* topographic and gravimetric surface modelling: based on the xyz ASCII data;

* cartographic data processing by 'psbasemap', 'psxy', 'grdcontour' for contour plotting;

* 3D-mesh modelling;

* automatically digitized set of the orthogonal cross-stacked profiles by 'grdtrack' module along the selected segment of the trench with track 400 km long, spaced 20 km, sampled 2km along;

* visualizing mean, median values and plotting trends for the profiles by 'trend1d' module;

* descriptive statistical analysis through the visualized histograms.

The use of the GMT was presented as mapping in some works, e.g. [255]. Examples of other advanced tools for geological data analysis include, e.g., PHASER diffractometer, direct observations received from the R/V cruises, observations from DEM [195]. The mapping in the presented dissertation is based on the available data sets and those embedded in the GMT [688].

GMT scripting toolset was selected as a geospatial tool for raster grids processing, mapping and automated digitizing of the geomorphological profiles. The GMT was initially developed by P. Wessel and P. Smith [743], maintained and improved to finally become a powerful open source scripting toolset for the geospatial analysis and cartographic visualization.

In contrast with classic GIS, GMT is completely based on the UNIX-based shell scripting approach with command-line usage without graphical GUI interface. Thus, it is not a classical GIS. Using the full powerful functionality of the GMT requires writing a script with numerous arguments for the fine-tuning and setting up of every cartographic element: projections, labels, north direction, major and minor grids, titles, color palettes, legends, vector elements, annotations, map elements (coastlines, directional roses, scale bars) etc.

Hence, the functionality of the GMT approach consists in its scripting when individual features for the cartographic objects can be defined in each line of the GMT code. Each small code line presents a reference to a specific GMT module. Combined together as a sequential script, the codes are used for mapping. As demonstrated further in the presented codes (see Appendix for the GMT scripts), a great advantage of the GMT is a high cartographic functionality of the toolset enabling to fully adjust and refine aesthetics of the cartographic elements: fonts, lines width and type, layout, a large variety of the visualization, plotting, modelling, advanced cartographic modules (mapping projections, modelling, grid raster data processing), color scales etc. [741].

However, the high cost of the steep learning curve of the GMT as well as its non-visual command line based operation pays off by its powerful functionality and excellent visualization. The GMT is compatible with Unix commands (cut, rm, echo) and Geospatial Data Abstraction Library (GDAL) raster
processing [393], it operates both with raster and vector data, produces high-quality maps and figures, has embedded functions for descriptive statistical analysis (histograms). All this makes it one of the most preferred data visualization and mapping toolsets by the cartographers.

<figure><img src="image_70.png" /><figcaption>**Map 3.3.1:** North Pacific Ocean: contour map in the Orthographic azimuthal projection, GMT</figcaption></figure>

As the utility of the seafloor-mapping technology, a GMT scripting was demonstrated with a case study for comparative analysis of the 20 trenches. Because the GMT is a scripting toolset consisting of modules, several main modules were used for preparing maps. The auxiliary GMT modules were used for mapping and cartographic visualizations. Methods of data visualization through GMT modules include cross-section profiles, topographic modelling, thematic geologic mapping, modelling geoid and gravity fields. A number of GMT techniques and modules (e.g. grdhisteq, pscoast, pshistogram, pstext, gmtset, gmtdefaults, grdcut, makecpt, grdimage, psscale, grdcontour, psbasemap, nearneighbor for KKT, Figure and gmtlogo, psconvert) were used to visualize, model and map seafloor characteristics, as well as texture and shape of the cross-section profiles. A 'gmtinfo' module was used to get information about data tables from the raster grids. Auxiliary GMT modules were used for technical data processing, for example, 'gmtwhich' to find the full paths to the files given on the command line, a 'gmt.conf' module was used to define configuration of the GMT parameters. Interpretation of the topographic profiles, geological sampling, gravity modelling and statistical analysis was also made by means of the GMT scripting toolset. A combination of GMT modules was applied for modelling, analysis, comparison and visualization of the variations in the geomorphology of the trenches.

ML approach of GMT was used for semi-automatic plotting, visualization, modelling, mapping and interpretation of the three hadal trenches located in southern Pacific Ocean. Sequential use of GMT modules enabled to gain further insights into their geomorphic shape, and to compare accuracy of the data grid interpretation approaches based on the scripting concept aimed to visualize cross-section transects.
Indeed, GMT cross-section stacking methodology as performed in this research 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. Indeed, using GMT enables to automatically digitize profiles, then to 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. Existing semi-automated methods of cartographic digitizing (e.g. Autotrace based) [661], [442] or sequential use of the specialized software for bathymetric mapping, Hydrosweep DS-2 [470], CARIS HIPS [443], GMT for mapping using CARIS HIPS data by grids and coastlines was driven largely by the developed algorithms of the ML. Application of the Artificial Intelligence (AI) in the cartographic routine, significantly minimizes the handmade routine and subjectivity. On the contrary, it increases the overall precision of the data processing and accuracy of visualization. The GMT proposes more advanced solutions in the automatization of the cross-sections through the ML based digitizing that results in a series of the profile transects across the 20 trenches of the Pacific Ocean.

<figure><img src="image_71.png" /><figcaption>**Map 3.3.2:** Landmask of the shorelines (land vs water areas), GMT</figcaption></figure>

Applying GMT techniques to mapping trenches enabled to model profiles and visualize geologic settings. Each trench was distinguished by specific geomorphic shape of the submarine landform, geologic substrate, complexity of the tectonic settings and historical development of the lithospheric plates movements in this area. Bathymetry data were modeled by GMT cross-section profiling for the trenches and compared. Maximal depths in the trench segments derived from digitized GEBCO and ETOPO1 grids were used for comparative analysis and visualized at histograms.
### 3.3.1 Bathymetric mapping

MAPPING oceanic trench geomorphology data from the bathymetric maps should be based on the high resolution data to highlight isolines enabling to interpret steepness and gradient which requires certain cartographic interpretation. High-resolution bathymetric data is important for mapping trench shape in a cross-section and to identify clearly its geomorphology. Pixel size at a course resolution bathymetric maps (e.g. ETOPO5) may reduce accuracy, therefore, ETOPO1 and GEBCO dataset grids were chosen for bathymetric mapping.

The ETOPO1 global relief model was modeled using `grdimage` module of GMT. Examples of trenches are presented in corresponding section (Map 2.2.1), KKT (Fig. 2.1.1). The `psbasmap`, `grd- contour` GMT modules were used to map base map elements (e.g. contours of the land, geoid and bathymetry, net grids and titles) on the bathymetric maps. Each module consists of the small code line. Combined together as a script, they were used to produce and visualize maps. For example, the GMT modules `grdimage` was main module used to map the bathymetric grid of the RT (Map 2.13.1) using GMT code:

```
gmt grdimage GEBCO_2019.nc -Cmyocean.cpt -R120/134/20/33 -JM6i -P -I+a15+ne0.75 -Xc -K > $ps
```

For example, in case of Tonga-Kermadec region, the -R command was used to cut off part of the image was by the given coordinates: -R140/195/-50/-5 from the whole grid in WESN way. As can be noticed, the southern hemisphere has negative coordinates.

The color palette was made using `makecpt` module from the available 'geo' palette adjusted according to the data range (for example, the Tonga-Kermadec area had topography from -11000 to 4500. The cartographic projection, for the bathymetric map was selected as traditional Mercator, which was set up using the following command: '-JM6i' means Mercator projection, with 6 inches width of the map; '-P' means portrait orientation; '-K' means continue of the script (not finalized). The legend on the left part of the bathymetric maps was added using 'psscale' module by following code snippet (here example is given for the Kermadec region, Map 2.5.1):

```
gmt psscale -Dg131/-50+w14.8c/0.4c+v+0.3/0i+ml Rtkt_relief.nc -J -Cmyocean.cpt -FONTLABEL=8p,Helvetica,dimgray -FONTANNOT_PRIMARY=5p,Helvetica,dimgray -Baf+l"Topographic color scale" -I0.2 -By+1m -O -K » $ps
```

For the case of Middle America Trench, the general bathymetric map (Fig. 2.3.1) was visualized using a sequence of GMT modules: 'grdcut' module for cutting area from the global grid earth_relief_01m.grd by selecting coordinates -R263/278/7/17, following by 'grdimage' module for visualizing the grid itself. Several auxiliary modules were then used both for cartographic purposes (color scale, scale bar, directional rose, grid ticks and lines) and for the map embellishment (GMT logo) using available techniques [688]. The shoreline vector layers were driven from the existing GMT data sets [740]: tectonics, bathymetry, geomorphology and geology.

The GMT modules `psconvert`, `psscale` were used for visualizing scale bar on the maps, GMT module logo to add GMT logo below the maps, `pstext` GMT module was used for adding annotations on the maps, `psxy` module - for plotting geologic lineaments and sample points. The bathymetric data are presented as plan and oblique shaded-relief grids to highlight the major and distinct seafloor landforms of the MAT (Fig. 2.3.1) which was mapped using the script provided in full in Appendix A.4.1.

Cartographic functionality of the GMT is assessed by technical visualization of the produced maps, models and graphs as well as tested variety of map projections. In this study, following various cartographic projections were used for plotting maps: Mercator (Fig. 2.19.1, 2.20.1, 3.3.7, 3.3.8, 4.21.1, 4.20.1), Mollweide pseudo-cylindrical homolographic equal-area, Eckert IV equal-area pseudo-cylindrical, Lambert Azimuthal Equal-Area (Fig. 2.19.2), Oblique Mercator (Fig. 2.20.2 by GMT, Orthographic azimuthal projection (Fig. 3.3.2), Polyconic projection, Cylindrical Equal-Area Gall-Peters, Transverse Mercator cylindrical, Fig. 2.16.1. Conic equal-area Albers projection for map with digitized lines of the
profiles (on Fig. 7. Location of the cross-sectioning 25 profiles in [419]) is plotted by GMT.

<figure><img src="image_72.png" /><figcaption>**Map 3.3.3:** Grid contour modelling by Nearest neighbor algorithm, KKT GMT</figcaption></figure>

<figure><img src="image_73.png" /><figcaption>**Map 3.3.4:** Grid contour modelling by XYZ2grd algorithm, GMT</figcaption></figure>
3.3.2 Seismic mapping

FOCAL mechanisms (Fig. 3.3.5) were drawn using 'psmeca' GMT module that reads dataset values from ASCII file and generates a PostScript code plotting focal mechanisms. Centroid moment tensor solutions are shown for shallow depth earthquakes of Mw < 10 along the PSB margins from 1976 to 2010. Earthquakes along the tectonic plate subduction zone are shown by GMT using data from Global CMT Project.

<figure><img src="image_74.png" /><figcaption>Map 3.3.5: Seismic map of the focal mechanisms: PSB margins</figcaption></figure>

Focal mechanisms depict a commonly accepted theory of the focal depths of the earthquake, showing
a depth from the surface to the earthquake's origin (hypocenter). Given dataset contains significant
amount of shallow earthquakes of the PSB with focal depths of several tens kilometers (mostly <50),
intermediate earthquakes with focal depths from 70 to 200 kilometers, and also a few earthquakes of PSB
with deep focus reaching depths > 500 kilometers, Table 3.3.1. The foci of the most PSB earthquakes
are concentrated in the crust and upper mantle, that is, originate in shallow parts of the Earth's interior.

The original dataset of the focal mechanism solutions are from the was taken from the Global CMT Catalog, formerly known as the Harvard CMT catalog [155]. Technical parameters are set up as follow-
ing:

* Moment magnitude (Mw) < 10;

* Surface wave magnitude (Ms) :< 10);

* Body wave magnitude (Mb) :< 10;

* Data time span: 1976-2010.

While magnitude is a widely understood concept, describing the energy of the earthquake release on a logarithmic scale, some technical details are presented in Table 3.3.1 and Table 3.3.2 used for mapping velocity (Fig. 3.3.6), assigning and interpreting magnitudes. Introduced in the 1930’s, the Richter Scale is the best known scale for measuring the magnitude of earthquakes [559]. L-waves (Surface waves) magnitude, creating the strongest disturbance within the upper layers of the Earth, is computed by the following algebraic Formula 3.3.2:

$$MS = \log_{1} 0(A/T) + 1.66\log_{1} 0(D) + 3.30 \quad (3.1)$$

where T is the measured wave period and D is the distance in radians.

On the contrary, earthquakes occurring deep in the Earth do not generate large surface waves. Therefore, a body wave magnitude (Mb) is scaled based on the seismic waves penetrating through the Earth’s interior (body) [559]. Body wave magnitude (mb) is measured based on the maximum amplitude A by a Formula 3.3.2:

$$Mb = \log_{1} 0(A/T) + Q(D, h) \quad (3.2)$$

where T is the measured wave period and Q is an empirical function of focal depth h and epicentral distance D [559].

Moment magnitude (MW) relies on an underlying robust physical and mathematical development through converting seismic moment using a formula calibrated to agree with MS over much of its range. Others parameters are taken as default ones.

Data are collected in both GMT ’psvelomeca’ and GMT ’psmeca’ input formats and tested. Final solution was for GMT ’psmeca’ (the reason is for the GMT version 5.4.5 compatibility). The focal mechanisms are derived from a solution of the earthquakes’ moment tensor, which is estimated by an analysis of the observed seismic wave forms.

Technically, the seismic moment tensor (Harvard CMT, with zero trace) was plotted (Fig. 3.3.5) by ’psmeca’ module using the following GMT code snippet:

```
gmt psmeca -R CMT.txt -J -Sd0.5/8/u -Gred -L0.1p -Fa/5p/it -Fepurple -Fgmagenta -Ft -W0.1p -Fz -Ewhite -O -K >> $ps
```

The scale (0.5) was adjusted to the scaling of the ’beach ball’ radius, which is proportional to the magnitude. Here the ’-Sd0.5/8/u’ implies the scale, label size and annotation placement below the ’beach ball’. Filling of the extensive quadrants was defined by -E parameter; -L parameter defines drawing the ’beach ball’ outline of 0.1 pt. Shaded compressional quadrants of the focal mechanism ’beach ball’ are colored by -Gred function. The -Fa/5p/it option was used to plot size, P_axis_symbol and T_axis_symbol, to compute and plot P and T axes with symbols (selected inverse triangle (i) and triangle (t)).

Focal mechanisms (Fig. 3.3.5) visualized as commonly accepted in geophysical mapping ’beachballs’ show a representation of the lower half of the focal sphere as viewed from above the focus at the earthquake epicenter and respective zones of compression and dilatation. Based on the recorded dataset of
<table><caption>**Table 3.3.1:** Earthquake events for 1976/1977</caption><thead><tr><th>Location code</th><th>Date</th><th>Centroid Time</th><th>Lat (°N)</th><th>Lon (°E)</th><th>Depth</th><th>HD*</th><th>CHT**</th></tr></thead><tbody><tr><td>021576A PHILIPPINE ISL</td><td>1976/2/15</td><td>1:54:30.3 GMT</td><td>13.12</td><td>126.02</td><td>15.6</td><td>4.2</td><td>7.2</td></tr><tr><td>060776A LUZON PHILIPPINE ISL</td><td>1976/6/7</td><td>7:36:0.3 GMT</td><td>14.23</td><td>125.07</td><td>24.4</td><td>4.8</td><td>4.9</td></tr><tr><td>081676B LEYTE</td><td>1976/8/16</td><td>16:11:58.7 GMT</td><td>7.07</td><td>123.75</td><td>33.0</td><td>21.7</td><td>51.4</td></tr><tr><td>081776A LEYTE PHILIPPINE ISL</td><td>1976/8/17</td><td>1:11:16.5 GMT</td><td>10.08</td><td>126.10</td><td>47.3</td><td>2.6</td><td>6.3</td></tr><tr><td>081776B MINDANAO</td><td>1976/8/17</td><td>4:19:39.8 GMT</td><td>7.14</td><td>123.01</td><td>15.0</td><td>8.7</td><td>12.5</td></tr><tr><td>092976A MINDANAO</td><td>1976/9/29</td><td>21:2:31.5 GMT</td><td>6.61</td><td>124.24</td><td>15.0</td><td>1.8</td><td>-1.2</td></tr><tr><td>110776B MINDANAO</td><td>1976/11/7</td><td>17:9:13.7 GMT</td><td>8.35</td><td>126.86</td><td>42.1</td><td>6.2</td><td>7.6</td></tr><tr><td>112276A MINDANAO</td><td>1976/11/22</td><td>4:22:20.7 GMT</td><td>7.03</td><td>123.58</td><td>60.0</td><td>3.0</td><td>-4.6</td></tr><tr><td>121476A RYUKYU ISLANDS</td><td>1976/12/14</td><td>16:6:49.2 GMT</td><td>28.12</td><td>130.64</td><td>15.0</td><td>3.8</td><td>4.8</td></tr><tr><td>010177C SOUTH OF HONSHU JAPAN</td><td>1977/1/1</td><td>11:33:45.9 GMT</td><td>30.62</td><td>136.80</td><td>476.5</td><td>1.8</td><td>4.3</td></tr><tr><td>010777B TAIWAN REGION</td><td>1977/1/7</td><td>19:36:48.6 GMT</td><td>20.78</td><td>120.01</td><td>10.0</td><td>2.5</td><td>1.7</td></tr><tr><td>010877A LUZON PHILIPPINE IS.</td><td>1977/1/8</td><td>6:41:10.9 GMT</td><td>15.77</td><td>122.73</td><td>13.3</td><td>1.6</td><td>6.8</td></tr><tr><td>011577A PHILIPPINE ISL</td><td>1977/1/15</td><td>10:49:8.8 GMT</td><td>12.80</td><td>126.06</td><td>32.9</td><td>2.2</td><td>3.0</td></tr><tr><td>011777A BONIN ISL REGION</td><td>1977/1/17</td><td>6:23:40.7 GMT</td><td>26.40</td><td>142.74</td><td>13.9</td><td>3.4</td><td>4.6</td></tr><tr><td>011977B MINDANAO</td><td>1977/1/19</td><td>13:54:9.8 GMT</td><td>5.01</td><td>126.56</td><td>50.1</td><td>4.8</td><td>5.3</td></tr><tr><td>012077A MINDANAO, PHILIPPINE IS.</td><td>1977/1/20</td><td>20:51:17.2 GMT</td><td>8.00</td><td>122.70</td><td>41.8</td><td>2.3</td><td>2.3</td></tr><tr><td>022777A MARIANA ISLANDS</td><td>1977/2/27</td><td>18:34:15.0 GMT</td><td>18.21</td><td>145.30</td><td>592.6</td><td>2.1</td><td>6.5</td></tr><tr><td>022877B MINDANAO, PHILIPPINE IS</td><td>1977/2/28</td><td>1:50:38.7 GMT</td><td>8.87</td><td>127.01</td><td>55.8</td><td>2.3</td><td>6.8</td></tr><tr><td>030277B MINDANAO, PHILIPPINE IS</td><td>1977/3/2</td><td>9:53:26.4 GMT</td><td>6.62</td><td>123.53</td><td>33.4</td><td>6.0</td><td>3.2</td></tr></tbody></table>

the P-waves, compressional longitudinal waves, a graphical model shows red areas of the 'beachballs' as upward moving, i.e. compressional motion, while white areas- a downward, dilatational motion. Determining geological information, such as rupture of the Earth's surface, it depicts two nodal planes, fault and auxiliary ones. For instance, auxiliary plane intersects the line of fault plane between the circle and its center of the 'beachball' with the general rule that the closer the intersection is to the center, the more dominant is the strike slip. Increase of the curvature of the crossing hemispheres denotes shallower dip. The strike of a nodal plane is measured in degrees around the circle from north to the line of the nodal plane. Based on existing classification [559], various focal plane solutions are shown:

1. normal fault striking north, dipping E/W at ca. 45°

2. oblique fault, right or left lateral and reverse slip

3. oblique fault, reverse slip dominant, with some right/left lateral

4. oblique fault, strike slip dominant, attitude ca. N°40W / 70°W (right/left lateral), or N35°E / 70°E (left/right lateral).

Earthquakes event map at PSB area (Fig. 2.16.3) shows prime hypocentres and magnitude values by color and size of the circles, respectively. The map was plotted on the basis of the ISC-EHB Bulletin, overlaid on the Google Earth map as a file in Keyhole Markup language Zipped (KMZ), which is a file extension for a placemark file used by Google Earth Pro. The geographic coordinates include x,y,z components in decimal degrees defined by the World Geodetic System of 1984 (WGS84). The overlay is visualized on Fig. 2.16.3. The explanation of the map is as follows. Circle size shows magnitude values. Circle color shows hypocenter depths, km (from red: 0 to dark blue/black: > 600).

A multi-source dataset applied for seismic mapping includes a compilation of GEBCO, ETOPO1, ETOPO5, CMT, EGM96, GVP and GSHHGD sources. Combined studies of seismic and tectonic settings (focal mechanism, velocity, geological lineaments), volcanic activity and general topography of the

<footnote>*Half duration; **Centroid Hypocenter Time (HcT)</footnote>
ocean seafloor are essential for understanding complicated bathymetry and geomorphic patterns of the trenches in the PSB. A multi-source GIS data analyze on the southern side of the PSB and discussion of the current situation and historical evolution of the geology of this area was performed and presented.

<table><caption>**Table 3.3.2: Magnitude and Moment Tensor**</caption><thead><tr><th>Location code</th><th>Moment Tensor</th><th>Mw*</th><th>Mb**</th><th>Ms***</th><th>SM****</th><th>FP<sub>1</sub>*****</th><th>FP<sub>2</sub>*****</th></tr></thead><tbody><tr><td>021576A Philippine Isl</td><td>Expo=25 -5.140 3.270<br>1.860 -1.560 0.510 -4.690</td><td>6.5</td><td>6.1</td><td>6.1</td><td>6.47e+25</td><td>322/39/-73</td><td>121/53/-104</td></tr><tr><td>060776A Luzon</td><td>Expo=25 -0.930 5.570<br>-4.630 -1.810 2.820 -5.740</td><td>6.5</td><td>6.1</td><td>6.4</td><td>8.35e+25</td><td>249/67/180</td><td>339/90/23</td></tr><tr><td>081676B Mindanao</td><td>Expo=27 9.330 0.240<br>-9.570 1.450 -3.340 4.330</td><td>8.0</td><td>6.4</td><td>7.9</td><td>1.09e+28</td><td>341/35/92</td><td>158/55/-89</td></tr><tr><td>081776A Leyte Philippines</td><td>Expo=25 -1.250 0.390<br>0.860 -0.500 0.070 -0.570</td><td>6.0</td><td>6.0</td><td>5.7</td><td>1.34e+25</td><td>340/39/67</td><td>131/55/-108</td></tr><tr><td>081776B Mindanao</td><td>Expo=26 -1.220 5.310<br>-4.100 -1.920 1.330 2.110</td><td>7.1</td><td>6.2</td><td>6.8</td><td>5.71e+26</td><td>214/64/-172</td><td>120/83/-26</td></tr><tr><td>081776B Mindanao</td><td>Expo=24 -4.620 2.220<br>2.400 -0.530 -1.810 1.870</td><td>5.7</td><td>6.0</td><td>5.4</td><td>4.78e+24</td><td>232/35/-76</td><td>34/56/-100</td></tr><tr><td>110776B Mindanao</td><td>Expo=26 1.540 0.230<br>-1.770 0.990 0.610 -0.170</td><td>6.8</td><td>6.0</td><td>6.8</td><td>2.03e+26</td><td>162/39/49</td><td>30/62/118</td></tr><tr><td>112276A Mindanao</td><td>Expo=23 2.030 4.260<br>-6.290 1.780 2.020 0.220</td><td>5.1</td><td>6.0</td><td>0.0</td><td>6.07e+23</td><td>136/59/14</td><td>38/78/148</td></tr><tr><td>121476A Ryukyu Isl</td><td>Expo=25 -4.160 2.940<br>1.230 -2.490 0.330 1.350</td><td>6.4</td><td>6.3</td><td>6.2</td><td>4.62e+25</td><td>230/33/-114</td><td>78/61/-75</td></tr><tr><td>010177C Honshu Japan</td><td>Expo=24 -0.324 0.799<br>-0.475 1.014 -0.356 0.401</td><td>5.4</td><td>5.2</td><td>0.0</td><td>1.34e+24</td><td>33/32/-163</td><td>289/81/-59</td></tr><tr><td>010777B Taiwan Region</td><td>Expo=24 -3.009 -0.044<br>3.054 1.185 -2.128 -0.023</td><td>5.7</td><td>5.7</td><td>5.1</td><td>3.88e+24</td><td>194/29/-62</td><td>343/65/-105</td></tr><tr><td>010877A Luzon Philippines</td><td>Expo=23 6.364 2.443<br>-8.807 1.626 -1.176 2.219</td><td>5.2</td><td>5.3</td><td>0.0</td><td>8.16e+23</td><td>8/43/118</td><td>152/53/67</td></tr><tr><td>011577A Philippine Isl Reg</td><td>Expo=24 -1.842 1.248<br>0.594 -0.737 0.570 -1.652</td><td>5.5</td><td>5.6</td><td>4.6</td><td>2.41e+24</td><td>311/34/-88</td><td>128/56/-91</td></tr><tr><td>011777A Bonin Isl Reg</td><td>Expo=24 4.813 0.539<br>-5.352 0.757 4.521 0.467</td><td>5.8</td><td>5.6</td><td>5.6</td><td>6.87e+24</td><td>164/25/71</td><td>5/66/99</td></tr><tr><td>011977B Mindanao</td><td>Expo=25 1.646 -0.866<br>-0.781 -0.990 0.267 1.834</td><td>6.2</td><td>5.8</td><td>5.9</td><td>2.45e+25</td><td>161/40/130</td><td>293/61/62</td></tr><tr><td>012077A Mindanao</td><td>Expo=24 0.044 1.453<br>-1.496 0.517 0.947 1.691</td><td>5.5</td><td>5.4</td><td>4.9</td><td>2.49e+24</td><td>111/65/3</td><td>20/87/155</td></tr><tr><td>022777A Mariana Islands</td><td>Expo=24 -1.502 -0.378<br>1.880 -0.671 -0.473 -0.742</td><td>5.5</td><td>5.0</td><td>0.0</td><td>2.04e+24</td><td>136/48/-132</td><td>10/56/-53</td></tr><tr><td>022877B Mindanao</td><td>Expo=24 2.750 -1.458<br>-1.292 0.444 -0.274 0.678</td><td>5.5</td><td>5.7</td><td>0.0</td><td>2.46e+24</td><td>313/39/92</td><td>130/51/88</td></tr><tr><td>030277B Mindanao</td><td>Expo=25 3.195 -3.570<br>0.375 0.728 -1.526 2.261</td><td>6.4</td><td>6.1</td><td>6.1</td><td>4.34e+25</td><td>271/41/53</td><td>135/58/117</td></tr></tbody></table>

<footnote>*Moment magnitude (Mw); ** Body wave magnitude (Mb); ***Surface wave magnitude (Ms); ****SM: Scalar Moment; *****Fault Plane (FP); ******strike/dip/slip.</footnote>
### 3.3.3 Velocity mapping

VELOCITY ellipses with 95% confidence in rotated convention (red) and in (N,E) convention (green), on a PSB region were plotted using 'psvelo' module of GMT, Fig. 3.3.6. Theoretical background of the best-fitting angular velocities and Mid-Ocean Ridge VELocity (MORVEL) is presented by [128] who described and modeled geologically current motions of 25 tectonic plates using geologically determined and geodetically constrained subsets of the global circuit. The datasets include Pacific VELocity (PVEL)

<figure><img src="image_75.png" /><figcaption>**Map 3.3.6:** Velocity ellipses at PSB margins</figcaption></figure>

Velocity ellipses at PSB margins from 1976 to 2010 in (N, E) convention (green) and rotated convention (red). Rotational wedges (purple 'gear wheels'): by 'psvelo' module of GMT. Data: Global CMT Project. Contour: GEBCO grid. Technical visualization of the velocity map (Fig. 3.3.6) was performed
using following GMT code snippet:

* Velocity ellipses in (N,E) convention: `rgmt psvelo CMT.txt -R -J -W0.3p,green -L -Ggreen -Se0.1/0.95/5 -A0.2p -O -K >> $ps'`. Here, the `Vscale` gives the scaling of the velocity in inches; the `CMT.txt` is the dataset; parameters `-R` and `-J` indicate the region (in this case coordinates: `-R120/152/4/35`)

* Velocity ellipses in rotated convention: `rgmt psvelo CMT.txt -R -J -W0.3p,red -Sr0.1/0.95/5 -Ggreen -A0.2p -O -K >> $ps'`. Here the `Sr0.1/0.95/5` is the main Keyword arguments (`kwargs`) showing the `Vscale` (0.1 inch), `Confidence` (0.95) and `Fontsize` (5pt).

* Rotational wedges: `rgmt psvelo CMT.txt -R -J -W0.05p -L -Sw0.5/1.e7 -D2 -Gslateblue2 -Elightgray -A0.2p -O -K >> $ps'`. Here, the `-Sw0.5/1.e7` indicates Rotational wedges with `kwargs` as `Wedge_scale` and `Wedge_mag`. The `-D2` parameter enables to rescale the uncertainties of velocities by `sigma_scale` (in this case 2 scale was set up).

The respective red and green circles with confidence ellipses are PVEL dataset estimates, which use a plate circuit to estimate subduction across the plate subduction zones in the west Pacific. According to [129], relative to the PSP, the Pacific Plate rotates conter-clockwise around a pole near the southern end of the plate boundary. The best-fitting angular velocity indicates rapid decrease of the convergence rates southward, from $49 \pm 0.7 \text{ mmyr}^{-1}(1\sigma)$ at the northern end of the IBT to $9 \pm 0.8 \text{ mmyr}^{-1}$ of orthogonal subduction along the southern YT, south off MT [128].

The PSP-Pacific plate velocities are estimated to fit the interval $1 \text{ mmyr}^{-1}$ and $2^\circ$ along the lithospheric plate boundary [668], which reflects the overlap in the GPS stations used to estimate the motions of the two plates. The existence of the Caroline plate was first proposed by [738] who estimated its motions from a synthesis of marine seismic, bathymetric, and seismologic observations from its boundaries with the Pacific and Philippine plates. Located in the western equatorial Pacific immediately south of the PSP (Fig. 2.11.2), it remains the poorly understood and enigmatic with uncertainties about the style and rate of its present deformation caused by scarcity of reliable kinematic data.
### 3.3.4 Histogram equalization

THE smoothness of the grids intensities for GEBCO and ETOPO1 was improved by passing the output of 'grdgradient' to 'grdhisteq' with a case example of the Hikurangi, Puysegur and Hjort trenches [387].

New Zealand region, Hikurangi, Puysegur and Hjort trenches:
Histogram equalization on topography grid
ETOPO1 DEM Global Relief Model 1 arc min resolution

<figure><img src="image_76.png" /><figcaption>**Map 3.3.7:** Histogram equalization on the topography grids, ETOPO1</figcaption></figure>
The histogram equalization was performed by key GMT module 'grdhisteq' to enhance visualization of the various ranges of the comparison of two grids (ETOPO1 and GEBCO) depending on their frequency distribution (Fig. 3.3.7 and Fig. 3.3.8). Four subplots were generated to find values dividing grid files into equal area patches. The grdhisteq module performs a histogram equalization of raster images. Using 'grdhisteq', the ASCII data values dividing the range of the initial raster data into cells segments was written to an output file. A -C16 argument defines 16 cells. The resulting raster grid has an equal area in the images. Using 'makecpt' GMT module this output was then colored according to the actual bathymetric values by following code: gmt makecpt -Crainbow -T- 9000/3000 > t.cpt. The explanations for the subplot with four maps (Fig. 3.3.7 for ETOPO1 and Fig. 3.3.8 for GEBCO) is as follows:

1. Top left: The original map (top left) visualizes initial raster (artificial cpt rainbow) to highlight changes in elevation;

2. Top right: equalized raster grid;

3. Low left: normalized raster grid derived using -N argument that stands for Gaussian output, used to to receive data with smooth Gaussian distribution. The default standard normal scores was used for grid normalization.

4. Low right: quadratic equalization

The subplot on the low right of the presented Figures 3.3.7 and 3.3.8 shows quadratic equalization plotted using following snippet: gmt grdhisteq hpt_relief.nc -Gout.nc -Q. Here, the -Q argument is output selecting quadratic histogram equalization, unlike the default linear one. The -T-9000/3000 argument means range of the topographic elevations. The image was then visualized using 'grdimage' module by code snippet:

$gmt grdimage hpt_relief.nc

-I+a45+nt1 -Ct.cpt -JM3i -Y6i -K -P -Bpxg5f5a10 -Bpyg4f2a8 -Bsxg5 -Bsyg4 -BWSne > $ps

In this code, the -Ct.cpt argument passes cpt created in the previous step for visualization; the '-JM3i' arguments explains the Mercator projection with 3 inches width; '-Y6i' argument plots the map with 6 inches distance from the previous one by Y axis; '-K' argument refers to the continuation of the code. The annotation was added using Unix 'echo' prog: echo 172 - 33Original"|gmp pstext -Rhpt_relief.nc -J-F +jBL+f12p-T-Gwhite@10-Dj0.1i-O-K >> $ps.

The 'grdhisteq' GMT module enables to write a raster grid with statistics based on cumulative distribution function. That means, after applying the 'grdhisteq' module, the output raster file has relative elevations in the same locations as the input file (that is, x, y). However, the values are modified to reflect their place in cumulative distribution with reference to the initial input file. This illustrates the principle of the equalization of the topographic grids by GMT module 'grdhisteq', as visualized on Fig. 3.3.7 for ETOPO1 and Fig. 3.3.8 for GEBCO, respectively. The difference in grids intensities for GEBCO and ETOPO1 raster grids can be clearly seen on the subplot 4 on quadratic equalization. The grid was smoothed prior to use by applying 'blockmean' module of the GMT: it interpolated the original values in the table by smoothing algorithm. Two other modules were tested as well: 'blockmedian' and 'blockmode' as pre-processors before running surface to avoid aliasing in a grid.

The effect of using different histogram equalization algorithms was examined using subplots for data on GEBCO and ETOPO1 covering study area for three trenches (Fig. 3.3.7 and Fig. 3.3.8). Equalization compression for bathymetry grids were tested as comparison for four grids: original, equalized, normalized, quadratic algorithms. The results for this step are demonstrated in Fig. 3.3.7 and Fig. 3.3.8 showing a visual summary of a grid equalization by different GMT algorithms. Consistent cpt across all raster dataset grids where applied for better comparison. Differences can be seen between the values calculated at quadratic equalization at ETOPO1 and GEBCO. From the analysis of histogram equalization
algorithms it is notable that there are variations in the grids values obtained depending on the method: equalized, normalized, quadratic algorithms. This in itself means that the topic is worthy of further elaboration. Through the demonstrated examples (Fig. 3.3.7 and Fig. 3.3.8) the effects of using different algorithms varying with raster dataset grids equalization can be notable.

<figure><img src="image_77.png" /><figcaption>Map 3.3.8: Histogram equalization on the topography grids, GEBCO</figcaption></figure>
### 3.3.5 Geological mapping

GEOLOGICAL mapping (example of the KKT shown on Map 3.3.9) is based on the use of the sequence of the GMT modules, the most important of which are the following. Example of the geological mapping of the hadal trench using QGIS with a case study of the MT, the deepest trench on the Earth includes application of R programming and various statistical approaches towards data modelling. The scope of the research is more focused on the geological assessment than statistical analysis. Therefore, the special attention was paid to mapping geological lineaments, faults, trench extent, location of the Large Igneous Provinces (LIPs), volcanic spots and earthquakes in the KKT (Map 3.3.9). Various depth of the earthquakes are presented in different colors, as well as fracture zones, ophiolites, magnetic anomalies and extension of the tectonic slabs (Map 3.3.9).

<figure><img src="image_78.png" /><figcaption>**Map 3.3.9:** Kuril-Kamchatka Trench: map of geologic settings, [409]</figcaption></figure>

The geological layers (lines and points) were added on the thematic geological maps. Plotting of the vector elements on the map is done by `-W` command, e.g. `-Wthinnest,red`. Each code is saved to the initially created file using »ps command. The 1st mention of this command goes with single bracket, e.g. `>` $ps'. Then all the elements added to the map overlay the same. In this sense, there is certain similarity with GIS layers and a sequence of codes in the GMT script. The initial files with extension gmt (e.g. ridge.gmt) are tables with attribute data (coordinates and values) in native GMT format. Example of the geological mappings for the Kermadec-Tonga region is on Map 2.6.1, right, KKT (Map 3.3.9). GMT module psxy was used for adding annotations on the maps, plotting geologic lineaments and sample points. The psbasemap and grdcontour GMT modules were used to map base map elements: bathymetry contours, net grids, titles and subtitles. The full GMT codes for geological mapping are presented in Appendix, for contour mapping A.4.3, regional mapping A.4.6 and global geologic map A.4.5.
### 3.3.6 Geoid modelling

THE geoid model is a complex function that depends on various geophysical phenomena, such as permanent Earth tide, postglacial land uplift, sea level changes, polar drift and other parameters [151], [76], [236]. The geoid model was analyzed to determine the oceanographic properties of the study area: a mathematical shape that of the ocean surface on the Earth with the sole impact of the gravity and rotation, with absent other influence (e.g. winds and tides). A modeled geoid presents a smooth yet irregular surface with visualized shape resulted from the uneven distribution of the mass on the Earth's surface.

The sequence of GMT modules (gmtset, grd2cpt, grdimage, pscoast, grdcontour, psbasemap, psscale, psimage, logo, pstext, psconvert) was used to map the geoid model for the trench by GMT code areas for plotting contours of the terrestrial and water areas, geoid and bathymetry, net grids and basic cartographic elements such as titles, scales, annotations. Each GMT module consists of the set of code lines, which, combined together as a script produce and visualize maps.

The geoid modelling of the RT (Map 2.13.3) was done using following code: gmt grdimage geoid.egm96.grd -I+a45+nt1 -R120/134/20/33 -JD127/24/27/30/6i -Cgeoid.cpt -P -K > $ps Example of modelling geoid of MAT area is on Map 3.3.10. The additional cartographic elements were added using GMT basemap module: grid, title, coastline and GMT pscoast module. The full script is presented in Appendix A.4.4.

<figure><img src="image_79.png" /><figcaption>Map 3.3.10: Middle America Trench: geoid model</figcaption></figure>

There are various existing gravity and geoid concepts that originate from differing numerical applications and approaches in processing of the permanent Earths tide [152]. Modelling geoid and gravity was plotted using sequence of several GMT modules (to be more specific: img2grd, grd2cpt, grdimage, pscoast, makecpt, grdcontour and psbasemap) based on the existing methodology [420], [418]. The data are global marine gravity model from CryoSat-2, Jason-1, free-air gravity grids, Scripps Institution of Oceanography [649]: grav/grav_27.1.img and grav/curv_27.1.img via topex.ucsd.edu.

The geoid data grids used for current study and retrieved from the website geoid data grids are based on the EGM96 with n=m=360 and WGS 84 reference system [528]: geoid.egm96.grd. The EGM96 has reference ellipsoid the same as WGS 84, but higher spatial resolution and more accurate geoid. The
model is based on the fully-normalized, unitless spherical harmonic coefficients. The data are made public by National Imagery and Mapping Agency (NIMA), the NASA Goddard Space Flight Center (GSFC), and the Ohio State University.

<figure><img src="image_80.png" /><figcaption>**Map 3.3.11:** Mariana Trench: color geoid image with grey-shaded topography of the clipped areas, GMT</figcaption></figure>

The EGM96 geopotential model is a spherical harmonic model of the Earths gravitational potential complete to degree and order 360 (Lemoine et al., 1998). The EGM96 contains fully-normalized, unitless spherical harmonic coefficients and their standard deviations for a gravitational model from 2°, to 360°. These coefficients are consistent with scaling values of gravitational models [613]. The data are based on the surface gravity data covering various regions of the globe including southern hemisphere.

### 3.3.7 3D mesh grid modelling

THE 3D mesh gridding aimed to visualize and highlight the topographic and bathymetric landforms and patterns of the region of the trenches, through dividing the input grid model into smaller shapes formed after discretization of the geometric model of the area.
Composite overlay of the 3D topographic mesh model
on top of the 2D geoid contour plot

Region: Guatemala Trench

Perspective view azimuth: 165/30°

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

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

2019 Jul 07 19:07:47

Data: 2 min World Geoid Image 9.2, ETOPO 5 arc min grid

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

THE GENERIC MAPPING TOOLS

<caption>**Map 3.3.12: Middle America Trench: 3D topographic mesh model map**</caption>

Hence, the purpose of the 3D-mesh gridding (example for the MT: Map 3.3.12) was to create an overlay that highlights the geomorphic surface of each of the trenches and adjusting continental steep slope. Therefore, unlike in case of bathymetric mapping (Fig. 2.3.1) where the data were taken as (ETOPO1, lower resolution grid was used (ETOPO5). This was done to avoid small-scale irregularities possible in case of the high-resolution data which tend to visually overload the derived topographic attributes, making 3D mesh overlay visually complex. This allowed a smoother visualization of the 3D-mesh polygons by removing unnecessary details in 3D plotting using GMT code.

The 3D modelling consists in the following steps. the data were initially extracted from the ETOPO5 grid using the following ( snippet for the Vanuatu ) [407] : gmt grdcut earth_relief_05m.grd -R160/170/- 25/-10 -Gvv3d_relief.nc. The ETOPO5 was selected to plot the 3D with aim to reduce the grid size and hence to save computer memory space. The next step of the 3D modelling is a calculation of the regular 3D grid mesh using GMT parameters of the 'grdview' module. The dimension of the grid was chosen in a way which ensures that every point of the trench axis is located within the grid boundaries. For example, the coordinates for Vityaz were selected as 163°173°-15°0°and for Vanuatu [407]: 160°170°-25°-10° . Visualization parameters of the 3D models can be summarized as follows. Grid extent was in 160/45 azimuthal direction rotation with the grid ranges from 160°to 173°for both grids, in latitudinal direction
from -25°to 0°, and in vertical direction from -9,000 km depths up to 1000 km. On 'Z' argument the grid
was extended by 4cm as a visual distance between the 3D and basemap 2D, to enable interpreted them
better. The 2D mesh was divided into small simple polygons with curved contour lines and color palette
highlighting the topography. The 3D map was then overlaid with 2D geoid.

The background geoid contour (a plane grayscale map below the map, Map 3.3.12) was plotted using GMT module 'grdcontour'. The geoid field data of the 2D map were taken from the ETOPO1 dataset are given by following code (here: Vanuatu) [407]: `gmtgrdcontoureharth_relief_01m.grd -JM6.0c -R160/170/-25/-10 -p160/45 -C500 -Gd6c -Wthinnest, dimgray -Y5.5c -P -K > $ps`. In order to ensure that grids are visually depict geomorphic shapes of the surfaces, the '-Qsm' argument was selected to specify 's' for the surface plot, and 'm' to have mesh lines drawn on top of surface. For the vertical component, an option '-N-9000+glightgray' was selected to draw a plane at the -9,000 m level of bathymetry. Final 3D model of the two subplots is shown on Fig. 4.9.2. The contour lines are given by every one with annotations at every five (-C1-A5). The rotation of the 3D mesh grid was chosen at the angle most corresponding to the visualization of the trenches, after several trials at different angles, depending on the oblique rotation and depth by GMT code. A plane was drawn at -7500 m s z-level with dimgray color provided via the +g modifier to highlight the frontal facade between the 3D plane and the data perimeter and -Wf0.1p,red setting used for the pen outline.

After the data were compared and optimal resolution selected, the rotation azimuth was set as -p165/30° with an overlay distance from the geoid contour at 4,5 cm. As a result of the combination of several GMT modules, 3D mesh model of the trenches (example of the Aleutian: Fig. 2.2.4) was computed and visualized. They show the 3D model rotated in 165/45°overlapping the contour map of the geoid. Several GMT modules were used, e.g. 'pscoast' module used for mapping coastal lines and contour, 'psscale' module was used for mapping color legend strip showing visual colors, 'pstext' module was used to add textual explanations on the map. The Network Common Data Form (NetCDF) grid format [618] was processed as a grid file. The cartographic routine include adding elements on the map: plotting ticks, selecting lines for grids, adding directional rose and legend, selecting proper projection. Topographic data were collected from the ETOPO5 covering seafloor bathymetry and coastal topography of the area of each trench. The modelling was performed using GMT module 'grdview' followed by a sequence of the GMT modules: grdcut, grd2cpt, grdcontour, pscoast, grdview, logo, psconvert. Additionally, some UNIX progs (rm, echo) were used. The output 3D model is presented at (Map 3.3.12). The full GMT code is available in Appendix, A.4.2.

3.3.8 Gravity modelling

GRAVITY model is used to support geophysical investigations, geological prognosis as well as geological resources management. Here the GMT script is presented for modelling gravity from m the 'grdimage' GMT module, and 'img2grd' module was used for the data conversion. The combination between the GMT modules and cartographic features enables to properly visualize thematic map by the GMT tools. The essential GMT elements include selecting map projections, fonts, map layout, grid, adding elements on the map. Gravity data showing free-air gravity anomalies covering study area including terrestrial and marine regions, were taken from the data sets based on the 1 arc-minute resolution CryoSat-2 and Jason-1 satellite-derived raster grids described in detailed in relevant studies [651], [648], [687], [649], for the trenches and surrounding regions. Gravity data revealed ground information for the geologic and geophysical settings of the study area, as free-air gravity anomalies reflect subsurface density variations [141]. The gravity map was plotted using GMT following code (here example for the RT):

`gmt grdimage gravRT.grd -I+a45+nt1 -R105/123/8/24 -JY114/12/6.5i -CgravRT.cpt -P -K > $ps`

gmtdrdimage gravRT.grd -I+a45+nt1 -R105/123/8/24 -JY114/12/6.5i -CgravRT.cpt -P -K > $ps
<figure><img src="image_84.png" /><figcaption>**Map 3.3.13: Middle America Trench: gravity model map**</figcaption></figure>

Here the -R function shows the selected region, the JY114/12/6.5i stands for the cylindrical equal-area Gall-Peters projection and 6,5 inches the overall dimensions of the map; -I+a45+nt1 shows the illumination parameters and -CgravRT.cpt selects the color palette created on the previous step.

Chosen color palette (gravRT.cpt) shows the transition from lower (red) values (loser than 20 mGal) to higher values (bright blue) gravity values (over 90 mGal). The gravity modelling (example: Fig. 3.3.17) was performed using the following script. The example of the RT is on Map 2.13.4, MAT: Map 3.3.13 and Map 2.3.4. The full GMT code of modelling gravity is available in Appendix, A.4.7.

Gravity and bathymetric data were examined to obtain a visualization and better understanding of the seafloor fabric and structures within the PSB. The EGM96 gravity grid was visualized based on the with the satellite altimetry data [688] producing gravity map of the PSB (Fig. 2.11.4). The subset of the img file was extracted through the following GMT shell script (example for the Kuril-Kamchatka Trench): Appendix, A.4.7, (Vertical gravity for Kuril-Kamchatka Trench).
3.3.9 Automatic digitizing of the cross-section transecting profiles

CROSS-SECTION profiling is a common technique used in geosciences to receive a sample of the selected segment in the study area and analyze variability of the studied elements along the track. Examples of the cross-section area presented: DEM cross-profile [772]. Semi-automated methods have been explored in cartography, driven largely by the advantages of using machine-learning algorithms in cartographic routine, thus minimising handmade routine, subjectivity in processes of the digitizing and vectorizing of the bathymetric data.

<figure><img src="image_85.png" /><figcaption>Map 3.3.14: Izu-Bonin Trench: Digitizing cross-section profiles</figcaption></figure>

Example of such works is presented by [661]. The GMT proposes more advanced solution in the automatization of the cross-sections through ML based digitizing that results in a series of the profile transects across the trench [739]. Two methods were tested in visualizing the profile by GMT arguments: median '-Sm' and the mean '-Sa'. Of these two, median was selected as the best one representing the empirical nature of the trench's shape. Although mean method shows fine interpolation accuracy, median is more favorable, since visual representation of the trench lies within a reasonable reflection of the submarine geomorphology.

The geomorphological modelling consisted in plotting sets of the cross-section orthogonal profiles in
perpendicular direction to the main trench axis, and then comparing their parameters for each of the 20
trenches. Plotting cross-sectional profiles, or transects along the selected segment of the study area is a
common technique in geological investigations. Examples of the geological cross-section profiles plot-
ting and data analysis can be seen in [639], [556]. In this study, the orthogonal cross-section bathymetric
profiles were drawn automatically using presented GMT script across the trenches. Because the scope of
the current work is rather large (20 trenches of the Pacific Ocean), the automatization of the process was
crucial point in the methodology and a key decision for selecting GMT.

To analyze the geomorphology of the trenches, the method developed by [741] was adapted in this
study where the original focus was set on the modelling terrains, since the oceanic trenches have a similar
geomorphic setting with reverse general shape form. In this method, the observed elevations, Z, along
the cross-section transects is read into a table, using GMT 'grdtrack' module along the selected segment
of the trench with track 400 km long, spaced 20 km and samples taken 2km along the corresponding line
[420]; [418]. These recorded bathymetric measurements are then read into the input table from which a
graph was plotted denoting the profile shapes (example for the Philippine and MT: Fig. 4.12.1 and Fig.
4.17.1).

The profiles have length of 400 km, that is from -200 to +200 km from the trench axis, spaced 20 km, sampled every 2 km to achieve evenness of the sampling for the trenches. The depths values were calculated by each profile, and their mean and median values were modeled thereafter. Mean values are notable for smoothing the splines rather as a mathematical functions which lesser reflects the geographic phenomena of the objects. The uncertainties of the actual geomorphology are expected to exceed the theoretical interpolation of the actual geomorphology represented on the graphs. Hence the shape of the trenches is highlighted as median using -Sm argument, example is visualized by red line on Fig. 4.21.1.

The combination of several GMT modules was used for mapping geomorphological profiles, as shown
in GMT. Plotting and stacking of the profiles was done using using 'grdtrack' module of GMT by sta-
tistical mean, through the following main lines of the code based on the 'tkt_relief.nc' raster in NetCDF
format. Minor GMT modules used for plotting cartographic elements included 'psxy', 'gmtlogo', 'pscon-
vert', as well as Unix 'echo' utility. For example, with a case of the Philippine Trench, the target area for
the was selected off Mindanao Island, the 2nd largest island in the Philippines, situated in the southern
part of the archipelago. The Area of Interest (AOI) square is 126.62°E 10.7°N 127.27°E 8.0°N. The AOI
for the MT is located in the southwestern part of the crescent (Fig. 2.11.1) and is bounded on the east
by the Guam Island and on the south by the YT (Fig. 2.11.2) encircled by the target square of 141.2°E
11.0°E 143.8°E 11.65°N, capturing the Challenger Basin.

A sequence of the GMT codes was used for plotting transects as follows (example for Puysegur). A
'cat' Unix 'cat' utility was used to plot target coordinates. The line of transect with start and end points
was plotted using 'psxy' module. Generating cross-track profiles 400 km long, spaced 20 km, sampled
every 2 km and saving them in a table (here: tableP.txt) was done using GMT code (Appendix A.4.8).
Other examples of modelling and digitizing a series of the cross-section profiles along the two trenches:
Hikurangi and Puysegur is presented on Fig. 2.19.2, performed using GMT and to plot the graphs using
retrieved tables (Fig. 2.20.2) [387]. Example is presented on Fig. 2.6.1, left. The trench segment and
end points were plotted and profiles written into table with modelling done at the next step.

The annotations were added using Unix utility 'echo'. The advantage of the GMT for the automated
profiles digitizing lies in the increased precision and speed of the plotting and highlighted submarine
geomorphic features of the trench. For the case of MAT (Map 4.4.1), the cross-section profiles were
taken through the key sea floor features at coordinates 268.4°W 13.0°N 271.9°W 11.6°N to show the
variety of the steepness angles. Full GMT Bourne Again SHell (bash) code and characterization of
the module options and elements in comments including UNIX progs (echo, rm, cat) are provided in
Appendix, A.4.8.
### 3.3.10 Modelling 2D statistical histograms

STATISTICAL data modelling in geology is an important tool for the data analysis [125], [116]. The statistical techniques of the data analysis are widely used in geospatial studies, applied and described in previous works: [563]; [98]; [645]; [633]. Proper data interpretation always include testing several algorithms, selecting appropriate approaches and statistical packages [299]. This thesis applied the histogram techniques showing data distribution (depths), along the profiles. Unlike traditional statistical data analysis that utilizes standard tools for this purpose, such as R and Python [55]; [196], MATLAB/Octave, Statistical Package for the Social Sciences (SPSS) IBM Statistics [645], [425], Gretl [426], this research demonstrated the possibility to perform statistics directly from the cartographic GMT toolset using embedded statistical modules.

<figure><img src="image_86.png" /><figcaption>**Figure 3.3.1:** Kuril-Kamchatka Trench: 2D histograms and rose diagrams of depths [409]</figcaption></figure>

Plotting statistical 2D histograms and rose diagrams was done using GMT code by 'pshistogram' and 'psrose' modules, respectively (KKT: Fig. 3.3.1). It shows statistical 2D histograms and rose diagrams for comparative analysis of depths distribution in the southern and northern part of the KKT. Rose diagram shows the aspect of the slopes. Modelling digitized cross-section profiles was followed by a comparison of the two sets either in southern and northern segments of one trench, or comparison of the two neighboring trenches. The statistical modelling and digitizing cross-section profiles were done using GMT modules pshistogram; and 'grdtrack'. The statistics by histograms was performed using a sequence of the GMT modules: gmtset, psrose, pshistogram, pslegend, logo and psconvert. The histogram plotting demonstrated the frequency of the data distribution against bathymetric variables. The full GMT code of modelling 2D statistical histograms and rose diagrams is presented in in Appendix, A.4.12.
### 3.3.11 Modelling 3D statistical histograms

EXAMPLE of the 3D modelling of the histograms is presented on Fig. 3.3.2 (here: the KKT). The full GMT code of modelling 3D statistical histograms is presented in Appendix, A.4.13.

<figure><img src="image_87.png" /><figcaption>Figure 3.3.2: Kuril-Kamchatka Trench and Sea of Okhotsk: 3D histograms: depths distribution [409]</figcaption></figure>

### 3.3.12 Modelling statistical trends of the profiles gradient

ARIOUS approaches in modelling trends are presented to show the variety of the mathematical algorithms embedded in the GMT tool set. The trend curves were modeled using trend1d module by the following code: gmt trend1d -Fxm stackRTn.txt -Np1 > model.txt Using generated table the statistical data processing were performed: modeling trend regression for the cross-section profiles (Fig. 3.3.3), computing trends for the profiles and plotting two types of the statistical histograms, that is rose diagram and standard box histogram (Fig. 4.4.1). A sequence of the codes was applied to model of the mathematical approximation of the slope gradient trends by profiles (example is visualized on Fig. 4.7.1). For example, the basic LS line shown on the Fig. 4.7.1, E (lower left) was plotted using following methodology combining GMT and Unix utilities. Four types for the mathematical models were used to analyze the curvature of the profiles, GMT regression analysis Formula 3.3.12:

$$m_2t = a + b * t \qquad (3.3)$$

GMT regression analysis Formula 3.3.12:

$$m_3t = a + b * t + c * t^2 \qquad (3.4)$$
<figure><img src="image_88.png" /><figcaption>Figure 3.3.3: Middle America Trench: trend regression models for the cross-section profiles</figcaption></figure>

GMT regression analysis Formula 3.3.12:

$$m_5 t = a + b * t + c * t^2 + d * \cos 2\pi * t + e * \sin 2 * \pi * t \quad (3.5)$$

GMT regression analysis Formula 3.3.12:

$$\varepsilon(t) = y(t) - m_5 t \quad (3.6)$$
Example of the resulting graph of the trend regression models for the MAT was plotted on Fig. 3.3.3. The GMT scripts used for the cartographic visualization and modelling of the statistical trends are based on GMT codes provided in Appendix, A.4.11.

### 3.3.13 Surface modelling

THE dataset of this research was augmented by free ASCII data on gravity and topography as xyz table, collected from the available source in a table format: USGS. The data range was set to the following coordinates square (here with a systems of 0-360; 0-90 rotation): Lon 263°/278°; Lat 7°/17°. According to the gravity data inspection (gmt info grav Mat.xyz), the range for the data is as follows: total number of points = 554115; minimal value is -157.1 mGal at point with coordinates 263.0083° 7.0076°; maximal value is 452.4 mGal at coordinate 278.0083° W/ 17.0032° N. Based on this, the color palette was adjusted accordingly to the data range: gmt makecpt -Ccyclic.cpt -V -T-160/455 > surface.cpt where -T-160/455 shows the min/max data. Thereafter, before modelling surface, the data were interpolated for the correction for the source data by generated 1 by 1 minute block for the mean values from the raw ASCII xyg table. Afterwards the data were processed by 'gmt surface grav MAT BM xyg -R263/278/7/17 -T0.25 -I30s -GSurface MAT nc -Vv' command, which modeled the surface from the xyz data. Both topographic surface modelling (Fig. 2.3.5) and gravity surface mapping (Fig. 3.3.15) are based on the modelling numerical ASCII data in xyz format to the NetCDF format [618].

<figure><img src="image_89.png" /><figcaption>Map 3.3.15: Middle America Trench: gravity surface map modeled by ASCII data</figcaption></figure>

Gradient illumination was added to intensify the surface with azimuth 45° using gmt 'grdgradient' module with shading azimuth gradient 45°. It highlighted topographic pits, furrows and submarine landforms on the seafloor of the MAT. The data were coded by the cyclic color palette (cyclic.cpt) with the resulting image: Fig. 3.3.15. Raster image output was made by GMT code 'gmt grdimage Surface MAT nc -Csurface.cpt -R263/278/7/17 -JM6i -P -ISurface MAT int -Xc -K > ps'. Additional cartographic elements (ticks, legend, scale bar, grid ticks, directional rose, etc) were added by GMT modules. The 1 m bathymetry and xyz grids were interpolated to 30-sec m resolution and grids derived from bathymetry and gravity were visualized. Gravity contours (mGal) were derived with tension factor
of the curvature splines at 0,25. Final data were exported as a 30-sec resolution grid spacing in NetCDF grid. Full GMT code for surface modelling (codes for gravity surface: codes for topographic surface of the xyz ASCII data is available in A.4.9 and A.4.10.

<figure><img src="image_90.png" /><figcaption>Map 3.3.16: Aleutian Trench: surface modelling from the ASCII xyz data. [370]</figcaption></figure>

<figure><img src="image_91.png" /><figcaption>Map 3.3.17: Aleutian Trench and Aleutian Islands: gravimetry map. [370]</figcaption></figure>
# 3.4 AWK

## 3.4.1 Advantages of AWK

AWK, a domain-specific programming language designed for the text processing, was used in this dissertation for table reshaping. The AWK programming language was created at Bell Labs in the 1970s [6] and its name is derived from the surnames of its authors: Alfred Aho, Peter Weinberger, and Brian Kernighan. The GNU implementation of AWK is called GAWK. The available documentation and information on AWK was used extensively in this research [697], [5], [494]. AWK language is highly suitable for sorting multi-field records in a table due to its specially designed syntax and embedded algorithms. Other possibilities include managing small databases, generate reports based on tables, validating data, producing indexes, performing document-preparation tasks, extracting pieces of data for further processing, sorting data [627].

AWK is a data-driven scripting language and a powerful tool for data extraction integrated with commands of Unix operating system and can be run from the console. It consists of a set of codes-actions that can be taken against to operate with streams of table data [630]. AWK enables to extract or transform text, produce formatted records [628]. It uses string datatype, associative arrays indexed by key strings, and regular expressions [629]. Based on this, AWK was selected as a powerful advanced tool for the processing of tables.

<table><caption>**Table 3.4.1:** GMT-derived table for the Cross profile #1</caption><thead><tr><td>No</td><td>Coordinates, lon</td><td>Coordinates, lat</td><td>Location, m</td><td>Azimuth, grad</td><td>Depth, m</td></tr></thead><tbody><tr><td>1</td><td>149.353224965</td><td>41.5059926426</td><td>-200.065204481</td><td>326.61376193</td><td>-5388.91250029</td></tr><tr><td>2</td><td>149.34000552</td><td>41.5210101811</td><td>-198.064782729</td><td>326.61376193</td><td>-5386.2095364</td></tr><tr><td>3</td><td>149.326779938</td><td>41.5360262051</td><td>-196.064356317</td><td>326.605000128</td><td>-5378.97817425</td></tr><tr><td>4</td><td>149.313548214</td><td>41.5510407131</td><td>-194.063925248</td><td>326.596231663</td><td>-5363.02735362</td></tr><tr><td>5</td><td>149.300310343</td><td>41.5660537036</td><td>-192.063489519</td><td>326.587456531</td><td>-5336.36810313</td></tr><tr><td>6</td><td>149.287066319</td><td>41.5810651752</td><td>-190.063049132</td><td>326.578674725</td><td>-5300.09852962</td></tr><tr><td>7</td><td>149.273816136</td><td>41.5960751262</td><td>-188.062604087</td><td>326.569886241</td><td>-5249.45260443</td></tr><tr><td>8</td><td>149.260559789</td><td>41.6110835552</td><td>-186.062154384</td><td>326.561091074</td><td>-5188.08489127</td></tr><tr><td>9</td><td>149.247297273</td><td>41.6260904606</td><td>-184.061700023</td><td>326.552289219</td><td>-5116.35098797</td></tr><tr><td>10</td><td>149.234028582</td><td>41.6410958411</td><td>-182.061241004</td><td>326.54348067</td><td>-5051.39489267</td></tr><tr><td>11</td><td>149.22075371</td><td>41.6560996949</td><td>-180.060777327</td><td>326.534665422</td><td>-5070.64221477</td></tr><tr><td>12</td><td>149.207472652</td><td>41.6711020207</td><td>-178.060308993</td><td>326.525843469</td><td>-5345.61405245</td></tr><tr><td>13</td><td>149.194185402</td><td>41.6861028169</td><td>-176.059836002</td><td>326.517014808</td><td>-5504.76112339</td></tr></tbody></table>

Sequential use of the AWK programming language has been integrated with GMT for geospatial data analysis with a case study of the KKT. Practical technical aim is to analyze and compare bathymetry in the southern and northern part of the trench using digitized cross-section profiles. The initial mapping and geospatial analysis was performed in GMT. The GMT was used for mapping based on ETOPO1 grid and automatic digitizing of the profiles crossing the trench perpendicularly.

The processed geodata were received in a numerical form as a complex multi-field table for each segment. These tables were generated by the GMT in its native format and could not be directly processed by MATLAB/Octave. Therefore, tables were exported to AWK. The table was restructured, sorted and
reshaped by the AWK script. Tab 3.4.1 illustrates Initial table derived from the GMT for the Profile #1, KKT. The geometric model of each profile in Tab. 3.4.1 is as follows: the trench valley corresponds to the zero value of the segment, that is the origin of the plot. Left and right parts of the segments stretch to the 200 and -200 m, correspondingly as shown on X axis. The direction of each profile is NW-SE, that is from the Sea of Okhotsk and Greater Kuril Chain directed to the Pacific Ocean.

The conceptual idea of the use of AWK consists in a chain of data processes by GMT, reshaped by AWK and further converted to Octave. Technical experiment aimed at the analysis of the geomorphic differences of the KKT located in the north-west Pacific Ocean (Fig. 2.1.1) through the analysis of its cross sectioning profiles. For this purpose, the profiles were automatically digitized by the GMT and then the data were converted to the AWK where the tables were reshaped. Further processing and plotting of tables was done in Octave. The geographic aim is to visualize the difference in the northern and southern parts of the trench separated by the Bussol Strait (Fig. 2.1.1).

The use of AWK was to enabled table reshaping for regional geomorphological variations in the trench in its southern and norther parts [410]. The table generated directly by GMT was not readable by Octave, therefore reshaping in AWK was applied. The experimental batch sequence of the data processing through three different tools is advantageous, because they were used sequentially for each step of the data analysis where they suit best: GMT for mapping, geospatial data analysis and generating series of profiles, AWK for reshaping complex tables, and Octave for plotting graphs of the selected profiles. The explanation of the methodology is as follows.

## 3.4.2 Data pre-processing

BATHYMETRIC mapping, modelling, geospatial analysis and plotting profiles by GMT modules using codes provided in Appendix. Map was visualized (Fig. 2.1.1) using ETOPO1 raster grid. Profiles were generated using GMT code. The enlarged view of the profiles is shown on Map 4.2.1 where numbering starts from S to N for both segments. Stacked profiles with a median as a red line and error bars are shown on Fig. 4.2.1. The profiles contained coordinates and depths in each point across the profile (by 2 km). 62 profiles were plotted for northern and 52 profiles for the southern parts of the trench.

These data were stored as two tables in GMT format. The fragment of the initial table is shown on Tab. 3.4.1. The geological conditions associated with the studied geomorphology were analyzed using available literature on the geology of KKT. The geospatial data visualization methods included the following steps: selecting cartographic projections, choosing fonts and cartographic layout, ticks, scale bar, legend, annotations, color palettes, etc. Because the table produced by the GMT automatically during digitizing process (Tab. 3.4.1) was not compatible with GNU Octave and could not be used directly for plotting, the AWK table reshaping was then performed in the 2nd step.

## 3.4.3 Table reshaping

THE second step consisted in reshaping the tables generated automatically by GMT in previous step. These tables were reshaped using AWK programming provided in Appendix. Only necessary columns and rows were selected from the initial tables using AWK script. As a result, two tables were reshaped into the new ones, as shown on Tab. 3.4.2 (fragment) where each column contained depth values by the profiles which could be processed by Octave. Tab 3.4.2 illustrates AWK-reshaped new table containing location and depths (m) for the first 15 rows.

The AWK part of the research followed the GMT aimed at table processing. The .csv table was run directly as a file following by a set of the temporary files (tempfiles) created to store the selected columns from the initial tables as part of a script pipeline. The AWK script (see Appendix) was used for reshaping
initial table (Tab. 3.4.1) that was received directly from the GMT as a result of the automatically digitized profiles.

<table><caption>**Table 3.4.2:** AWK-based reshaped table (fragment): location and depths (m) for the first 10 profiles</caption><thead><tr><th>location</th><th>profile1</th><th>profile2</th><th>profile3</th><th>profile4</th><th>profile5</th><th>profile6</th><th>profile7</th><th>profile8</th><th>profile9</th><th>profile10</th></tr></thead><tbody><tr><td>-200,07</td><td>-5388,91</td><td>-5066,99</td><td>-5096,23</td><td>-5007</td><td>-5095,57</td><td>-5089,61</td><td>-5111,86</td><td>-5105,17</td><td>-5230,93</td><td>-5387,38</td></tr><tr><td>-198,06</td><td>-5386,21</td><td>-5056,08</td><td>-5087,92</td><td>-4990,77</td><td>-5111,1</td><td>-5113,8</td><td>-5126,26</td><td>-5116,11</td><td>-5242,93</td><td>-5395,11</td></tr><tr><td>-196,06</td><td>-5378,98</td><td>-5056,92</td><td>-5076,23</td><td>-4991,77</td><td>-5127,27</td><td>-5143,67</td><td>-5136,49</td><td>-5127,79</td><td>-5263,13</td><td>-5400,87</td></tr><tr><td>-194,06</td><td>-5363,03</td><td>-5068,51</td><td>-5063,31</td><td>-4998,02</td><td>-5125</td><td>-5174,22</td><td>-5137,67</td><td>-5145,92</td><td>-5286,53</td><td>-5405,93</td></tr><tr><td>-192,06</td><td>-5336,37</td><td>-5084,85</td><td>-5052,82</td><td>-5000,48</td><td>-5102,27</td><td>-5197,73</td><td>-5131,9</td><td>-5174</td><td>-5300,04</td><td>-5412,8</td></tr><tr><td>-190,06</td><td>-5300,1</td><td>-5099,18</td><td>-5047,83</td><td>-5001,55</td><td>-5077,16</td><td>-5207,11</td><td>-5117,67</td><td>-5204,43</td><td>-5304,92</td><td>-5417,57</td></tr><tr><td>-188,06</td><td>-5249,45</td><td>-5109,15</td><td>-5047,52</td><td>-5005,99</td><td>-5057,1</td><td>-5203,18</td><td>-5092,88</td><td>-5222,58</td><td>-5303,37</td><td>-5410,98</td></tr><tr><td>-186,06</td><td>-5188,08</td><td>-5114,52</td><td>-5049,58</td><td>-5008,41</td><td>-5041,25</td><td>-5189,31</td><td>-5067,21</td><td>-5232,12</td><td>-5295,12</td><td>-5394,34</td></tr><tr><td>-184,06</td><td>-5116,35</td><td>-5117,15</td><td>-5039,51</td><td>-5017,05</td><td>-5031,61</td><td>-5169,68</td><td>-5049,59</td><td>-5232,31</td><td>-5281,9</td><td>-5380,4</td></tr><tr><td>-182,06</td><td>-5051,39</td><td>-5123,59</td><td>-5024,97</td><td>-5029,04</td><td>-5029,31</td><td>-5146,61</td><td>-5043,71</td><td>-5220,18</td><td>-5265,19</td><td>-5367</td></tr><tr><td>-180,06</td><td>-5070,64</td><td>-5133,68</td><td>-5020,5</td><td>-5038,68</td><td>-5034,43</td><td>-5125,47</td><td>-5044,13</td><td>-5190,62</td><td>-5245,44</td><td>-5364,2</td></tr><tr><td>-178,06</td><td>-5345,61</td><td>-5147,98</td><td>-5028,46</td><td>-5045,35</td><td>-5040,97</td><td>-5109,78</td><td>-5043,65</td><td>-5147,64</td><td>-5224,38</td><td>-5350,02</td></tr><tr><td>-176,06</td><td>-5504,76</td><td>-5156,31</td><td>-5042,72</td><td>-5048,1</td><td>-5047,01</td><td>-5102</td><td>-5033,99</td><td>-5108,91</td><td>-5201,76</td><td>-5320,17</td></tr><tr><td>-174,06</td><td>-5511,42</td><td>-5150,97</td><td>-5057,82</td><td>-5048,89</td><td>-5051,76</td><td>-5108,89</td><td>-5018,05</td><td>-5082,59</td><td>-5179,36</td><td>-5282,81</td></tr><tr><td>-172,06</td><td>-5475,67</td><td>-5137,74</td><td>-5073,48</td><td>-5057,95</td><td>-5052,26</td><td>-5132,15</td><td>-4991,71</td><td>-5070,96</td><td>-5155,37</td><td>-5246,77</td></tr></tbody></table>

This table was then reshaped by AWK into the second table (Tab. 3.4.2) that now had structured view with columns containing only necessary data: depths (m) and profile length (from -200 to 200 m, according to the stretch). The main task of the reshaping was to select necessary data from the initial huge table that contained 12726 rows 5 columns. Therefore, during this procedure the unnecessary rows were removed, and only necessary columns were picked. The AWK research steps were done consequentially as explained below stepwise.

Initially, the working directory was set up to store all the files (both temporary and main ones) in this folder (cd /Users/pauline/Documents/Octave) which was then controlled by pwd check. Then the input file (table2.csv) was read in into AWK environment in units (records) processed by the embedded syntax rules of the AWK one record at a time line by line [143]. Each record was automatically split into pieces, that is fields of the table making it more convenient for AWK to work on the parts of a table [631].

1. At the first step the series of the temporary files was created where all the preliminary results were stored using Unix: mktemp utility. Several kinds of tasks occurred repeatedly when working with table file received from the GMT: extracting one by one columns from the initial file and then storing them as separate files. Therefore, certain lines were extracted from the table and the rest was discarded at this step.

2. The initial GMT tables for two segments of the profiles contained 62 and 52 profiles for northern and southern parts, respectively. Therefore, at this step only two necessary columns were selected using AWK predefined variable FS (field separator) as space. That means, in both input and output data (output field separator as OFS) the items were separated by single spaces.

3. Now only profile 1 was selected and save as a separate table with X and Y (location and depths).
4. Only Y column (field 5) was taken from the initial table. This column contained depths for all profiles, for example for the northern segment these are 12726 rows.

5. The profile I was selected consequentially from the Y column. This is done using command printing section of file between two regular expressions (regex). Patterns in awk control the execution of rules when its pattern matches the current input record. This table now contained information for the depths in all points along the profile nr. 7. Then the same procedure was repeated for all other profiles. The regular expression that matches the text of the table (`'/-L0-07/`, `'-L0-08/'`) is put between the slash lines, like this: `/regular expression/.` Consequently, this pattern was repeated for the northern segment for the profiles of i number: 1, 7, 14, 21, 28, 35, 42, 49, 56, 62. The output data are stored as tables.

6. At this step the single Y column of the profile i was added to the table. Then, the output was saved to table. Making these modifications in the table was done using patterns (see AWK script), but leaving the rest of the table file alone. The same procedure was repeated for the southern segment using the same regular expression as a template.

7. Now the header and tail were deleted from the output table using the AWK expression that prints only lines which do not match regex.

8. The auxiliary files were removed by the rm.

9. The number of lines was counted for initial table. In this case initial table had 12726 lines.

10. Finally, the number of lines was counted for the output table. In this case the final table had only 201 lines (against previous 12726).

As can be noticed, the presented AWK script consists of expressions that are the basic building blocks of AWK patterns and actions. An expression performed in each step of the script evaluates to a value that was printed and consequentially passed to the next expression at the following step in the script. For the case of the temporary files (Step-5), a series of expressions assigned a new value to the table variables using an assignment operator. AWK expressions extensively used in the script served as a pattern for the selected geomorphological profiles along the whole length of the trench, thus facilitating the repeatability of the research steps.

Final output table was reshaped in view, that is columns became rows. Azimuth and absolute coordinates were deleted from the table. As a result, the output table contained 201 rows and 11 columns and initial long table has been reshaped. The new structure of the table is as follows: column 1 contains cross-section line 400 km long (-200:200 km). Columns 2-11 contains depths for the cross-section profiles. Unix programs and utilities used additionally included: `mktemp` for creating a series of temporary files and `rm` for removing preliminary temporary files. The final table was then processed in GNU Octave language for the visualization of the cross-section profiles.
**3.5 Octave/MATLAB**

**3.5.1 Advantages of Octave**

GNU Octave is an open source high-level interactive language, designed for numerical computations by J. B. Rawlings and J. G. Ekerdt. It is compatible with Matlab and is suitable for a variety of scientific tasks: linear and nonlinear equations, numerical linear algebra, statistical analysis, numerical experiments [150]. A significant well-known advantage of the Octave consists in its compatibility with Matlab: although there are some minor differences in their syntax, GNU's Not Unix (GNU) Octave is known as a free open-source analogue of the commercial licensed Matlab which is widely used in both academia and industry by thousands researchers and engineers worldwide.

<figure><img src="image_92.png" /><figcaption>**Figure 3.5.1:** Kuril-Kamchatka Trench: Octave graph based on AWK-reshaped table, north [410]</figcaption></figure>

Therefore, Octave scripts may be converted to the Matlab environment using the syntax and power of Matlab. Technically, as Matlab, GNU Octave has various built-in functions, operates with arrays, vectors, strings, booleans, scalars and matrices, enables plotting of the high-quality graphics, performs mathematical computations, operates with variables and different data types, as demonstrated further in this paper (Appendix: GNU Octave codes). The Octave was used in this work for plotting several geomorphological profiles.

The experimental novelty of using Octave is integration of this typically statistical tool with geospatial mapping in GMT via AWK. Combined together for the geospatial data analysis, these tools enabled to model submarine geomorphology of the KKT. Therefore, the main goal is to develop technical methodology for the batch processing of the geospatial data by the smooth stepwise GMT-AWK-Octave workflow chain. This technical workflow also supported geomorphic studies helping to visualize and better understand structure of the trench. Usually, similar research is maintained in the GIS, but applying advanced data processing scripting tools facilitates research through the automatization. Automatically digitized transect profiles are precise, finely plotted and statistically correct comparing to the hand-made vector drawing using GIS (e.g. ArcGIS, MapInfo, or QGIS).
## 3.5.2 Data processing

**B**ECAUSE the total amount of the KKT profiles overstepped 100 (62 and 52 for the northern and southern trench segments), only selected profiles were visualized. For this purpose, the modified tables were converted to GNU Octave language for visualizing and plotting selected profiles, as further described in this section. Finally, the geomorphology was analyzed and two segments compared.

<figure><img src="image_93.png" /><figcaption>**Figure 3.5.2:** Kuril-Kamchatka Trench: Octave graph based on AWK-reshaped table, south [410]</figcaption></figure>

Plotting graphs of the selected profiles by the GNU Octave was based on the code provided in Appendix. After the table was converted to the necessary format through AWK manipulations, it was then processed by the GNU Octave scripting and two graphs showing northern and southern segments of the trench were plotted (Fig. 3.5.1 and Fig. 3.5.2) by Octave scripts. The GNU Octave scripting was used for visualizing and plotting graph and geomorphological analysis for the selected profiles.

Trench segments were visualized on graphs (Fig. 3.5.1 and Fig. 3.5.2) plotted based on the tables reshaped by AWK in the previous step. The profiles were selected taking into account qualitative and quantitative assessment of data. Quantitative comparability consists in the repeatability of the sampling intervals: every 7th profile in the N segment and every 5th profile in the S segment until profile 40 and profile 52. Qualitative selection was based on the analysis of the bathymetric shapes: similar profiles were removed to avoid repeatability and only varying profiles were compared.

Visualization of the selected graphs of the geomorphology was done by plotting using available technical documentation [518]. Fig. 3.5.1 and Fig. 3.5.2 show selected graphs of the profiles made on basis of the table reshaped by AWK and visualized in GNU Octave using script of m-file, example of N segment. S segment was repeated using the same script template in corresponding table and necessary minor amendments (numbers of profiles, depths for the S segment).

## 3.5.3 Data analysis

**U**SING Octave, a comparative analysis of the geomorphology in selected profiles crossing the KKT was performed visualizing the V shape of the KKT which can be notice in a series of the cross-section profiles (Fig. 4 and 5). The Octave graphs (Fig. 4 and Fig. 5) highlighted
geomorphic patterns in northern and southern parts of the trench visualizing variations in depths and geomorphic shapes. Southern part of the trench proved to have deeper bathymetric ranges and steeper trench gradient slope which is caused by the intensive seismicity and volcanic activities in this particular area of the trench.

The series of the SW-NE-aligned orthogonal bathymetric profiles (Fig. 2) show gradual changes in the geomorphology of the **KKT** along the islands from the Greater Kuril Chain northwards from the Bussol Strait: Simushir, Ketoi, Ushishir, Rasshua, Matua, until it finally approaches the north-eastern end of the Sea of Okhotsk and Bering Sea. GMT/AWK/Octave sequential data processing, automatic digitizing of the orthogonal profiles and graphical plotting based on various algorithms of the data analysis have been tested as an experimental design of this research, aimed to visualize two morphological patterns along the southern and northern segment flanks of the **KKT**.

Trench migration rates are also mainly controlled by the lower continental tectonic plate velocity [350], which in turn depends on the tectonic slab age buoyancy [92]. More detailed study on the tectonic impact goes beyond the scope and size of this manuscript which aims at development of the technical methodology of the batch geodata processing by GMT, AWK and Octave as a chain methodological workflow. As a suggestion for further studies on this phenomena, a comparative analysis of the subduction speed of the tectonic plates in northern and southern parts of the trench could be performed in a separate research based on the availability of the precise geophysical and seismic data.

Combination of various data analysis tools, such as GMT scripting toolset, AWK and Octave programming languages for sorting and reshaping data, plotting graphs and mapping enabled to create geospatial analysis of the selected specific areas of the Pacific Ocean: the **KKT**. Successful application of the programming languages and advanced statistical tools towards geospatial data analysis has been demonstrated in various papers: Python [770], [424], Fortran language [4], [94], R language [633], [434], [429], [422]. However, only few examples were given so far specifically for AWK language and its application is traditionally limited by the IT domains [603], [726]. Current work contributes to filling the gap in the need for using such effective language as AWK [492] in the geological and Earth sciences and absence of the practical examples and full codes that are always useful besides technical manuals.

Understanding complex interrelation and spatial variation of the geological and geophysical factors affecting the **KKT** is crucial to analysis of its geomorphic patterns. Various factors contribute to the trench geomorphology and shape: motion of the Pacific plate subduction to the adjacent small Okhotsk plate, seismicity and location of the earthquakes. Special geophysical settings, geological parameters and special sedimentation processes induced by active volcanism and seismicity in the region are among the crucial factors affecting trench geomorphology. Geological settings of the Kamchatka region in general are notable for the mountains chains in the adjacent land, recorded active seismicity and magmatism, steep gradients of the mountains Hokkaido and Sakhalin islands, and Kamchatka Peninsula, as well as frequent earthquakes and erodible lithology. Necessarily, these factors contribute towards the active sedimentation processes which in turn affect geomorphology of the trench. Large amounts of the sediments are being transported to the Sea of Okhotsk and then driven in SE direction from the Hokkaido and NW from the Kamchatka Peninsula to the neighboring **KKT**.
## 3.6 Quantum GIS

QGIS part of this project was performed in version 3.6.1 Noosa. Submarine geomorphology was studied through the spatial and statistical analysis of the cross-section profiles digitized across the trenches. QGIS data processing included various aspects of the data processing: data search from the web resources, data storage in QGIS platform, digitizing and geospatial processing using GIS plugins, geospatial analysis, data capture and export to the csv tables, transfer data sets to R and Python libraries, statistical processing, modelling and visualizing graphical plots.

<figure><img src="image_94.png" /><figcaption>Map 3.6.1: Mariana Trench: map of the profiles and fragment of the coordinates</figcaption></figure>

The QGIS plugins splits several complex tasks into smaller and more manageable components: reading coordinates, crossing profile lines. The GIS digitizing workflow are illustrated on the Fig. 3.6.2. The fragment of the map and coordinates are shown on Fig. 3.6.1: map of the digitized profiles (left) and coordinates for the selected profile (Nr. 21), Table (right): (1) fragment in geographic coordinate system European Petroleum Survey Group (EPSG) Geodetic Parameter Dataset: EPSG 4326 WGS84. (2), right: fragment in Cartesian coordinate system, re-projected The QGIS methodology includes following steps:

* Geospatial data were collected as thematic layers and uploaded to GIS project including their attributes and coordinate system parameters: bathymetry, geology, topography, tectonics, sediment thickness, slope angle, aspect class, igneous volcanic zones, ophiolites, etc.

* A QGIS project has been projected into the Universal Transverse Mercator coordinate system (UTM) cartesian coordinate system (square N-55).

* 25 cross-section bathymetric profiles were drawn manually across the trench by 'Add Line Feature / Digitizing Toolbar'. Each profile has a length of 1000 km and a distance between each 2 is 100 km. The profiles geometry was digitized along the trenches, attributes (profiles names, coordinates) entered (Fig. 3.6.1)

* The line has been copied and re-entered for the next 24 profiles.

* Numbering of the profiles goes consequently from the north to the south, from 1 to 25.

* Intersections between the lines and bathymetry was done by 'Vector / Analysis Tools / Lines Intersection'. This resulted in the array of the points along each profile. Each point has depth value and a reference number.

* Each profile contained 518 observation sample points, XY coordinates, elevation depth data and thematic geological data (e.g. thickness of the sediment layer) stored in an attribute table. The total data intakes consisted in pool of 12.590 points.
* 'Lat Lon Tools' plugin adding projected coordinated was installed to activate QGIS option of adding XY values to the Spreadsheet.

* Lat/lon coordinates were added to the attribute table of each intersection line by Lat Lon Tools plugin. This plugin was selected as external map support, point digitizing tool to capture and zoom to profiles coordinates. Attribute tables contained data in various segments.

* The coordinates were added to attribute tables of each corresponding profile using plugin function Lat Lon Tools / Conversion / Point layer to fields. The initial 25 layer's latitude, longitude (Y, X) were added and copied into the text fields in the new layers.

* The topography of the 25 cross-section profiles (Fig. 4.12.1) was plotted using ProfileFromPoints plugin of QGIS. The topography data were plotted along a centerline of the cross-section profiles by elevation Z (depth). The resulting profiles are shown on Fig. 3.6.2.

* The attribute information has been derived from each profile and stored in a table containing coordinates, depths and thematic information.

* The tables were exported from QGIS (Export / Save Features As) into csv format for each of the 25 profiles, respectively. The csv data were read into R and Python for statistical analysis using Pandas and NumPy libraries.

<figure><img src="image_95.png" /><figcaption>Map 3.6.2: Philippine Trench: map and digitized cross-section profiles. QGIS</figcaption></figure>

Cross-sectionig transects (Fig. 4.12.1) were digitized in QGIS and their coordinates and attribute values collected as a csv table. A GIS methodology consists in a stepwise illustration of the creating a series of 25 bathymetric profiles, Mariana Trench by digitizing routing of the profiles. Using queries of the QGIS, the data on geology of the underlying surface was retrieved: local geomorphic conditions, variation in the lithology and sedimentation transport processes. Geology and tectonics as attributes were used to analyze their effects on variations in trench geomorphology. Prior to statistical analyses, the bathymetry of the cross-sections on the joint facetted plot was visually inspected (Fig. 4.12.1).
## 3.7 Statistical Algorithms

THE methodology of the current research includes a range of the mathematical algorithms designed for statistical analysis [609]. The functionality of these methods has been tested using statistical libraries of the scripting languages R, Python, Octave and GMT. The statistical methods includes various types summarized in a Table 3.7.1. Various statistical methods were developed for application at geospatial research [733]; [585]; [632]. Existing powerful technologies provided by R, Python and Matlab/Octave enable to perform precise computations and statistical analysis of the geodata and to create data frames.

<table><caption>**Table 3.7.1:** Types of the statistical analysis used by Python and R libraries</caption><thead><tr><td>No</td><td>Type</td><td>Example</td></tr></thead><tbody><tr><td>1</td><td>Data grouping & partition</td><td>Clustering (k-means method), hierarchical clustering, dendrogram trees, clustermaps with marginal dendrograms</td></tr><tr><td>2</td><td>Correlation analysis</td><td>Covariance scatterplot matrices; Correlation ellipses; Pearson, Spearman, Kendall correlograms, 3D correlation plots</td></tr><tr><td>3</td><td>Regression analysis</td><td>Linear regression; Locally Estimated Scatterplot Smoothing (LOESS)</td></tr><tr><td>4</td><td>Data distribution</td><td>Boxplots; QQ, Violin/Letter-value plots; KDE for analysis of the probability of data distribution</td></tr><tr><td>5</td><td>Exploratory analysis</td><td>Ranking dot plots with areas; PCA; EFA; ANOVA</td></tr><tr><td>6</td><td>Quantitative variables</td><td>Pairwise double Y-axis plots; Ridgeline plots, Dumbbell charts, Diverging bars</td></tr><tr><td>7</td><td>Data sorting</td><td>Circular Bar plots, Stemplots; Radar charts</td></tr><tr><td>8</td><td>Set theory</td><td>Euler-Venn diagrams; Conceptual analysis Mindmaps, Word clouds</td></tr><tr><td>9</td><td>Compositional charts</td><td>Categorywise bar charts; 'Waffle' charts</td></tr><tr><td>10</td><td>Descriptive statistics</td><td>Histograms, Heatmaps; line plots, circular pie charts, area charts, strip plots, hexagonal plots, facet grids, stacked area charts for data range across trench segments; stacked bar plots for data distribution by tectonic plates; stacked bar charts;</td></tr><tr><td>11</td><td>Qualitative variables</td><td>Mosaic plots; Silhouette plots; Hierarchical edge plots</td></tr></tbody></table>

The ML and data mining for the marine geological research is nowadays one of the most important technical approaches. Using programming based ML algorithms ensure the precision, effectiveness and objectiveness of the big data set processing, which is always the case for marine geological research. This enables to process big data dividing them into compatible data sets for statistical analysis. Thus, one can select attribute tables specifically for magmatism of the nearby areas to analyze the possibilities of the earthquakes, or to focus on the bathymetric properties and distribution of various depth values across the selected profiles. Statistical analysis of the geospatial data by histograms, KDE Plots, boxplots, regression analysis and other statistical methods has been presented by Python. Tab. 3.8.1 illustrates a summary of the approaches of the statistical analysis used in current dissertation.
## 3.8 R language for marine geology

R programming language is especially attractive to the geologists for the scientific computing due to its logical and straightforward syntax and multi-functional standard libraries [380]. Using R statistical libraries, the unevenness of various factors affecting trench geomorphic structure has been studied. These include such factors as sediment thickness, slope steepness, angle aspect, depth at the basement and magmatism of the nearby areas. Methods includes following major R libraries: `{ggplot2}` for regression analysis, Kernel density curves, compositional charts; `{ggalt}` for Dumbbell charts for data comparison by tectonic plates, ranking dot plots for correlation analysis; `{vcd}` for mosaic plots, silhouette plots for compositional similarities among the bathymetric profiles, association plots; `{car}` for ANOVA.

<table><caption>Table 3.8.1: Algorithms of the descriptive statistical analysis</caption><thead><tr><th>No</th><th>Type</th><th>Method</th><th>Method</th></tr></thead><tbody><tr><td>1</td><td>Variables Relationship</td><td>Ternary diagrams</td><td>Link</td></tr><tr><td>2</td><td>Advanced methods</td><td>ANOVA Hypothesis testing</td><td>SARIMA</td></tr><tr><td>3</td><td>Dataset dimension</td><td>EFA</td><td>PCA</td></tr><tr><td>4</td><td>Advanced methods</td><td>Association plots</td><td>Level plots, Hexagonal binning</td></tr></tbody></table>

Methodological scheme of the statistical analysis of the bathymetric profiles includes four general steps of the statistical analysis of the submarine geomorphology. The geospatial analysis aims to extract a more regular impact factors of the geologic morphology development, whereas cluster analysis and correlation matrices extracts groups and classes from a total cloud of observation points. Statistical analysis performed by means of R scripting [65] is composed of separate written codes and algorithm for data processing. R codes presented below call for various statistical formulae and algorithms used by the ML for data processing.

### 3.8.1 Analysis of bathymetric data distribution

The data distribution for the bathymetric profiles has been computed using various existing methods of the statistical analysis [702] of R, e.g.: plotting histograms using standard libraries or `{ggpubr}` library for plotting common legend, QQ faceted wrap script for cowplot with labels, 3D scatterplot, hierarchical network by `{igraph}` library, animated approach to run animated sequence for 25 bathymetric profiles showing regression analysis, Mariana Trench, hierarchical network of the environmental factors affecting trench by libraries `{ggraph}` and `{igraph}`. Methodological approach of the assessment of data distribution may change depending on the structure and size of the data frame. However, the most essential statistical analysis in analysis of the marine geological data start by question of data distribution. Computing and plotting Kernel distribution curves, boxplots, regression lines (probability of bathymetric data across profiles), quantile statistics, empirical distribution density function and other methods can be distinguished as the most useful for data processing.

#### 3.8.1.1 Box plots

The boxplot of depths distribution show summary statistics: range and quartiles of the observation points across the cross-section profiles. A boxplot (Fig. 3.8.1) graphically displays frequencies of the bathymetric observations data set of the 25 profiles of MT. It plots frequency, or count, on the depth on y-axis (vertical) and measured variable (number of profile) on the x-axis. The boxplot was generated by the default R library (stat) graphically representing five most important descriptive values for a bathymetric data set: A.6.2. These values include the following bathymetric values: minimum value, the 1st quartile
(0.25), the median, the 3rd quartile (0.75), and the maximum value (absolute depth). When graphing this five-number summary, only horizontal axis displays values, while vertical line is placed above each of the summary numbers showing depth values. A box is drawn around the middle three lines (1st quartile, median, and 3rd quartile) and 2 lines are drawn from the boxs edges to the 2 endpoints (minimum and maximum depths values) using R code: A.6.2

<figure><img src="image_96.png" /><figcaption>**Figure 3.8.1:** Mariana Trench: boxplots on data distribution by profiles. R code: A.6.2</figcaption></figure>

### 3.8.1.2 Histograms

The histogram plots (Fig. 4) [422] show normal distribution of the depth observation data by the cross-section profiles, their outliers, skewness, median, mean, average, maximal, minimal and quartile distribution values. For each of 25 bathymetric profile various colors are taken for the following statistical data of depth values along the profile: black curves stands for normal distribution, blue curves for density distribution. Vertical dashed lines: purple median values, green mean values. The histograms was drawn by R library {ggplot2}. On the 1st step a single histogram for each of 25 profiles was created. The full R script to create one histogram is as follows (here: for the profile Nr. 1, further applied for every one from 25 profiles by changing the name of the file from 01 to 02, and so on to 25). Using this script further 25 profiles were created, consequently, p01, p02, ... p25. The combination of histograms of the 25 profiles on one layout was done using following R script: A.6.1. The histograms illustrate the frequency of the score occurrences in a continuous data set of the bathymetric depths that were divided into classes showing variation in the samples of the depths, assumptions of the underlying statistical distribution of observation points values.

### 3.8.1.3 Violin plots

The violin plots (Fig. A.6.25) show Kernel probability density distribution of the bathymetric observations, as multimodal distributions with multiple peaks. The {violinplot} package was installed and
activated prior to the current work. Kernel density distribution plot as shown on the Fig. 3.8.2 was cre-
ated using library {violinplot} of R in a combined plot, which includes calculated quantiles for 0.25
and 0.75 of the data pool.

<figure><img src="image_97.png" /><figcaption>**Figure 3.8.2:** Mariana Trench: violin plots on depth distribution by profiles. R code: A.6.25</figcaption></figure>

3.8.1.4 Notched boxplots

Another way of representing boxplots is 'notched boxplots'. These were computed and visualized in R using code A.6.3 aimed to show depths distribution with majority of the data from -3000 to -5000 m and outliers (the deepest points) show that from profile 1 to 16 general depth increase with the deepest values at profiles 10 and 11. The notched boxplot on Fig. 3.8.3 shows the shallowest parts of the trench located on profiles 24 and 25 where MT crosses YT, near Yap Island in west Pacific Ocean, [422].

<figure><img src="image_98.png" /><figcaption>Figure 3.8.3: Mariana Trench: notched plots on depth distribution by profiles. R code: A.6.3</figcaption></figure>

3.8.2 KDE

The KDE is a non-parametric way to estimate the probability density function of a random variable of the samples. The KDE was performed by smoothing data points across the sample points in the profiles 1:25 of the trench. The KDE can be endowed with properties such as smoothness or continuity of the
bathymetric data by using a suitable kernel. Technically, besides {violinplot} library, a KDE function can be implemented in R through density, {bkde} function in the {KernSmooth} library, and through {kde} function in the {ks} library: Fig. 3.8.4. To use {kde.R} function, it is not required to install special packages. The KDE shows the majority of the depths of the trench between 3000 and 6000 m. The graph was made as a combined plot (Fig. 3.8.4) enabling to compare overlapping and max aptitude of the density curves for the profiles and depth distribution by tectonic plates.

<figure><img src="image_99.png" /><figcaption>**Figure 3.8.4:** Mariana Trench: KDE for the 25 bathymetric profiles, R code: A.6.19</figcaption></figure>

### 3.8.3 Regression analysis

Regression analysis was done to estimate the relationships between variables using existing methods [301]. The principle of the regression analysis is based on the analysis of the probability of the data values according to their actual distribution.

<figure><img src="image_100.png" /><figcaption>**Figure 3.8.5:** Example of the regression analysis: Mariana Trench, selected profiles</figcaption></figure>
A regression analysis represents outlying depths by each profile, as shown on Fig. 3.8.5, code (3 profiles) or 1 profile and faceted plot of 25 profiles.

<figure><img src="image_101.png" /><figcaption>**Figure 3.8.6:** Philippine Trench: Regression analysis of 25 profiles. Plotting: R code A.6.5</figcaption></figure>

For the assessment of the variation in the geomorphology of the trench, it is of significant advantage to visualize the data set using regression analysis (Fig.3.8.6). Regression analysis determines the range, general trend in data distribution and scale of the changes in the observed bathymetry by analyzing changes in the values of the elevation and (Fig. 3.8.6). The following R codes were used for the four
methods of the regression analysis visualized on Fig. 3.8.6: A.6.5. Advantages of plotting a facetted graph consists in visualizing a stream of the data at one look. The algorithms and approaches for the regression analysis use existing methods [80]; [247]; [125] and embedded in {ggplot} R library.This methodology utilizes the relation between two or more quantitative depth variables so that one variable can be predicted from another. Thus, one can estimate the probability that the depth values will be located in this or that interval of values profile by profile according to the machine-based algorithm of the probability of their distribution.

Regression analysis and facetted plot wrapping for 25 bathymetric profiles were executed by R library {ggplot2}. Using non-parametric regression analysis the conditional expectation of the dependent variable (that is depth values) was estimated given the independent variables. The confidence interval for each of 25 profiles was done by 4 methods, visualized by the facetted plot wrapping for the series data comparison. The regression analysis reveals how typical value of the criterion variable changes are, when any of the independent variables change, while other independent variables are fixed. The Formula 3.8.3 describing regression equation describing the relationship between two variables used in this research is as follows:

$$Y = a + bX + e(1) \tag{3.7}$$

The literature refers [301] to three major applications of the regression analysis including general loess method ('glm'), loess method ('lm') and LOESS, all of which were tested in this research. A smooth curve through a set of data points obtained with LOESS statistical technique is called a LOESS curve, Fig. 3.8.6 and 3.8.5. LOESS data smoothing is based on the mathematical assumption that a function can be approximated in a neighborhood by a low-order polynomial and simple models can be fit to data easily. LOESS uses the tri-cube weight function computed by Formula 3.8.3:

$$w(x) = (a - |d|^{3})^{3} \tag{3.8}$$

The subsets of data used for each Weighted Least Squares (WLS) fit in LOESS are determined by a nearest neighbors algorithm using Savitzky-Golay filter, applied to a dataset for smoothing the data 3.8.3

$$Y_j = \sum_{i=\frac{m-1}{2}}^{\frac{m-1}{2}} C_i y_{j+i}, \quad \frac{m-1}{2} \le j \le n - \frac{m-1}{2} \tag{3.9}$$

Where dataset includes n $x_j, y_j$ points ($j = 1, ..., n$), where x is an independent variable and $y_j$ is an observed value, treated with a set of m convolution coefficients, $C_i$.

Hence, the quantiles were drawn by calling (geom-quantile) function of R 3.8.5. The Generalized Linear Models (GLM) method was used to give a symbolic description of the linear predictor of depth variability within each bathymetric profile, and a description of the error distribution in a range of the depth observation points. Example of statistical quantile-quantile (QQ) plot for the 25 cross-sectioning profiles is given on Fig. 3 in [426]. Quantiles graphically depict groups of the depth values through their distribution mode, according to the methods of the descriptive statistics.

The spacings between the different parts of the box indicate the degree of the dispersion, i.e. spread, skewness in the bathymetric data, and show outliers. The data outside the upper and lower quartiles are the outliers that are plotted as individual points. The outliers are the non-usual values of the depths indicating specific features of the geomorphology with extra deep values. The histograms show probability of the density distribution of the bathymetric observation by the profiles, i.e. how likely is that the data will belong to this or that part of the group values within certain depth range.
**3.8.4 Comparative geospatial analysis by tectonic plates**

Mariana Trench is located in the collision of four tectonic plates: Mariana, Philippine Sea, Pacific and Caroline. Therefore, a comparative analysis of multidimensional data was performed to analyze inter-relationships between the geographic location of the specific parts of the trench and factors that may impact changes in the geomorphic systems. Using existing methodology [73], various dependences between the impact factors were tested: Fig. 7 in [422].

First, the distribution of various environmental factors affecting the morphology of the trench were analyzed: slope angle, sediment thickness, depth and closeness of the igneous volcanic zones to the selected bathymetric profiles.

Second, changes in the trench deepest angle, calculated as tg°(A/H), where A is the width, and H is the maximal depth of the profile, were calculated for each of the 25 bathymetric profiles and compared, respectively.

<figure><img src="image_102.png" /><figcaption>**Figure 3.8.7:** Four tectonic plates, trench slope angle in the deepest point</figcaption></figure>

Third, the depth distribution by four tectonics plates was visualized using multiple scatter function of R. A pairwise correlation of the environmental factors was performed by the scripting approach of R, which is convenient not only for looping over parts of the geomorphic model of the Mariana Trench, but also for pairwise iterating over bathymetric and geomorphic variables for the comparative analysis. Called functions 'geom-smooth', 'geom-point' and 'geom-quantile' were used to draw smoothed conditional means and to perform regression analysis. The multiple panel strip plots were generated in a combined plot using `{LatticeExtra2}` package of R by strip plot function.

The distribution of the slope angle, sediment thickness, depth and closeness of the igneous volcanic zones towards Mariana Trench, by four tectonic plates is shown on Fig. 7 in [422]. It shows distribution of slope angle, sediment thickness, depth and igneous volcanic zones, by four tectonic plates. The R code used: `A.6.17` (for several strip plots on one layout) (single strip plot with multiple parameters). The multiple panel by groups analyzing four tectonic plates, trench slope angle in the deepest point, depth
scope (min-max) by profiles using overlapping points, are computed and visualized in R code: A.6.18, scatterplot with overlapping points and brackets using library {pBrackets} A.6.20.

<figure><img src="image_103.png" /><figcaption>**Figure 3.8.8:** Multiple panel by groups: tectonic plates, depth by profiles</figcaption></figure>

Variation of the slope angle in the deepest point and depths by 25 bathymetric profiles of the Mariana Trench according to four tectonic plates that is crosses is illustrated on Fig. 3.8.7 and Fig. 3.8.8 (overlapping points approach), [422]. The graphs show variation of the slope angle in the deepest point and depths by bathymetric profiles of the Mariana Trench for four tectonic plates: Mariana, Caroline, Pacific and Philippine, respectively.

### 3.8.5 Ridgeline, categorywise and circular plots

The ridge plots (R: A.6.35 combine information about statistical data on 25 bathymetric profiles: arithmetic density distribution, mean and standard deviation along the profiles in a facetted way, (Fig. 3.8.9. They have lines extending horizontally from the boxes indicating variability of the depths distribution across all the profiles. The code was written in `{ggrids}` and `{ggplot2}` R libraries A.6.35. To create categorywise plot, first, the table was read into R: Second, the tables columns were merged calling `{merge}` function from `library(data.table)` by environmental variables: depth, tectonics, slope angles. Third, the new data frame (DFDT) was created. Fourth, a plot function was called by the code using data frame created in previous step: A.6.13 and a graph was visualized (Fig. 12 in [422], left.

To better understand the behavior of the spatial distribution in elevations in a stream view, a ridge line technics of the landforms data visualization was applied by R [609]. From the graph showing ridge lines (Fig. 3.8.9) one can see comparative visualization of the density distribution of depths by profiles. Ridgeline profiling provides insight into the statistical distribution of the depth and elevation ranges of the trench and archipelago. From these profiles (Fig. 3.8.9), it is possible to comparatively analyze and observe distribution of an elevation values varied for 25 cross-section profiles. The elevation curves illustrate ranges across the profiles transversal to the main channel of the trench (Fig. 3.8.9).

A `{tidyverse}`, another R package, was used to visualize a circular plot showing the comparative analysis of the bathymetric ranged by the sediment thickness. The algorithms and syntax of R coding was applied from the available manuals [609]; [329]. The circular bar plot shows distribution of depths
by tectonic plates (Fig. 12 in [422] right, and the distribution of the sediment thickness (%) according to
the slope steepness by profiles. To create circular barplot divided by groups, following code was used:
R: A.6.16.

<figure><img src="image_104.png" /><figcaption>Figure 3.8.9: Mariana Trench: Ridgeline plots</figcaption></figure>

**3.8.6 Dumbbell charts for data pairwise comparison by tectonic plates**

Since Mariana Trench is located in the collision of four tectonic plates, the comparative testing of the
values distribution by the tectonic plates has been completed through the Dumbbell chart plotting.

<figure><img src="image_105.png" /><figcaption>Figure 3.8.10: Mariana Trench: ranking dot plots by data grouping, by R code: A.6.12</figcaption></figure>

Dumbbell chart is a visualization aimed to give an insight of how the margin tectonic plates constitute
to the Mariana Trench pairwise. It is one of the new statistical methods that initially was widely used in
biosciences, yet can be applied to spatial analyses when it is necessary to take pairwise comparison of the
data distribution from any thematic layer. Plotting was carried out by calling R libraries {ggplot2} and
`{ggalt}` by code: A.6.12. Dumbbell chart demonstrates pairwise distribution of the bathymetric points constituting the continental plates, to show the composition of the sedimentary coverage by four tectonic plates, Fig. 3.8.10 [429].

### 3.8.7 Ranking dot plots by data grouping

Pairwise correlation was visualized by ranking dot plots, an effective tool to perform data grouping by variables (tectonic plates). Plotting distribution of the sample points by igneous volcanic areas aimed at visualizing areas affected by magmatism. Large igneous volcanic areas contribute towards localization of possible earthquake zones. As can be seen (Fig. 3.8.11), the volcanic areas increase along the Mariana Trench SW. Alike to the volcanic areas, variation of the steepness slope was performed by ranking dot plot: profiles with the steepest angles across the trench are located in the NE with a slight decrease towards the SW. The tg° of the profiles were calculated by a standardized formula: a relation of the max depth by the profiles divided by the width of each corresponding profile.

<figure><img src="image_106.png" /><figcaption>2. 马里亚纳海沟. 剖面1-25. Mariana Trench, Profiles Nr.1-25.
   统计图表. 地貌聚类分析. 圆点图. Ranking Dot Plot.
   Variation of Steepness Angles
   tg°(A/H) by profiles 1:25</figcaption></figure>

<figure><img src="image_107.png" /><figcaption>**Figure 3.8.11:** Mariana Trench: ranking dot plots by data grouping, by R code: A.6.14</figcaption></figure>

Crucial points (Fig. 3.8.11) were selected using library `{ggalt}` by calling following R code: A.6.14. The distribution of the steepness values across the trench were shaped by the components of the igneous volcanic areas near the profiles, as seen from the Fig. 3.8.11: profiles # 20, 22, 24, 8 and 11 with notable amount of igneous volcanic observation points (over 180) correlate with steepness angle.

### 3.8.8 Ternary diagrams, radar and pairwise double-Y axis charts

Ternary diagrams and radar charts show the triple correlation of several factors across tectonic plates: slope morphology classes, tectonics, igneous volcanic areas, aspect degree, slope angle, sediment thickness [416]. The programming R code for the ternary diagrams is R: A.6.29. Radar chart was plotted by `{fmsb}` library: A.6.21 (plotting radar chart for depth distribution within bathymetric profiles, Mariana Trench (modelling distribution of the profiles across the trench).

The pairwise double-Y axis was enabled to analyze the effect of the *bi*-factor correlation between the slope angle and other variables (Fig. 3.8.12). The following script runs the model in each case, calculates the X-Y dependences and places both graph on one plot, which is very convenient for comparative studies. The following is a code for the double-Y axis with correlations "Slope angle", "Sediment thickness"
by library {LatticeExtra} in R in consecutive steps: A.6.15. At the end of each run, the final results are
copied to the p1, p2 and p3 plot, and afterwards combined three graphs on one plot, [422].

<figure><img src="image_108.png" /><figcaption>**Figure 3.8.12:** Pairwise double-Y axis charts showing correlation between factors</figcaption></figure>

**3.8.9 Ranking data for variation of the geomorphic steepness**

Geomorphic variation distinguishes a certain segment of the profile of the trench from neighbor areas.

<figure><img src="image_109.png" /><figcaption>**Figure 3.8.13:** Mariana Trench: ranking data for steepness variation, R code: A.6.22</figcaption></figure>

Variation in geomorphology changes in shape that varies from one segment to another depending on
the tectonic properties and geospatial location. Correlation between the geomorphic shape with geologic
settings and geographic location revealed that the submarine geomorphology of the Mariana Trench
is influenced by the geospatial settings of the underlying tectonic plates that affect its shape and cause variation of the steepness across the trench segments. The variation of the shape form is strongly affected by the geospatial location of the trench crescent.

Variations in the presented bathymetric profiles shows (Fig. 3.8.13) trench varying in shape across its crescent: five divided regions varying in steepness degree, ranging from the strong slope, very strong slope, extreme slope and very steep slope. The unsorted spatial series of the profiles (Fig. 3.8.13, left) only shows the consequent distribution of the slopes, so we cannot judge their variations, but only general range. Upon sorting and ranking profiles (Fig. 3.8.13, right), the segments with similar geomorphic slope steepness are ranked and grouped. Ranking data was based on the geomorphic steepness using R libraries and existing mathematical algorithms [609]. The methodology for plotting Fig. 3.8.13 consists in the application of the R code A.6.22.

### 3.8.10 Compositional charts of the determinants variations

The constitution of the system composition has been visualized by the categorical plotting. One of the methods enabling categorical plotting is presented by compositional charts, sometimes referred to as waffle charts.

<figure><img src="image_110.png" /><figcaption>**Figure 3.8.14:** Mariana Trench: compositional chart of the geometry of the profiles, R</figcaption></figure>

The principal approach of this method is to demonstrate the division of the whole system by parts in percentage. Several R libraries were used to perform technically compositional plotting, of which the {ggplot2} was selected as the most effective (A.6.37) comparing to the {waffle} library (A.6.38): it enables to control the appearance of the plot by adjusting details (colors, font size, types, facetted plotting). The compositional chart is aimed to compare distribution of the data across the area, and composition of the aspect class and slope angle for the trench. As can be seen (Fig. 3.8.14), the category very steep slope is the dominant among all other geomorphological types of the slope degree of the Mariana Trench. Likewise, N and SW aspect of the slope direction are best describing the geometry of the northern and southern part of the trench, respectively.

### 3.8.11 Mosaic, silhouette and association plots

Mosaic chart aims at categorical comparison of the geomorphological features of the trench by four tectonic plates (Fig. 3.8.15). Visualizing mosaic plot is a statistical method that involves a subdivision of a rectangular tiles into the areas that represent the conditional relative frequency for a cell in the contingency table. The algorithms and approaches used in mosaic plotting vary slightly. Upon examination
of possible packages, the R library {vcd} was selected for this research (R code: A.6.31. Using {vcd} library, each tile is colored to show the deviation from the expected frequency, that is residuals, from a Pearson X2 or likelihood ratio G-test showing maximum likelihood statistical significance by Formula 3.8.11:

$$G = 2 \sum_{i} O_{i} * \ln\left(\frac{O_{i}}{E_{i}}\right) \qquad (3.10)$$

The algorithm set includes execution of the R code. The convenience of the mosaic plot for the geomorphological analysis of the ocean trenches consists in its visual representation of the association between the environmental variables. Thus, it gives an overview of the data structure and enables to recognize relationships between the different environmental variables. Conversely, boxes across categories constituting Pacific and Mariana plates have similar areas. The area of the tiles, the bin size, includes the identification of the sampling data, giving the proportional value to the number of observations within that category.

<figure><img src="image_111.png" /><figcaption>Figure 3.8.15: Mosaic and silhouette plots for the Mariana Trench</figcaption></figure>

The interpretation and validation of consistency within clusters of geomorphologic and geological variables have been tested using silhouette plotting (Fig. 3.8.15). This technique provides a succinct graphical representation of the suitability and fitness of each from 518 observation points across 25 bathymetric profiles within the clusters. The silhouette value measured similarity of the points derived from the thematic layers (sediment thickness, slope steepness) to the cluster by cohesion and comparison them with other clusters, and separating them from the distinct points. Lower sediment thickness area lie within one class while a class of areas with high level of sediment level fits to another.

The silhouette was calculated with Euclidean distance metric by calling library {cluster} using following R code: A.6.30. The Euclidean distance metric is computed by Formula 3.8.11:

$$\|a - b\|_2 = \sqrt{\sum_{i} (a_i - b_i)^2} \qquad (3.11)$$

The R script plots the silhouettes, returns the values in the n-by-1 vectors. A deviations from the independence of rows and columns containing bathymetric data and environmental variables is presented in a two-dimensional contingency table of the Cohen-Friendly association plots (Fig. 3.8.16). Extended asso-
<table><caption>**Table 3.8.2:** Statistical algorithm of clustering, R code: A.6.9</caption><thead><tr><th>No</th><th>Step</th></tr></thead><tbody><tr><td>1</td><td>Computed hierarchical clustering and cutting the tree in the k-clusters</td></tr><tr><td>2</td><td>Computed pairwise standard scatter plot of k-means cluster correlation</td></tr><tr><td>3</td><td>Computed pairwise standard scatter plot of k-means cluster correlation</td></tr><tr><td>4</td><td>Computed centers (2 to 7 were tested, in total 6 centers) of each cluster</td></tr><tr><td>5</td><td>Created groups of the correlation matrix tested by distributions on the plates</td></tr></tbody></table>

ciation plot in R reveals relationships between the variables by the assoc() function in the {vcd} package. In this case, it is tectonic plates and geomorphological values (bathymetric depths, slope steepness).

<figure><img src="image_112.png" /><figcaption>**Figure 3.8.16:** Association plots for the Mariana Trench</figcaption></figure>

Pearson's residuals on the association plot report test for the statistical goodness of fit: it enables to assess the model fit at each of samples for each profile for the generalized linear models within the context of the chi-square, as shown further in [429]. It has been computed as a raw residual divided by the square root of the variance function. Technically, association plot method includes R code: A.6.31.

### 3.8.12 Clustering (k-means method), correlation ellipses and data partition

Finding correlations among the geomorphic factors that impact the morphological formation of the profiles was done by cluster analysis and correlation ellipses. The k-means clustering is a machine learning technique for classification performs data partition into set up number of clusters [488], [444]. K-means clustering was done by R library {factoextra} using codes for (pairwise cluster analysis) and cluster analysis.

The goal of the k-means clustering was to perform partition of the dataset of the observations across the trench into clusters, in which each depth observation belongs to the cluster with the nearest mean, serving as a prototype of the cluster. Clustering is sensitive to the initial random selection of the cluster centers. Therefore, several clusters were tested starting from 2 to 7. This function provides a solution using a hybrid approach by combining the hierarchical clustering and the k-means methods. The workflow procedure include following steps: Tab. 3.8.2.
The k-means of the clustering of the trench profiles was computed and visualized in R by code A.6.9, Fig. 14, the k-means clustering with different k-values is shown on Fig. 15, [422]. The k-means cluster analysis has been performed using set of cluster centers as the initial cluster centers. The algorithm worked iteratively to assign each profile to one of the groups (2 to 7) based on the morphometric features of the trench across these profiles that were clustered based on their geomorphic feature similarity. The results of the k-means clustering are shown on the Fig. 3.8.17, further details are given in [423]. The pairwise standard scatter plot of k-means cluster correlation distributed by four tectonic plates is shown on Fig. 3.8.17.

<figure><img src="image_113.png" /><figcaption>Figure 3.8.17: Mariana Trench: classification of the k-means clustering by plates, R code A.6.9</figcaption></figure>

The k-means cluster analysis was done using following R code: A.6.9. The mathematical background is as follows. k-means clustering aims to partition the n observations into k (≤ n) sets S = {S1, S2, ..., Sk} so as to minimize the within-cluster sum of squares (WCSS) by the Formula 3.8.12:

$$
\operatorname{arg\,min} \sum_{i=1}^{k} \sum_{x \in S_i} \|x - \mu_i\|^2 = \operatorname{arg\,min} \sum_{i=1}^{k} |S_i| \operatorname{Var} S_i \quad (3.12)
$$

Where $\mu_i$ is the mean of points in $S_i$. The advantages of the k-means clustering applied for Mariana Trench consists in the algorithm nature: rather than defining groups before looking at the data, clustering enabled to find and analyze the groups in the profiles that have formed organically. As the structure of the tectonics and trench geomorphic properties are rather complex, clustering facilitated data partition and grouping. The correlogram ellipses were generated to visualize correlation between the factors affecting the morphology of the trench [422]. The methodological approach includes calling {ellipse} library and executing following script: A.6.32.

3.8.13 Hexagonal plots and data partition

Plotting hexagonal plot was performed using R library {hexbin}. This plot shows how bathymetric
data are distributed by the profiles in terms of frequency: how often the values are being repeated by
profiles. Hexbin plots serve as a useful alternative to scatter plots, as the data are dense. Therefore, plotting each point individually in this case is not effective: the data would be indistinguishable on the graph. Hence, the {hexbin} library was used to plot the range of the depths values within Mariana Trench using the R code: A.6.27. Studies revealed that major depth observation points of the Mariana Trench are located in between the -3000 and -5000 meters (Fig. 3.8.18), as also discussed by [212]. The significant point is that sediment thickness changes notably both within the trench by profiles (1:25) and between four adjacent tectonic plates: Philippine Sea, Pacific, Mariana and Caroline. Data grouping and partition was performed by dividing angles across the trench on to classes (Fig. 19 left [422]). The level plot heatmap (Fig. 19 right [422]) shows classes of the sediment thickness by Mariana Trench and modelling correlation between sediment thickness, slope angle and bathymetry, was built by R code: A.6.33.

<figure><img src="image_114.png" /><figcaption>Figure 3.8.18: Mariana Trench: hexagonal data distribution, by R {hexbin} library.</figcaption></figure>

### 3.8.14 Correlation analysis of impact factors

Correlation coefficient is a numerical measure of direction and strength of the linear correlation between variables, i.e. how strong a relationship is between geomorphological, bathymetric and geological factors is across four different tectonic plates. Three methods of computing correlation analysis were tested: Pearson correlation, Spearman correlation, Kendall correlation [434]. Additionally, numerical correlogram, correlation and cross-correlations matrices were generated to analyze environmental impact factors for the trench. The correlogram matrices are computed and visualized by R scripting libraries.

#### 3.8.14.1 Pearson correlation

Developed by Karl Pearson from a related idea of Francis Galton in the 1880s, a Pearson product moment correlation coefficient, a. k. a. the bivariate correlation, is a measure of the linear correlation between two variables by Formula 3.8.14.1:

$$r = \frac{\sum (x - mx)(y - my)}{\sqrt{\sum (x - mx)^2 \sum (y - my)^2}} \quad (3.13)$$

In this case, the environmental variables and geomorphological factors were tested by Pearson correlation: sediment thickness, slope angle and class, aspect degree, hill shade in values, location of igneous
volcanic zones across the abyssal valley, tectonics by 4 plates, bathymetric values, that is max/min depths, median, means for each profile, tan angle°. Pearson correlation has a value between +1 and −1, where 1 is a total positive linear correlation, 0 is no linear correlation, and −1 is total negative linear correlation. It is presented on Fig.3.8.19, computed and visualized by R script: A.6.41.

<figure><img src="image_115.png" /><figcaption>Figure 3.8.19: Mariana Trench: Pearson and Spearman correlation matrices, by R code A.6.41</figcaption></figure>

### 3.8.14.2 Spearman rank correlation coefficient ρ

The environmental variables were tested by Spearman rank correlation coefficient ρ: sediment thickness, slope angle and class, aspect degree, hill shade in values, location of igneous volcanic zones across the abyssal valley, tectonics by four plates, bathymetric values, max/min depths, median, means for each profile, tan angle°. Named after Charles Spearman, a Spearman non-parametric ρ coefficient is handling ties between the variables less properly comparing to the previously computed Pearson coefficient (Fig.3.8.19). The Spearman's ρ was computed on sets of paired rankings between the factors. On Fig.3.8.19 the more correlated variables are represented by the colors of green. The size of the circle also depends on the correlation value. Spearman's rank correlation coefficient is denoted by ρ or r_s. It is given by the following formula [563]:

$$r_s = 1 - \frac{6 \sum (d_i^2)}{n(n^2 - 1)} \qquad (3.14)$$

where

$d_i = rg(X_i) - rg(X_y)$ is the difference between the 2 ranks of each observation, n is the number of observations.

Here $d_i$ represents the difference in the ranks given to the values of the variable for each data item received and processed from the csv table via library {data.table}. The formula of Spearman's rank correlation coefficient ρ is applied for case with no tied ranks. Technically, R script for Pearson correlation was changed for case of Spearman's ρ by following code: A.6.41. On Fig.3.8.19 the measure of rank correlation, that is a statistical dependence between the rankings of two variables has been computed using {GGalt} library by R code: A.6.41. Spearman correlation shows (Fig.3.8.19) the relationship between every pair of two variables described as a monotonic function and comparing following factors: tectonics, geology, geographic, bathymetric and magmatism. Spearman's correlation assesses monotonic relationships, whether linear or not, between the available factors. If there are no repeated data values, a
perfect Spearman correlation of +1 or 1 occurs when each of variables is a perfect monotone function of the other.

### 3.8.14.3 Kendall τ coefficient

Similar to Spearman, Kendall τ coefficient to measure rank correlation, invented by M. Kendall, is a non-parametric test for statistical dependence between environmental ordinal (or rank-transformed) variables. Kendall rank correlation distance is defined as follows:

$$ \tau = \frac{2}{n(n-1)} \sum sgn(x_i - x_j)sgn(y_i - y_j) \qquad (3.15) $$

Kendall rank correlation τ coefficient is statistically used to measure ordinal association between two measured quantities: geological factors in this case. Based on non-parametric hypothesis τ coefficient test for statistical dependence, it was computed by R code A.6.41. It measures rank correlation: the similarity of the orderings of the data when ranked by each of the quantities. Unlike Spearman's, Kendall can handle ties. There are three Kendall τ statistics: τ-a, τ-b, and τ-c, of which τ-b is specifically adapted to handle ties.

By τ-b statistics, each pair of members of the environmental variables (tectonics, geological factors, geographic location, magmatism, sediment thickness) was tested against igneous volcanic spots. It has the same ordinal value by a divisor term. The latter represents geometric mean between the number of environmental factor pairs that are not tied on factor 1 and the number not tied on 2. Visualized Kendall τ-correlation (Fig. 7) for geomorphological factors of the trench and numerical scatterplot matrix are computed by R code: A.6.41, [433], libraries {Matrix}, {mclust}, {correlplot}. Numerical scatterplot matrix of correlation between the factors and cross-correlation matrix of the depths in profiles of the trench (Fig. 9, [433]) are plotted by R code: A.6.7 and visualized by library {ggimage}.

### 3.8.15 Calculation of the normalized steepness angle of the Mariana Trench

The normalized steepness was computed as a steepness of dominance hierarchies of the slope angle (Fig.3.8.20).
<figure><img src="image_116.png" /><figcaption>**Figure 3.8.20:** Mariana Trench: normalized steepness of the slope angles by profiles, R</figcaption></figure>

Steepness is defined as the absolute slope of the straight line fitted to the normalized scores. The normalized scores were obtained on the basis of dominance indices corrected for chance or by means of the proportions of wins. Given an observed data set, it computed hierarchy's steepness and estimated statistical significance by means of a randomization test for the steepness of the trench (here: example of the normalized steepness angles by the bathymetric profiles, Mariana Trench, single-colored (here: plotting 'lollipop charts' and diverging bars of the Mariana Trench bathymetry).

### 3.8.16 PCA

The PCA (Fig. 10, [433]) performed using R code A.6.8 enabled to visualize eigenvectors as showing major direction and vector length for the principal components affecting categorical values. Through the PCA a statistical procedure using an orthogonal transformation to convert a set of bathymetric depth observations of possibly correlated variables, was performed. Thus, the direction of the eigenvectors shows depths, highly influenced by its geomorphic settings: location by 25 profiles and similarities between the profiles.

### 3.8.17 Hierarchical cluster analysis

In this dissertation, the main method of the determining clusters of the environmental variables affecting trench formation is the algorithm provided by cluster analysis for recognition of the data patterns, executed by {dendextend} R library, code A.6.6, [402]. R enabled to perform stepwise cluster analysis, data partition, sorting and grouping. Comparative analysis of the multi-dimensional data enabled to analyze inter-relationships between the factors impacting changes in the morphology and dynamics of the trench system. The most efficient method was computing hierarchical dendrogram clustering with p-values called by R libraries {dendextend}, {gridExtra}, {magrittr}, {methods} and {gridGraphics} [402]. Used hierarchical clustering method is default {hclust}.
The results of the algorithm are visualized as a plotted dendrogram. Working flowchart include following steps: creating unsorted and sorted dendrograms (Fig. 3.8.21), cross-comparison of the both, hierarchical clustering using p-values, and dendrogram by clustering using bootstrap probability method (Fig. 3.8.22). Bootstrapping is based on using a scheme that creates new data sets through reweighting and resampling the initial dataset. Given a set of N data points, the weighting assigned to data point i in a new data set D^J is $w_i^J = x_i^J x_{i-1}^J$, where x^J is a low-to-high ordered list of N - 1 uniformly distributed random numbers on [0, 1] preceded by 0 and succeeded by 1. Using this approach, a small amount of normally distributed random noise is added onto each resampled data. This is equivalent to sampling from a KDE of the data by Formula 3.8.17:

$$ \hat{f}_h(x) = \frac{1}{ng} \sum_{i=1}^{n} K\left(\frac{x - X_i}{h}\right) \qquad (3.16) $$

where h is the smoothing parameter, K is a symmetric KDE function with unit variance, $\hat{f}_h(x)$ is a standard kernel estimator. This sampling process was repeated in a loop several times until the data were sorted and grouped (Fig. 3.8.22).

The initial step included creating unsorted dendrogram. Several adjustments were applied thereafter: changed and colored labels based on real profiles groups category, colored branches based on the procedure of cutting dendrogram tree into the clusters. Then dendrogram was sorted using ML algorithm. After the reflection of the unsorted dendrogram on the sorted according to the data distribution by sample values, the model proceeded to the comparison of the two dendrograms (Fig. 4, [434]). Compared dendrograms formed cluster kernels. Further processing was completed on the machine by an embedded mathematical algorithm of R that filter out less significant profiles and compares them to the kernel ones. The data are stored to make final step: clustering using bootstrap probability. Thus, there are 4 modifications of the initial model consisted of cross-section bathymetric profiles of the trench.

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

<figure><img src="image_118.png" /><figcaption>**Figure 3.8.21:** Mariana Trench: hierarchical dendrogram tree of the profiles, by R code: A.6.6</figcaption></figure>

Hierarchical dendrogram clustering method was selected due to its flexibility: this type of the ML is used specially for unlabelled data, i.e. without defined categories or groups, aimed to understand data possible allocation into groups. As the structure of the tectonics and trench geomorphic properties are rather complex, the dendrogram clustering facilitated data partition through defining the resulting groups, with similar morphology, distance to the volcanic areas, slope angles and sediment thickness. Used R script for cluster analysis: A.6.6.
<figure><img src="image_119.png" /><figcaption>Figure 3.8.22: Mariana Trench: hierarchical clustering, R code: A.6.6</figcaption></figure>

To summarize methodology of clustering, the advantages consist in usability of R libraries and algorithms at one time to understand the best option for the data set. The approach is dendrogram by the {dendextend} library aimed at data partition according to the geospatial properties [402]. A table in (.csv) format was read as data frame into R, and a script was called based on mathematical algorithm of data processing. It calculated variables (geology, bathymetry, geomorphology, sedimentation and tectonics), and divided the dataset into clusters according to their similarity. Once the data are analyzed by the ML, the visualization of the groups and data partition were performed.

### 3.8.18 Scatterplot Matrices

Scatterplot (covariance) matrix is an advantageous approach for the analysis of the big datasets containing multi-source data. Covariance matrix is a square matrix giving the covariance between each pair of elements in a dataset. It illustrates determinant for linear correlation between multiple variables, are efficient way to highlight the environmental variables that might have similar correlations to the geomorphic structure of the trench. Mathematical definition of the covariance matrices is given in Formula 3.9.6.9:

$$K_{XX} = cov[X, X] = E[X - \mu_{X}](X - \mu_{X})^{T}] = E[XX^{T}] - \mu_{X}\mu_{X}^{T} \quad (3.17)$$

The matrix is used to solve the problem of selection of the crucial impact factors affecting the hadal formation. The application is supported by purely practical algorithms of data partition that is, inspection of the possible determinants. The data with multiple variables congregating a geo-system are particular suitable to be used for the scatterplot matrix. The matrices correlation for tectonic plates were executed by {Ggally}, {data.table} and {ggplot2} libraries.
<figure><img src="image_120.png" /><figcaption>Figure 3.8.23: Mariana Trench: scatterplot matrix of the correlation between variables. R.</figcaption></figure>

The matrix of Pearson correlation coefficients between each of the random variables in a given dataset, closely related to the covariance matrix is defined by Formula 3.8.18:

$$corr(X) = (diag(K_{XX}))^{-\frac{1}{2}} K_{XX} (diag(K_{XX}))^{-\frac{1}{2}} \quad (3.18)$$

where diag(K$_{XX}$) is the diagonal matrix X$_{i}$ of the variances (diagonal elements) of K$_{XX}$ for i = 1, ..., n. The R script, based on this mathematical definition, was used to assess correlation between factors affecting trench geomorphology: A.6.41. The bathymetric determinant of the 4 distinct groups of the cross section profiles cause certain variability in the sediment thickness of the basement, slope angle steepness degree, angle aspect, depth at the basement and statistical values (means, median and absolute minimum). A magmatism of the nearby area is notable among other factors. The closeness of the igneous volcanic areas contribute towards earthquake frequency across tectonic plates. The results demonstrated that cross-section profiles of the trench can be divided into clusters according to their properties, e.g. bathymetry, geographic location (latitude and longitude), geological (width of the sediment thickness), magmatism (large igneous polygon areas) and tectonic (location of the tectonic plates nearby).

The algorithms consist of taking variables to show correlation of their values with others across the range of the dataset. Scatterplot matrix shows similarity at the diagonal view, with increasing dissimilarity as the differences in the geologic and tectonic values increase. It analyzes correlation between the factors. Correlation function shows similarities between the characteristics of the tectonic plates. When depth increases (Fig. 1, [430], right center), correlation with sediment thickness increases to a maximum.
As sample points change their position from plates, correlation increases to a maximum, because geo-logic settings vary accordingly. A notable peak in the correlation between the magmatism of the igneous volcanic areas and slope°indicates the dependancies between the geomorphology and tectonics.

In a resulting covariance mapping the values of the cov(X,Y) or pcov (X, Y | I) matrix are plotted as a 2-dimensional map (Fig. 3.8.23) by following Formula 3.8.18:

$$[X_1, X_2, \dots, X_n] = \begin{bmatrix} X_1(t_1) & X_2(t_1) & \dots & X_n(t_1) \\ X_1(t_2) & X_2(t_2) & \dots & X_n(t_2) \\ \vdots & \vdots & \ddots & \vdots \\ X_1(t_m) & X_2(t_m) & \dots & X_n(t_m) \end{bmatrix} \quad (3.19)$$

The method of the scatterplot matrix shows appropriate visualization for the large dataset overstep-ping thousands of the observation points. Scatterplot matrix displays values for typical variables for a dataset: depths, aspect, slope° of the transects, sediment thickness, volcanic areas, topographic settings (max/min/mean depths, by profiles). The diagonal view of the scatterplot matrix gives an inside in to the relationships between these factors. Based on the scatterplot matrix it can be roughly determined that there is a linear correlation between several geological variables: sediment thickness have correlations to the slope angle degree and a depth values.

Scatterplot matrix enables to visualize several approaches (Fig. 3.8.23): scatter plots, boxplots with residuals, margin plots and surface plots, R code A.6.41. The effects of factors that increase sedimentation (location on tectonic plates, steepness) may hidden other environmental triggers underlying sedimentation responses. This is also induced by the biomass accumulation in the deep layers of the trench. Changes in patterns of the deep currents and intensity affected by submarine volcanism may cause increase sedimentation. This proves a feedback mechanism between the marine geological and biological factors. The submarine sedimentation in the oceanic trenches spatially correlate to the geological and tectonic factors and to the deformation features that induce submarine volcanism and closeness of the igneous volcanic areas. In such cases, sedimentation response rapidly.

The results indicate (Fig. 3.8.23) that sediment sources at the Mariana Trench are connected with trench tectonics. The relationships and links between the sediment thickness, slope angle steepness and profiles depth of the trench are shown on (Fig. 3.9.11. The sediment thickness gives a response towards the steepness and the depth. The rate of the submarine sediment accumulation is critical due to the preconditioning by these factors. It is also sensitive to the triggering factors: ocean deep currents contributing to the sedimentation and topographic variability observed in the cross-section profiles. R proved to be an effective tool for modelling scatterplot matrices that show distribution of the factors affecting geomorphic properties of the oceanic trench by distinct unevenness of the factors [609].

### 3.8.19 EFA

Factor analysis and PCA are based on R script: A.6.26 using dataset [431]. The fa() function from R was used (A.6.26), with the following primary arguments: r: the correlation matrix; Nfactors: number of factors to be extracted; rotate: one of several matrix rotation methods, such as varimax or oblimin. R libraries for scripting simulations were used to model principal factors affecting the trench morphology. Mathematical definition of the statistical model is as follows.

There is a given set of p observable random variables, $x_1, \dots, x_p$ with means $\mu_1, \dots, \mu_p$, unknown constants $l_{ij}$ and $\kappa$ unobserved random variables $F_j$ (common factors, as they influence all the observed random variables), where $i \in 1, \dots, p$ and $j \in 1, \dots, k$ where $k < p$. Then the definition can be described by Formula 3.8.19:

$$x_i - \mu_i = l_{i1}F_1 + \dots + l_{ik}F_k + \varepsilon_i \quad (3.20)$$