Seafloor and Ocean Crust Structure of the Kerguelen Plateau from Marine Geophysical and Satellite Altimetry Datasets

Abstract

The volcanic Kerguelen Islands are formed on one of the world’s largest submarine plateaus. Located in the remote segment of the southern Indian Ocean close to Antarctica, the Kerguelen Plateau is notable for a complex tectonic origin and geologic formation related to the Cretaceous history of the continents. This is reflected in the varying age of the oceanic crust adjacent to the plateau and the highly heterogeneous bathymetry of the Kerguelen Plateau, with seafloor structure differing for the southern and northern segments. Remote sensing data derived from marine gravity and satellite radar altimetry surveys serve as an important source of information for mapping complex seafloor features. This study incorporates geospatial information from NOAA, EMAG2, WDMAM, ETOPO1, and EGM96 datasets to refine the extent and distribution of the extracted seafloor features. The cartographic joint analysis of topography, magnetic anomalies, tectonic and gravity grids is based on the integrated mapping performed using the Generic Mapping Tools (GMT) programming suite. Mapping of the submerged features (Broken Ridge, Crozet Islands, seafloor fabric, orientation, and frequency of magnetic anomalies) enables analysis of their correspondence with free-air gravity and magnetic anomalies, geodynamic setting, and seabed structure in the southwest Indian Ocean. The results show that integrating the datasets using advanced cartographic scripting language improves identification and visualization of the seabed objects. The results include 11 new maps of the region covering the Kerguelen Plateau and southwest Indian Ocean. This study contributes to increasing the knowledge of the seafloor structure in the French Southern and Antarctic Lands.

1. Introduction

1.1. Background

In cartography and spatial data processing, the task of plotting maps (also known as mapping layouts) is widely used as a common practice. This involves the visualization and representation of spatially defined objects using cartographic techniques [1,2,3]. One of the generally accepted methods for mapping and quantitative and qualitative cartographic visualization is implemented through Geographic Information Systems (GIS), using algorithms for raster and vector data processing embedded within these programs [4,5,6]. However, the high complexity and time-consuming nature of GIS are not conducive to large-scale mapping applications in various Earth science tasks. To address this, programming and scripting algorithms are proposed as an alternative to the GIS-based approach, aiming to reduce data processing complexity in computer graphics and time costs associated with mapping [7].

The application of programming methods in cartography has seen extensive development, resulting in the release of several programs that utilize scripts as mapping tools. Programming-based cartographic approaches can be categorized into three general types. First, there are partial uses of scripts, such as plugins or alternative tools in addition to the existing Graphical User Interface (GUI) in classical modes. Examples include ESRI’s ModelBuilder [8], modules in the Geographic Resources Analysis Support System (GRASS GIS), processing script editors in QGIS or ArcGIS [9,10], and the Python-based ERDAS Macro Language (EML) used in Erdas Imagine to create spatial models. Second, there are spatial libraries of programming languages, including selected packages in R and Python specifically designed for satellite image processing and geospatial data analysis [11,12,13,14]. Third, there are programs that are completely based on using scripting languages without any Graphical User Interface (GUI), such as the Generic Mapping Tools (GMT). All these examples of using scripts aim to automate cartographic data processing using advanced scripting tools.

The script-based cartographic programs are founded on the principle of utilizing the syntax of a programming language, which includes a number of key commands and recognizable expressions by the system [15]. Following the rules of the embedded language, it becomes possible to compose scripts for data modeling and cartographic visualization. In contrast to traditional GUI-based GIS methods, cartographic scripts do not directly generate maps. However, they provide a series of commands that encompass crucial information about the map’s appearance, governing specific features on a plotted cartographic layout [16]. In fact, scripts and programming commands used in cartography define the elements present on the map produced by the script during execution. Specifically, it is possible to define key map concepts such as mathematical definitions of the map (projections, resolution, grid, coordinate systems, extent, and scale), design of symbols and legends (colors and object sizes, palettes for continuous fields, and transparency), and exposition of the elements (overlay, topology, generalization) [17].

1.2. Problem Formulation

The analysis of the seafloor structure and ocean crust concerns mapping, detecting, and recognising the bathymetric structures and variations in geophysical fields. Among these problems, the integrated geophysical interpretation of the satellite altimetry, marine gravimetry, and magnetic anomalies, as well as acoustic and seismic surveys present one of the most active research areas that have attracted research interest in recent decades [18]. Previous studies [19] provided a broad overview of approaches for analysis of the satellite gravimetry data to model the Earth beneath the sea. Marine geophysical and satellite altimetry data provide information on crustal density and processes in the upper mantle [20,21,22], sediment thickness and basement depth [23]. Moreover, the analysis of the geophysical data enables to reveal key characteristics regarding the lithosphere thickness such as Moho discontinuity and elastic features of the lithosphere [24], viscosity and flexural rigidity. Such data are essential for investigations on the Earth’s structure.

In terms of the data-driven analyses, previous works investigating the oceanic seabed can be roughly classified into studies focused on the bathymetric mapping and geophysical analysis. Mapping the seafloor characterizes the bathymetric patterns using data visualization and cartographic methods [25,26,27]. Marine geophysical methods investigate the geology of the continental margins using methods of deep-sea drilling [28,29], the analysis of features extracted from seismic survey and remote sensing data [30], investigating the structure of the deep ocean basins [31,32,33], and modelling the mid-ocean ridge system [34,35,36]. Although traditional methods of bathymetric survey supply novel information on origin, evolution, structure, stratigraphy, and tectonic features of the oceanic crust, they require costly hydrographic equipment such as mutlibeam echo sounder systems [37].

1.3. Related Work

Seafloor bathymetry and geophysical setting of the oceanic crust are two important topics for modelling lithosphere [38]. However, the integrated capture and processing of the multi-source data—such as bathymetry, marine geologic data (development on the Earth’s crust including age, spreading rates and symmetry) and marine geophysical data (sediment thickness, gravity and magnetic anomalies)—has always been a challenging task, and few attempts have been made to explore it. Since the variety of data sources is increasing in modern cartography, it results in an exponential growth in the volume of data serving cartographic applications [39]. Therefore, the use of the remote sensing data, satellite altimetry and gravimetry has proven to be effective in geophysical studies. Numerous examples exist in the literature that focus on estimating the geoid values, evaluating gravity anomalies, analysis of the tectonic and crustal structures, glacial and hydrological modelling and habitat mapping.

Satellite altimetry data can be used to jointly represent the patterns of the oceanic currents and to draw conclusions on connectivity between the habitats [40]. For instance, a study by [41] uses the hydrographic and acoustic data for the analysis of the vertical layers of the ocean to identify the community distribution. An example of the tide modelling of bathymetric gradients was presented using the satellite altimetry data [42]. The retreat of glaciers using mass balance measurements was estimated using remote sensing data to analyze spatial distribution of the surface mass balance [43]. Other studies used the Gravity Recovery and Climate Experiment (GRACE) mission for detecting trends and variations in the ice-sheet mass balance by evaluating gravity signals [44]. Furthermore, Mathieu et al. [45] combine the SRTM DEM, remote sensing data and geological sampling to clarify seafloor structures.

These and similar examples illustrate that a combination of bathymetric and geophysical data for the detailed analysis of Earth’s crust structure and topography outperforms their use separately. This suggests that seafloor mapping should be based on the complex analysis of the heterogeneous features of the seabed. Such information can be derived from various Earth observation datasets, as pointed out earlier [46,47,48,49]. Since such methods can be regarded as feature-level combinations of seafloor features showing the local appearance of prominent seabed elements such as fracture zones, ridges, seamounts, deep-sea trenches, canyons, or rifts, the use of a programming approach to merge different integrated geophysical datasets enables matching multi-source data with varying resolutions and origins. Such data can be used for the detailed analysis of seafloor structure and geodynamic processes [50,51,52].

1.4. Objection and Motivation

In this study, several bathymetric and marine geophysical datasets were used for the analysis of the Kerguelen Plateau (Figure 1) within the south-western segment of the Indian Ocean to visualize, describe and analyze the structures of the seafloor in the context of the geophysical settings. Ten new maps are presented and described to visualize spatial variations, reveal correlation and matches between the geophysical and topographic setting of the Kerguelen region to continue the existing similar studies [53,54,55]. Given the remote location of the Kerguelen Archipelago, one of the most isolated places on Earth, and associated difficulties and high cost of geological and geophysical sampling, the integrated use of the satellite-derived data enables to us have better insight into the structure of the Kerguelen Plateau.

2. Study Area

The Kerguelen Plateau presents a large topographic elevation in the south-west Indian Ocean, extending over 6500 km${}^2$ [56], Figure 1. The Kerguelen Archipelago, formed on the plateau, constitutes a shallow submarine plateau. The structure of the Kerguelen Archipelago is asymmetric with notable variations in the two distinct parts: the northern part (Heard, McDonald and Kerguelen Islands) is more shallow (<1000 m) compared to the southern segment (depths of 1500–2000 m) [57]. Morphologically, the plateau is oriented in a NNW direction toward the Antarctic continental margin where it is constrained by Princess Elizabeth Land. On the northwest, the archipelago is limited by the Crozet Archipelago, located 1250 km away [58]. To the north, it is bordered by the African-Antarctic Basins; to the northeast, by the Australian-Antarctic Basin; and to the south, by the Antarctic. Adjacent lands include Enderby Land, Kemp Land, Mac Robertson Land, Princess Elizabeth Land, Wilhelm II Land, and Queen Mary Land, all located in the Antarctic, as shown in Figure 2.

The origin of the Kerguelen Plateau is related to the hotspot associated with the Gondwana breakup in the Cretaceous period (120–110 Ma), and seafloor spreading between India and Antarctica in the south-western Indian Ocean [59]. The initial surface manifestations of the Kerguelen Plume in the southern region began during the Mesozoic era around 120 Ma [60]. Active submarine volcanism occurring on the aseismic ridge, coupled with sedimentation processes since the Cretaceous, led to the creation of the basaltic basement of the Kerguelen region [61]. This volcanic activity persisted within the Kerguelen hotspot until the Tertiary period (<40 Ma) and is evident in the geological and mineralogical composition of the archipelago. The archipelago is covered by up to 85% basalts originating from the mantle plume of the hotspot. These basalts result from the crystallization of raised and cooled basaltic magma [62,63,64]. Presently, the Kerguelen Plateau stands as one of the largest volcanic plateaus globally, spanning a total length of 2300 km. Remnants of volcanic activity are still visible on the Heard and McDonald Islands [65].

The Kerguelen Plateau belongs to the French Southern and Antarctic Lands, known for its protected environment [66,67,68]. Its high environmental value is explained by the unique wildlife structure, which includes vulnerable species. The remote location of the archipelago—4000 km from South Africa and Australia [69,70]—has created ideal conditions for preserving unique flora and fauna of the Kerguelen Archipelago. The presence of rare Antarctic plant species and high biodiversity is closely linked to the soil developed on the basaltic basement and the specific geological substrate of volcanic origin. Additional factors contributing to this uniqueness include the influence of the Polar climate [71,72,73]. The distribution of fish communities and phytoplankton is affected by seasonal changes of the Antarctic circumpolar currents [74,75,76,77,78,79,80,81]. Furthermore, the steep bathymetric slopes and the exposure of the Kerguelen plateau amplify the speed and intensity of the ocean currents’ circulation [82]. The combination of all these factors makes the ecosystems of the Kerguelen Archipelago unique and deserving of protection as an environmental heritage of the French Southern and Antarctic Lands [83,84,85].

3. Materials and Methods

All the maps have been plotted using the scripting toolset Generic Mapping Tools version 6.4.0 [86], https://www.generic-mapping-tools.org/. A key aspect of the GMT algorithms is that they consider each cartographic element by its parameters and add new layers on the map regarding the target location, which can be adjusted using refined flag options in special modules. Thus, a modular scripting approach of GMT distinguishes it from GIS.

3.1. Data

The materials used as input cartographic grids include the following datasets. The marine geological data on age, spreading rates, and asymmetry of the ocean crust were derived from the NOAA high-resolution dataset [87], Figure 3. The units of the dataset on spreading asymmetry are expressed as % of crustal accretion of the seafloor, i.e., symmetric spreading results in values of 50% on conjugate ridge flanks. The bathymetric mapping was based on the ETOPO Global Relief Model [88]. The gravity grids were obtained from the EGM-2008 and EGM96 [89]. Gravity data are grounded on the concept of gravity anomaly across different Earth surfaces, which correlates with the interplay between topography, mass distribution within local and regional relief structures, marine gravity values, and the establishment of gravitational potential across various Earth locations.

The free-air and Bouguer gravity anomaly grids are derived from high-resolution grids [90,91]. Specifically, the gravity field varies significantly over the oceans due to density fluctuations. The most pronounced anomalies emerge from fluctuations in density, such as those occurring at the heterogeneous seafloor or at the crust-mantle interface (Moho discontinuity) [92]. Such data facilitate the modelling of gravity fields for the analysis of Earth’s structure. The latest updates in gravity grids enhance the precision of altimetry-derived gravity anomalies. While gravity grids do not directly indicate topography, they offer crucial insights into the Earth’s relief, providing an approximate relationship between topographic and geophysical data and illustrating the correspondence of surface gravity values with local geoid undulations [93].

The data on sediment thickness are collected from the GlobSed dataset, which models distribution of the sediment thickness in the World’s Oceans [94]. The magnetic data are derived from the two available information sources compiled originally from satellite and marine magnetic measurements: the Earth Magnetic Anomaly Grid (EMAG2) [95], which has a resolution of two arc minutes, and the World Digital Magnetic Anomaly Map (WDMAM) with a resolution of three minutes [96]. The gravity data were derived from the free-air global gravity grids [97].

The analysis of the geodetic data is useful in such complex terrains since it facilitates the identification of major tectonic structures and topographic patterns. In addition, the EGM grids, the identification of lithospheric structure and tectonic blocks can also be performed using satellite gravity data such as World Gravity Model (WGM) and GRACE [98]. Furthermore, the EGM and astrogeodetic vertical deflections are useful for modelling gravity and geopotential differences [99]. Finally, processing terrestrial and satellite geodetic data can be applied to address theoretical geodetic challenges, such as the altimetry–gravimetry boundary-value problem [100]. Other instances of the implications of geodetic datasets encompass the utilization of satellite-derived data from Cryosat-2 and altimetry in conjunction with ship-measured gravity to estimate marine gravity anomalies [101].

3.2. Methods

In this study, the cartographic methodology is based on using the Generic Mapping Tools (GMT) programming suite version 6.4.0 [102] by the developed scripting workflow, explained in detail in earlier works [103,104]. The most prominent feature of the GMT is a scripting approach that principally distinguishes it from the conventional software due to the embedded programming language [105]. The traditional GIS-based methods either employ the GUI for mapping with existing standard menu or allow scripting as a complimentary workflow.

In contrast, the GMT is a completely console-based software that operates entirely using scripts. In this way, it captures the data by running a script written using the embedded syntax that operates similarly to programming. A script consists of a sequence of predefined commands with parameters that control the appearance of cartographic elements and features. The most well-known cartographic tools that employ scripts for data processing are GMT [106] and GRASS GIS [107]. Other software allows scripts as an additional functionality alongside the standard GUI, as seen in ArcGIS or QGIS. While the latter is predominantly used for image processing, cartographic, and ecological studies [108,109,110,111], GMT finds application in geophysical and geological research [112,113,114].

3.2.1. Topographic/Bathymetric Mapping

The script for plotting the topographic/bathymetric map is provided in Listing A1. The code’s most crucial elements are as follows. The mapping process initiates by selecting the study area from the global ETOPO1 grid using ‘grdcut’. Isolines are chosen through ‘grdcontour’ with a 1,000 m gap. In our case, these lines denote significant locations that offer a high representation estimate for seafloor bathymetry. ‘grdimage’ is utilized to generate an image plot of the raster grid, colored as defined by the ‘makecpt’ module, with values -T indicating the extreme (min/max) data for the topographic grid, as extracted by the ‘gdalinfo’ module. The remaining code details are outlined in Listing A1 along with key comments.

3.2.2. Mapping Seafloor Age, Spreading Rates, Spreading Asymmetry and Age Uncertainty of the Ocean Crust

Mapping the age, age uncertainty, spreading rates, and spreading asymmetries of the Kerguelen Plateau is executed using raster grids with a 2-min resolution, as shown in Figure 4. These grids are reprojected to the Lambert Azimuthal Equal-Area projection utilizing data from major features within the Indian Ocean basin. The ‘psxy’ module is employed to add vector lines representing mid-ocean ridges, tectonic plates, and lithospheric plates. This provides a depiction of the borders of the Indian, African, and Australian plates, as seen in Figure 4. Spreading half rates of the ocean crust are visually presented using the 2-arc-minute netCDF grid v.3.6 from NOAA. Similarly, tectonic slab contours are added via the ‘psxy’ module of GMT. Multiple grids are combined and displayed as plotted raster grids to illustrate age-depth relationships in the seafloor around the Kerguelen Plateau and southwest Indian Ocean. The complete scripts are available in Appendix A of this study, with Listing A2 utilized for mapping the age of the oceanic lithosphere over the Kerguelen Plateau, Listing A3 for mapping asymmetries in crustal accretion on conjugate ridge flanks over the Kerguelen Plateau, Listing A4 for visualizing spreading half rates of the oceanic lithosphere in the Kerguelen Plateau and SW Indian Ocean, and Listing A5 for mapping crustal age uncertainty of the oceanic lithosphere over the Kerguelen Plateau, as depicted in Figure 5.

3.2.3. Mapping Sediment Thickness

A representation of the GlobSed 5-arc-minute total sediment thickness grid is extracted from the global raster grid on a target location of the Kerguelen Archipelago area and surrounding regions using ‘grdcut’ command, which selects the file. The amplitude of values is checked using the Geospatial Data Abstraction Library (GDAL) library embedded in GMT using the ‘gdalinfo’. Based on the range of the values, the colour palette table is defined accordingly and presented using the ‘psscale’ module. This module is configured with parameters specifying the geographic location for each element on the map. The ‘-Dg’ command is employed to define the position of the colour scale on the map using coordinates. Furthermore, grid annotations and graticules are added utilizing the ‘psbasemap’ GMT module. The remaining part of the script, complete with added comments for concise explanations, is outlined in Listing A6 within the Appendix A of this study. This particular listing was utilized to create the map featured in Figure 6.

3.2.4. Geophysical Mapping

To emphasize features and subsurface characteristics that might not be apparent from geological data alone, maps of gravity anomalies were generated using gravity datasets projected onto the Lambert Equal-Area Azimuthal Projection for consistency with other maps. The complete scripts for creating these geophysical maps can be found in Appendix A, with Listing A7 detailing the process for free-air gravity anomalies and Listing A8 for vertical gravity gradient. The crucial lines of code are summarized below.

For free-air gravity anomalies, the ‘img2grd’ module was initially employed to convert the image from the .IMG format into GRD format. The projection was defined using the ‘-JU43/6.0i’ flag within the ‘grdimage’ module. The extent of the region was specified using the ‘-R’ command, and this information was subsequently passed to subsequent modules using the ‘-R -J’ flags.

In the case of the vertical gravity gradient map, a satellite-derived grid sourced from altimeters CryoSat-2 and Jason-1/2 [115] was used. This dataset facilitated the identification of ridge propagation on the seafloor.

Due to GMT’s ability to process each feature only once and utilize relevant modules (‘psbasemap’, ‘grdcontour’, ‘psscale’), it boasts a high processing speed. This technical advantage empowers GMT to efficiently handle grids encompassing extensive geographic areas in a matter of seconds. With an understanding of the seafloor’s geological and geophysical structure, GMT demonstrates rapidity and proficiency in mapping.

By adjusting projection parameters, the algorithm concentrates on key seafloor features within the desired scale and projection. Additionally, the GMT-based framework excels in determining the location, scope, and boundaries of mid-ocean ridges and significant basins. Through cartographic generalization involving the gap in isolines, redundant bathymetric features can be eliminated. This process identifies pertinent and nonrepetitive features, underscores correlations between geophysical settings, gravity, and magnetic anomalies, and facilitates their comparison with geodetic anomalies. Consequently, GMT offers an effective approach to map and visualize seafloor structures.

3.2.5. Mapping Magnetic Anomalies

To map the magnetic anomalies, two different grids were used—the WDMAM and the EMAG2. The first was a subset from the global grid embedded in the GMT using the ‘grdcut’ module, and the second was obtained from the existing raster file of EMAG 2-min resolution grid. As a result, GMT offers flexibility in data processing through two approaches: utilizing pre-existing stored files and accessing embedded grids available within the system remotely. Subsequently, the images were rendered using the ‘grdimage’ module, accompanied by color palettes tailored to the data distribution range across the Kerguelen Plateau and the adjacent southwestern Indian Ocean. The ‘makecpt’ module was employed for this purpose. For a comprehensive view of the GMT scripts for both maps visualized with WDMAM and EMAG2, refer to Listings A9 and A10 in Appendix A.8.

3.2.6. Mapping Geoid Anomalies

The map of geoid inundations was created using the EGM96 dataset through the GMT script, with each module defining specific elements on the map. To achieve effective mapping, the scripts were executed from the console, tailoring the geoid plotting to target concepts such as projections, grids, colour palettes, raster extents, and annotations. The stepwise definition of these cartographic elements is evident on the geoid maps, as demonstrated in the script provided in Listing A11. This approach to data visualization ensures high levels of automation in data processing through GMT, surpassing conventional tools.

Significant features such as geoid undulations, seafloor isochrons, basement depth isolines, and asymmetries in crustal accretion are influenced by the asthenospheric flow from Kerguelen’s mantle plumes to the spreading Southwest and Southeast Ridges. Mapping these data reveals insights into seafloor attributes and enables comparisons with the geological structure and geodynamics. Consequently, mapping multi-source topographic and geophysical datasets allows for the comparison of features using consistent projections. However, it is important to note that real-world seafloor features are seldom identical due to various factors impacting bathymetry, although they may exhibit strong correlations. To address this, GMT scripts can be adapted to different scales to map features in more detailed views. In order to assess the influence of geophysical and geological settings on seafloor structure using GMT’s capabilities, numerous cartographic re-projections were generated. Each map represented selected seafloor features at varying levels of detail. The Lambert azimuthal equal-area projection was ultimately chosen as the optimal projection for its applicability across all maps, ensuring comparability and coherence.

4. Results and Discussion

Automatic matching of topographic and geophysical grids with high accuracy is essential for complex geologic–tectonic investigations. This paper demonstrated the use of scripting algorithms of GMT. These algorithms were applied for plotting eleven thematic maps covering the south-west segment of Indian Ocean and East Antarctic using bathymetric, geodynamic and high-resolution geophysical datasets on gravity (Figure 7 and Figure 8). Cartographic scripts by GMT, as demonstrated in this research, provide visualized information on the geodynamic and geophysical setting of remotely located areas such as Kerguelen Plateau. During the cartographic process and workflow, the scripts enable us to save time through the increased speed of mapping due to the high level of automation.

The programming concept of GMT enables us to better tune and adjust the layout of the cartographic plots in various scales and focus on a specific area for comparability of maps in a series. Such compatibility facilitates the evaluation of correlations among various geophysical and bathymetric features that have developed over the extensive history of the Kerguelen Archipelago. Furthermore, by utilizing geospatial analysis as a complementary technique, it becomes possible to compare and analyze the relationships within the geodynamic setting, such as age, spreading half rates of the ocean crust, asymmetry in crustal accretion on conjugate ridge flanks, and other variables. In this manner, GMT scripts offer an advanced cartographic method for visualizing datasets and extracting information through efficient data processing and modelling of geophysical properties of the seafloor.

These maps reveal the details of the structure of the seafloor in the Kerguelen Plateau. The results confirm that the dynamics of seafloor development, as reflected in the maps of oceanic crust age, asymmetries, and spreading rates of the south-west and southeast Indian ridges, are closely related to the major geophysical setting as depicted on the topographic, magnetic, and gravity grids. Additionally, the volcanic activity of the Kerguelen hotspot has a significant impact on the distribution of magnetic anomalies, which aligns with previous studies [116,117,118]. Overall, the results demonstrate that GMT scripting is a powerful and stable cartographic method that efficiently performs geophysical and bathymetric seafloor mapping. In the sections below, we discuss the obtained results on relevant maps, providing comments on the essential features and characteristics of the seafloor around the Kerguelen Plateau, south-west Indian Ocean, and East Antarctic.

4.1. Ocean Floor Formation

The advanced methods of visualization by GMT constitute an important element of the map content through detailed plotting of the depicted objects, which enables us to indicate qualitative and quantitative geophysical specifications for reference and analysis. Thus, the general physical-geographic structure of the Kerguelen Plateau is visible in Figure 1, which shows a morphological orientation of the archipelago in the NW-SE direction and an extent surpassing 2000 km in length. As can be seen, the northern and southern parts of the plateau are asymmetric, where the less expressed southern part is older and lies in deeper water in the topographically downlifted areas. Age, spreading rates, and spreading symmetry of the ocean crust indicate the gradual evolution steps of the ocean floor formation in the southwest Indian Ocean and around the Kerguelen Plateau. The geodynamic setting of the oceanic crust (Figure 2, Figure 3, Figure 4 and Figure 5) shows a strong relationship between the Kerguelen Plateau and the two mid-ocean ridges, which can be revealed from the analysis of the relief.

Relief is the main element of the seafloor since it reflects its geological structure and geodynamic history. Accordingly, the relief of the seafloor surface around the Kerguelen Plateau forms a continuously changing field of bathymetric heights. There are also sharp changes in altitude around the archipelago and mid-oceanic ridges. To depict the relief, GMT enables the modeling of isolines using the ‘grdcontour’ module and adjusted color gradients using methods of qualitative coloring of background and gradients according to the actual heights. At the same time, there are specific benefits of GMT techniques for mapping the hypsometric maps. Thus, the quantitative values of the relief make it possible to obtain absolute heights and elevations from the raster grid; the characteristics of the curvature and steepness of inclination can be obtained using the GMT module ‘grdtrack’ through cross-sectioning [119,120,121]. Moreover, GMT enables the modeling of the plasticity of the relief, that is, to depict a nonlinearity of the landform irregularities that form a visual image of the submarine terrain. This enables an analysis of the morphological conformity of the relief, which highlights major seafloor features, specific landforms, and their structure.

The age of the ocean crust was determined by interpolating the adjacent seafloor isochrons oriented towards the direction of seafloor spreading, as shown in Figure 2. This correspondence highlights the unique geophysical setting of the Kerguelen Plateau, underlain by the oceanic crust, which is strongly associated with its tectonic origin associated with volcanic hotspot and geologic history. The formation of the Southwest and the Southeast Indian Ridges is related to the uplift of the Kerguelen Plateau as a remnant of the Mesozoic oceanic basin existing after the separation of Gondwana. The comparison of the bathymetric map with geodynamic maps shows that seafloor heterogeneity around the Kerguelen Archipelago correlates with the seafloor spreading rates, where rougher basement is formed in the areas with the low half-spreading rate threshold (30–35 mm/year and lower). This correlation can be revealed by comparing the maps in Figure 1 and Figure 4. Such heterogeneity in seafloor patterns varies significantly in various basins of the Indian Ocean, depending on the geodynamic setting and the geologic development of the oceanic crust.

The gridded map of age of the ocean crust around the Kerguelen Plateau (Figure 2) shows a correlation with the observed spreading half rates of the lithosphere (Figure 4). Originally formed as a single structure, Kerguelen was then split by the seafloor spreading in the south-west sector of the Indian Ocean which resulted in the formation of the two segments of the archipelago [122] which are visible on the maps. Furthermore, the difference in volcanic activity between the northern and central parts of the Kerguelen Plateau, underlying Heard Island, indicates that it is located on a hotspot, with various parts of the islands experiencing the effects of the mantle plume on different scales [123,124,125]. This unstable position has an impact on the geodynamic patterns over the Kerguelen Islands, leading to a higher uncertainty in the age of the oceanic crust (Figure 5) compared to the adjacent areas. Additionally, the structure of the oceanic crust beneath the Kerguelen Plateau is similar to that beneath aseismic ridges such as the Crozet Rise and the Madagascar Ridge, providing evidence that it originates from active volcanism associated with a hot spot [126].

4.2. Sediment Thickness

Mapping sediment thickness using the 5-arc-minute GlobSed grid relies on the approximated modelling, which highlights the sedimentation trends around the Kerguelen Archipelago, Figure 6.

The analytical map showing sediment thickness (Figure 6) displays patterns and key characteristics of sediment distribution over the seafloor. Using GMT-based techniques of mapping, it is easy to highlight the variations in sediment accumulation in various parts of the ocean using an adjusted colour table and the actual range in the data on sediment thickness. Hence, the level of cartographic details depends on the depth of the analysis with regard to the geologic formation of the seafloor and sediment accumulation. The objects on the map show the main regions of accumulated sediments, the structure and trends in distribution, and special features and properties compared to the closeness of coastal areas. Thus, the distribution of the sediment thickness correlates with the age of the underlying oceanic lithosphere and its latitude, which can be noted by the comparison of maps in Figure 2 and Figure 6. Such correspondence is especially visible for higher values of the sediment thickness near the shorelines of the Antarctic, Amery Ice Shelf (3000–4000 m), Enderby Land (over 4000 m) and the Kerguelen Plateau (2000–3000 m).

A higher level of sediment thickness in these areas may also indicate earlier processes of subaerial erosion that occurred before subsidence and associated sedimentation. The rifting process that took place during the Late Paleogene resulted in changes in the Tertiary sediment facies of Kerguelen, which were influenced by the evolution of the Antarctic environment [127]. Sediments covering the Kerguelen Plateau include pillow lavas, tuffaceous sediments, and marine siltstones that were deposited since the late Miocene [128]. These sediments continue as a thick sequence of Cenozoic sediments (over 5000 m) within the Enderby Basin to the southwest of the Kerguelen Plateau (Figure 6).

The correlations observed between sediment thickness (Figure 6) and the age of the oceanic lithosphere (Figure 2) demonstrate the role of ocean floor formation in influencing the pattern of distribution and accumulation of sediments.

Moreover, the analytical map of sediment thickness reflects smaller features and details compared to the bathymetry of the southwest Indian Ocean. Hence, comparing the map of sediment thickness with the bathymetry enables us to detect associations, for example, high values of sediment thickness in the region of Dronning maud Land in the Antarctic, which can be associated with the effects from the processes of weathering and coastal erosion, factors of higher curvature in slopes and topographic variations in heights in the coastal areas. Other important factors increasing the accumulation of the sediments around the Kerguelen and the adjacent area include glacial processes and the turbidity of the ocean currents. Hence, intense circulation results in the accumulation of the large sediment fields with values over 3000 m.

These effects can be attributed to the Antarctic circumpolar currents that started around the Eocene-Oligocene periods and have continued until the present time. These currents constitute the strongest current system in the oceans, directed clockwise around the South Pole, and they significantly influence the adjacent sub-Antarctic regions, such as Kerguelen [129,130].

The general orientation of areas with maximal sedimentation is consistent with the sediment-filled troughs stretching in a NW-SE direction. These troughs are associated with the overall NW-SE orientation of Kerguelen and the axes of tectonic faults. The high values in sediment thickness around the Kerguelen Plateau, particularly contoured by the ridge isolines along its eastern margin, are associated with depositions resulting from bottom currents directed westwards.

4.3. Free-Air Gravity Anomaly and Vertical Gravity Gradients

The visualized marine gravity field over the Kerguelen Plateau region and the adjacent areas of the south-west Indian Ocean are shown in Figure 7. The comparison of the gravity roughness with the map of the half-spreading rates (Figure 4) and sediment thickness (Figure 6) shows the relationship between the speed of the spreading of mid-ocean ridges and roughness of the seafloor basement. Such phenomena are explained by the effects from the process of mid-ocean ridge formation. Other factors include the associated magma flows, spreading directions in Mesozoic and isochron orientations of the age of the oceanic crust, which affect current bathymetric and gravity patterns [131]. Furthermore, the gravity highs around the Kerguelen Plateau and Heard Island correspond to the maximal bathymetric elevations. These gravity highs indicate the presence of seamounts formed by Miocene basalts erupted during volcanic activity in the southern Indian Ocean. This volcanic activity contributed to the formation of a large igneous province [132].

The GMT-based geophysical maps enable us to determine the location and spatio-temporal structure of gravity phenomena that indicate on geological processes, their mutual relationships, and connections with topography. Such analysis supports the identification of trends and dynamics in seafloor development. It helps obtain quantitative characteristics from geophysical data and estimate both the highest and lowest values in gravity grids. In turn, zoning and classification of gravity variations helps to forecast changes in gravity anomalies over the seafloor of the Indian Ocean. Hence, the analysis of maps shows that the values of the free-air gravity over the Kerguelen are higher than in the surrounding areas and reach up to 80 mGal (Figure 8). In contrast, lower values are associated with the bathymetric depressions and have values of −40 to −60 mGal. Furthermore, the abyssal plain is characterized by the medium values of 0–20 mGal, Figure 8. This well illustrates the existing correlation of the free-air gravity anomalies with the distribution of topographic highs and depressions on the seafloor since they are strongly influenced by a gravitational effect of the distributed topographic masses that are caused by the differences in elevation.

Figure 8 displays the mapping outputs for the vertical gradient over the Kerguelen and the southwest Indian Ocean, showcasing the effects of different locations. The visualized map demonstrates a crucial property of gravitational systems, such as free-air gravity, which is not only subject to the effects of geographic location and the latitude of the selected measurement regions but also the altitude of the Earth’s surface. This is because greater altitude implies a greater distance from the Earth’s center, which in turn affects gravity values. Moreover, the vertical gravity gradient identifies variations in gravity with changes in topographic elevations, as depicted in Figure 8. The comparison of gravity datasets provides additional information on the distribution of major geological and seafloor structures, considering the variations of geophysical fields.

Furthremore, the comparison of the vertical gradient and free-air gravity map (Figure 7 and Figure 8) with the topographic map (Figure 1) illustrates the effects of seafloor structure and the distribution of the oceanic bed on gravity, which shows that the highest values correspond to the Kerguelen Plateau and other rises, while lower values are generally associated with topographic depressions.

4.4. Magnetic Anomalies over Kerguelen Plateau

The anomaly of the magnetic intensity at an altitude of 5 km above mean sea level over the Kerguelen Plateau is shown in Figure 9.

Here, high heterogeneity in the geophysical data is related to the past volcanism over the Kerguelen Plateau, including the voluminous basaltic flooding originated from a deep hot spot as an asthenospheric source of mantle plume product. This resulted from the processes of slab dynamics and tectonic plate movements in the southwestern segment of the Indian Ocean. In this regard, combining the data from the WDMAM and EMAG2 (Figure 9 and Figure 10) data on terrestrial gravity fields (Figure 7 and Figure 8) for comparison with maps on oceanic crust development (Figure 2, Figure 3, Figure 4 and Figure 5) presents an integrated GMT-based geophysical analysis.

Crustal volume contributes to the decreased amplitude of the magnetic anomaly around the Kerguelen, as can be seen in Figure 10, with lower values of around −500 mGal. The analysis of the magnetic anomaly patterns in the SW Indian Ocean supports the hypothesis of the spreading seafloor with variations in the oceanic crustal block movements, as reported earlier [133]. This phenomenon is evident from the different magnetic patterns observed over the mid-ocean ridge. Moreover, the comparison of Figure 10 with Figure 5 reveals that seafloor age uncertainties for grid cells coincide with the marine magnetic anomaly identified around the Kerguelen Islands. This correlation is also observed for the conjugate ridge flanks (Figure 3 and Figure 10).

The Earth Magnetic Anomaly Grid (EMAG2) offers the opportunity to assess magnetic-gravity field relationships as descriptors, going beyond the traditional analysis of gravity and magnetic anomalies. Magnetic anomalies arise from geological and topographic features that alter local magnetic fields, making it crucial to comprehend their correlation with geophysical phenomena and topography. For instance, scrutinizing local magnetic anomaly patterns in the southwest Indian Ocean reveals associations with oceanic crust formation, seafloor spreading, and subduction zones. Moreover, the age of oceanic crust and spreading rates, resulting from land accretion and extensive volcanism of the Kerguelen Plateau, are linked to the historical geological development of the region, as mentioned previously [134].

The EMAG-2 and WDMAM grids employed for plotting magnetic anomalies exhibit varying levels of grid detail, allowing insight into the subsurface structure of the seafloor around the Kerguelen Plateau and the composition of the Earth’s crust in the southwest Indian Ocean. The magnetic fabric data correlate with hotspot activities and active volcanism, particularly prominent over the central and northern sectors of the archipelago [135,136,137]. Furthermore, intermediate crustal thickness values within the oceanic crust beneath the Kerguelen Plateau and large-amplitude magnetic anomalies across the archipelago point to the plateau’s oceanic origin, attributed to plate volcanism resulting from tectonic plate activity, as previously reported [138,139]. Therefore, the distinctive magnetic patterns evident in Figure 10 maps correspond to heightened hotspot activity and the associated lava flows.

Deeper masses, asthenospheric upwelling, and mantle plume-driven convection are geodynamic processes that influence the magnetic properties of the Earth’s surface. Moreover, through data analysis, a deeper understanding of the impact of geophysical settings on the distribution of positive and negative magnetic undulations emerges, with the former situated over the Kerguelen Islands and the latter in the eastern regions and southwest of Australia. This analysis enables the assessment of variations in geophysical grids through comparative map analysis. Thus, the cartographic depiction of geophysical and magnetic datasets offers advanced methods for extracting information about seafloor formation and interconnected geophysical processes.

4.5. Geoid Models

The geopotential model over the Kerguelen, based on the EGM96 dataset, is illustrated in Figure 11. The variations in the geoid across the Kerguelen Plateau highlight the ongoing isostatic compensation of the archipelago due to its low-density mantle. Consequently, the high anomalies in the geoid level above the Kerguelen Plateau can be attributed to the significant volcanic activity associated with the formation of ridges on the hot lithosphere. This volcanic activity is reflected in the exceptionally thick crust beneath the Kerguelen, resulting in geoid values exceeding 40 m (Figure 11), surpassing the normal thickness of the oceanic crust. These findings corroborate previous studies investigating the geoid in the Kerguelen Archipelago, which documented anomalous thickness in this region [140]. This isostatic compensation, linked to anomalously high geoid values, corresponds to the rugged elevated terrain in regions experiencing active tectonic uplift. These observations shed light on processes occurring in the upper mantle [141].

A comparison of the geoid map with topography (sf. Figure 1 and Figure 11) implies an existing correlation between the continued geodynamic processes in the south-west Indian Ocean and the topographic structure of the Kerguelen Archipelago. Moreover, this proves a high positive correlation between the geoid height and deep structure of the seafloor topography, as also noted earlier for the regions of large plateaus and swells [142].

5. Conclusions

As demonstrated in this study, the GMT-based mapping approach offers a wide range of cartographic functions for comprehensive spatio-temporal modelling and data visualization through an automated scripting approach. Utilizing GMT for cartographic tasks provides various modules and methods for representing different types of data, making it applicable in diverse fields of geomatics. The advantages of employing GMT in cartographic workflows are manifold. It enables highly automated plotting, facilitating rapid visualization of complex elements and features such as geographical, geological, oceanological, and geomorphological characteristics.

Furthermore, GMT supports multiple formats, encompassing both raster and vector data, and accommodating various classes and types of information. Leveraging the technical capabilities of GMT within cartographic workflows allows for common modelling and basic statistical analysis of spatial data, enhancing the understanding of their properties.

Analyzing a multitude of maps generated using GMT scripts reveals the consistency in depicted objects, facilitating their recognition and interpretation. The broad spectrum of GMT modules integrates scientific and technical methodologies for cartographic visualization and geospatial analysis. This unification aids in feature detection, recognition, and related research support.

Consequently, GMT-based mapping enables the amalgamation of maps for spatial analysis of intricate processes, objects, and phenomena, such as seafloor structures. This approach proves invaluable in addressing scientific and practical challenges within the realms of geophysics and geodynamics. The flexibility of GMT’s syntax, the quality of its cartographic outputs, and its compatibility with various operating systems and computing devices all contribute to its effectiveness.

Given the success and applicability showcased in the executed GMT scripts–characterized by syntax flexibility, high-quality cartographic outputs, and compatibility across different platforms–it is foreseeable that this GMT-based cartographic method can be extended to study seafloor structures in other oceanic regions, considering varying geologic conditions and geodynamic evolutions.

Global surveillance of the seafloor through the use of altimeter satellites and gravity measurements has unveiled significant geophysical anomalies. Integrating data on magnetic field intensity, bathymetry, and deep seafloor geodynamics allows for an evaluation of the interrelations among these processes. This ongoing global surveillance continues to generate extensive high-resolution datasets. However, effectively analyzing these datasets in ever-higher spatial resolutions demands advanced tools for automated analysis. The toolset of GMT scripts has demonstrated its effectiveness in the realm of seafloor bathymetric and geophysical mapping. This approach facilitates the visualization and mapping of diverse seabed features across varying scales and resolutions, aiming to detect correlations between magnetic anomalies, geophysical patterns, and their connections to present bathymetry. Such visualization offers comprehensive coverage of seafloor features across the near-global scale.

Direct seafloor surveys for observations are resource-intensive, involving the use of complex and costly equipment, such as multi-beam echo-sounding systems, for data collection, generation, and storage. Yet, the need persists for efficient datasets that can be readily visualized and analyzed. Utilizing high-resolution geophysical datasets provided by NOAA and USGS, processed through the advanced scripting capabilities of GMT, establishes a cartographic processing pipeline for swift, automated, and accurate seafloor mapping. Notably, GMT’s flexibility plays a pivotal role. A variety of GMT modules can be harnessed to process diverse geospatial data types, catering to different cartographic tasks. This flexibility allows for the adaptation and expansion of the proposed cartographic workflow, addressing larger-scale or smaller-scale mapping needs. The equilibrium between topographic gradients and geophysical grids vividly illustrates the links between seafloor patterns, the structure of the oceanic crust, and the processes within the lithospheric mantle.

From a cartographic perspective, this study underscores that the effective analysis of geologic-geophysical datasets within Earth sciences extends beyond utilizing isolated parameters (e.g., topographic maps) to encompass the selection of multiple datasets. The deep mantle processes, as reflected in geophysical data, intricately shape the seafloor formation. The methodology showcased in this study demonstrates how data from diverse sources (geophysical, topographic, geodetic, geodynamic) can be harnessed for comprehensive cartographic analysis using standardized workflows supported by scripts. This approach unveils additional insights into seafloor variability and the factors influencing ocean crust formation, strongly correlated with topographic patterns.

The cartographic approach exhibited here enables data assimilation and extension, not only across the Indian Ocean but also to other regions of the global ocean. For example, the Pacific Ocean boasts a rich tapestry of seafloor features, including vast abyssal plains, mid-ocean ridges, oceanic trenches, numerous seamounts, and continental shelves. These features present a fertile ground for investigating potential correlations between geophysical and magnetic anomalies and the heterogeneous seafloor patterns. In this context, the application of GMT for seafloor mapping in diverse oceanic regions serves as an ideal scenario for validating the cartographic scripting approach outlined in the methodological framework.

The series of maps presented, along with the comparison of the geophysical settings over the Kerguelen Plateau, underscores the superiority of script techniques over the GIS approach in terms of cartographic workflow automation. The compactness of GMT’s syntax allows for code reusability with modifications. However, there are limitations to consider. GMT necessitates parameter tuning in advance when handling map elements and adjusting projection parameters, as it lacks the ability to preview maps before script execution, being a console-based program. Furthermore, GMT cannot remove redundant features once plotted, requiring the script to be run again for corrections. In contrast, the GIS approach permits real-time adjustments to map layouts, enabling the correction of colour palettes and bathymetric details on the fly. GMT’s console-based nature mandates direct modifications to the script’s code to address cartographic challenges.