Correlations between the Topography-Induced Gravity, Terrain Structure and the Seismicity in the Gulf of Panama

This study presents new maps of the topographic and geophysical setting and seismicity in the region of the Gulf of Panama. The spatial analysis is based on the comparative analysis of the datasets on geoid, free-air gravity anomaly, topography and earthquakes. The cartographic framework is developed using the Generic Mapping Tools (GMT) scripting toolset. Seismic activity in the Central America is high due to the complex geologic setting, tectonic activity and lithosphere plate subduction. The data include the Earth Gravitational Model (EGM2008), the General Bathymetric Chart of the Oceans (GEBCO) and gravity grids. The seismicity data were collected from the Incorporated Research Institutions for Seismology (IRIS) catalogue on 1970–2021. The variations in data were compared to analyse correlations between the geophysical, seismic and topographic parameters. Free-air gravity, geoid and topographic data derived from the high-resolution datasets were used to investigate their effects on the main seismic sources in the region. The comparison of the maps showed that the distribution of the shallow earthquakes in the Pacific segment of Panama coincides with negative free-air anomalies and lower geoid values. The results revealed high values of geoid in the high mountainous regions of Panama (Cordilliera de Talamanca, southern coast of Peninsula de Azuero and eastern Panama, 77.5–78.5°W), which correspond to the topographic roughness in the highlands. Negative values of geoid are found over the Caribbean Sea basin (−4 to 0 m). The analyses of seismicity showed 1740 earthquake events varying by magnitudes from 2.9 to 7.8 at the depths up to 225 m (near the west coast of Colombia). A high concentration of the earthquakes is in the western region of the Panama’s shelf waters (~82–83.5°W), and on the border with Colombia (~77–78.5°W). High gravity anomalies (over 220 mGal) are found in the mountainous regions which match the geodynamic processes associated with the Earth structure and geodetic and geophysical effects. The regions of the high seismicity were defined in the Gulf of Chiriqui and eastern part of the Gulf of Panama.


Introduction
This study presents an application of cartographic scripting techniques for geographical data processing of Panama with a series of new maps prepared using scripts of the Generic Mapping Tools (GMT). The aim is to visualize high-resolution geophysical, topographic and seismic data for analysis of their correlations. Novel cartographic technologies have been rapidly developed over the last decades to support spatial data analysis in Earth sciences. Mapping exploits a variety of the theoretical principles and practical techniques to create Earth models. Geophysical mapping of the seismically active region of Central America and hazard risk assessment require advanced tools of visualization as an essential part of research. As a response to such needs for high-quality mapping, our approach provides a scripting approach of the GMT-based cartographic techniques for mapping seismicity, topography and geophysical grids.
The main technical advantage of GMT scripts consists in the repeatability since they can be reused for any regions of the Earth through a sequence of the codes using GMT modules. As recent papers have shown (Lewis et al., 2008;Eakin et al., 2011;McNamara et al., 2011;Kobayashi et al., 2014), GMT is an excellent instrument for semi-automatic mapping of the geophysical variables. Besides creating maps, GMT enable to manipulate, organize, and display spatial information. Its cartographic functionality includes terrain analysis, topologic analysis, generalization, spatial analysis, 2D and 3D modelling, data formatting and extracting, projecting, gridding, clipping, overlay, logical queries, raster to vector conversion, map design, etc. The algorithms of mapping by GMT are improved substantially. Thus, from the geospatial perspective, GMT provides a novel approach in spatial data processing to design a series of maps with identical parameters for comparability and compatibility. The reuse of scripts applied from modular techniques of GMT and adjusted to particular datasets is presented within the framework of this project. Consequently, mapping accuracy by GMT reached high cartographic standards.
In this study, the GMT methods were applied for mapping Panama, to present a previously unpublished new series of maps on geophysical, topographic and seismic parameters of the country. This study proposes using GMT scripts to plot visualizations of Panama using high-resolution data. Thus, this paper integrates the programming methods of data processing in geographical studies to present a novel perspective of cartographic workflow. The GMT syntax was exploited to efficiently represent complex geophysical models of the terrain through automated procedure of mapping. Several geophysical variables were visualized using different datasets in a united series of maps including 2D and 3D data visualization. The data were analyzed to analyze patterns in the continuous fields and correlations between variables for data comparison and analysis of links among the geophysical and topographic parameters of Panama.

Study Area
Panama is an important region of Central America and a transatlantic transport hub with a key economic role in the world maritime trade (Fig. 1). The Panama Canal has one of the busiest traffic density in Latin America (Rodrigue, 2017;Rodrigue and Ashar, 2016). The topographic and environmental settings of the region are affected by natural factors such as high seismicity (Hopp and Waite 2016;Ambraseys and Adams 1996) and anthropogenic effects from industrial activities which results in river damming, forest flooding, species invasion, salt intrusions, chemical pollution (Salgado et al., 2020) and deforestation (Oestreicher et al., 2009). Regional development includes construction works related to hydroelectric and mineral resource exploration (Riva Palacio and Hernandez Lopez, 1998).
The Panama Isthmus is notable for the geological complexity where volcanic fronts of the eastern Central America are interrupted by a topographic low of the tectonic and magmatic origin (Buchs et al., 2019). Active tectonic faulting explains the deep links between the seismicity, volcanism and topography in Panama, which resulted in high orogeny in the area of canal. Deep associations between the geophysical and topographic setting and environmental sustainability of Panama are scarcely explained in the existing literature, except for selected studies. For instance, Dieter et al. (2010) reported on control of geology and topography on soil phosphorus in tropical rainforests of Panama, which results in biochemical changes and affects the biodiversity. Tectonic plates subduction is a major cause of earthquakes and high seismicity along the Pacific margins. Besides, seismic activity reflects the signals of the ground motions of the Earth (Lemenkova et al., 2023). Panama forms a part of the tectonic microplate resting on the Caribbean and Nazca plates that subduct under the Cocos and South America plates (Demets, 2001;Bergoeing, 2015). Such active tectonic processes result in a series of repetitive earthquakes at various depths, reported by the Incorporated Research Institutions for Seismology (IRIS) seismic center. The links between the tectonic activity, submarine volcanism and seafloor geomorphology in the Panama Basin are detected by faulted volcanic ridges created by the voluminous eruptions in the Galapagos and geomorphic processes: grabens, active deep-sea erosion and submarine slopes (Lonsdale and Fornari, 1980). The tectonic subduction of the Panama Fracture Zone is an evolutionary process that creates variations in convergence rate, obliquity, and crustal thickness at its intersection with the Middle America Trench. This process is maintained by the southeast migrating of the Panama Triple Junction (Cocos-Nazca-Caribbean). Another factor of high seismicity is the development of the inner forearc thrust belt in the colliding Cocos Ridge and Nazca-Caribbean convergence, east of the triple junction. This results in the thrust faults in the Panama Fracture Zone and contributes to the high seismicity in the region (Morell et al., 2008). Further details on the tectonic and geologic setting of Panama are described earlier (Wolters, 1986;Vergara Muñoz, 1988;Fuchs, 2012;Hardy, 1991).
Besides the unique geologic setting, Panama is an environmentally vulnerable region with a high degree of connectivity between the coastal and terrestrial ecosystems. Various approaches to spatial data analysis can be used for effective environmental mapping including geomorphological data analysis. For instance, forested mountains in western Panama (Fig. 2) provide sources of water, materials, and hydroelectric energy (Larsen, 2017). The interplay between the geochemical and geologic setting affects the sustainability of ecosystems (Gardner et al., 2014). Besides, geochemical trends are associated with multiple factors: soil, topography and land use (Neumann-Cosel et al., 2011). Therefore, linking the multi-source data through visualization by the advanced mapping methods is essential for interpreting landscape determinants and providing new environmental information.
Marine protected areas of Panama include rare species, benthic and fish communities and resources, including mangroves (Mair et al., 2009). Mangroves are one of the most productive ecosystems in the world, providing existence for coastal livelihoods (López-Angarita et al., 2018; Lemenkova and Debeir, 2023a). At the same time, natural ecosystems are linked with regional topography and experience the effects from geological setting (Benfield et al., 2005;Ruane et al., 2031;Lee and Su, 2008). For instance, the example of the environmental effects on ecosystems is described earlier (Sadatzki et al., 2019). The geochemical studies in Panama have identified the contribution of dissolved organic carbon (DOC) to the oceans from the mountainous rivers, which shows the influence of lithology and topography on the environment. The correlations between land cover and bedrock lithology show associations with sedimentation in the shelf (Goldsmith et al., 2015). The analysis of such environmental links requires spatial analysis by the advanced cartographic methods.

Data
The data include the IRIS seismological dataset with 1740 records (http://ds.iris.edu/ds/nodes/dmc/data/) collected as CSV table covering 40 years . These seismic events are located around the Panama coastline and vary according to their magnitude, depth of earthquake focus, date and time of event and coordinates. The topographic data used as a basis map for Fig. 1 and Fig. 4 are obtained from the General Bathymetric Chart of the Oceans (GEBCO) (GEBCO Compilation Group, 2020) and ETOPO1 grid (Amante and Eakins, 2009) for a 3D model in Fig. 2. The geoid data were retrieved from the Earth Gravitational Model EGM-2008 (Pavlis et al., 2012). Gravity anomalies (Figs. 5 and 6) were taken from the data available from Scripps Institution of Oceanography based on the CryoSat-2 and Jason-1 satellite observations (Sandwell et al., 2014). These data were imported into the GMT environment for processing. Apart from these layers, we used the builtin data of coastlines, riverine networks, boundaries and polygons of countries available by GMT and DCW.

Methods
This study applies a GMT (Wessel et al., 2019) to achieve the unification of maps, integration of data and implementation of GMT modules in a cartographic practice. GIS are the standard tools for cartographic visualization in a variety of applications in environmental sciences (Pelletier et al., 2019;Harper et al., 2010;Rucker et al., 2011). However, the limitations of mapping in the traditional GIS are essential, due to the lack of automation in data processing. Another problem in the commercial GIS is data format handling. In contrast, scripting methods of cartographic approach present more in-depth ways of data processing and visualization (Lemenkova and Debeir, 2022b, 2023c, 2023d. Scripting ensures an effective solution to automation in data processing. It can efficiently and rapidly manipulate with a wide range of geographic datasets and materials, including both vector and bitmap data. In this study, we present the application of GMT syntax as an effective cartographic tool. It has a significant advantage over the traditional GIS due to its scripting functionality that differs significantly from the GIS. The scripting approach of GMT builds upon command line arguments required by its modules. By changing the options of commands in the GMT, mapping can be performed semi-automatically for similar layouts. In this study, the GMT was used for 2D and 3D mapping using scripts for plotted six maps available in the authors' GitHub repository: https://github.com/ paulinelemenkova/Mapping_Panama_GMT_Scripts. The scripts were run in the console of the GMT using available methods (Lemenkova and Debeir, 2023b). The metadata of the spatial datasets were checked and the data were visualized in the Mercator projection. Apart from geographic data manipulation, we used the available raster data processing tools, including the subsetting and gridding. The cartographic projecting was selected from the available 31 projections, adjustable for the region of Panama, using the -J flag. To examine the elevation topographic component, this study used the GDAL-based option in the 'gdalinfo' utility. To visualize raster grids and to add isolines, the following modules of the GMT were used: 'gmt grdimage pa_relief.nc -Cpauline.cpt -R276.5/283/6.5/10.5 -JM6.0i -P -I+a15+ne0.75 -t30 -Xc -K > $p', and the 'gmt grdcontour pa_relief1.nc -R -J -C200 -W0.1p -O -K >> $ps'.
The isolines of the geoid curvature plotted every 0.5 m are notably curved according to the topographic elevations (Fig. 3). It is shown that the isolines of the gravity anomalies and gradient gravity grid curvatures (Figs. 5 and 6) form a set of curved lines on the maps, which well corresponds to the topographic elevations and oceanic depression in the margins of the Gulf of Panama and the continental slope of the Caribbean Sea. Subsetting data from the raw dataset is particularly useful in case of big data (millions of cells, Gb of size), where spatial extent is well known beforehand. To this end, the 'grdcut' module was used as follows: 'gmt grdcut GEBCO_2019.nc -R276.5/283/6.5/10.5 -Gpa_relief.nc'. Here the '-R' flag signifies the study extent in the westeast-south-north convention: '-R276.5/283/6.5/10.5' (0 to 360° from W to E, 0 to 90° from S to N).
Highly adjustable GMT enable performing print-quality maps and graphical illustrations through a large variety of color palette tables (CPTs), lines and polygons. The legends of color palette scales on the maps determine a possible range interval for the data extent, within which the variations of the data can be observed (Figs. 1 to 6), to define the concentration and distribution of the extreme values of the geoid undulations, geophysical grids and seismic magnitude with respect to the topography. The example of the legend (here for the map in Fig. 1) in the GMT script is as follows: 'gmt psscale -Dg276.5/6.05+w15.3c/0.4c+h+o0.0/0i +ml -R -J -Cpauline.cpt -Bg500a500f50+l"Color scale 'wiki-celtic-sea': ... [R=-3954/2426, discrete, RGB, 25 segments]" -I0.2 -By+lm -O -K >> $ps'.
The lines of the code can be broken into a series of commands with options. The most important lines of the code have the following meaning. The '-Bg500a500f50' flag controls the frequency, annotation, tick, and gridline interval for the color legend; the '-Dg276.5/6.05' defines the reference start point on the map for the color scale using the '-Dg' for map coordinates (in this case the Mercator projections); the '-I0.2' flag adds the illumination effects; the '-Cpauline.cpt' defines the selected color palette, pre-defined using the command 'gmt makecpt -Cwiki-celtic-sea.cpt > pauline.cpt'.

Spatial analysis
Modelling complex topographic surfaces is enabled by the GMT through the 3D perspective mesh-models, plots and geographic maps. This is useful to illustrate the shapes of the terrain for detailed geomorphological and bathymetric analysis. Therefore, the advanced data processing step included 3D modelling using a GMT-style syntax. The procedure of clipping the data to insert a small global map in Fig. 1 was used by the 'psclip' module (gmt psclip -R276.5/283/6.5/10.5 -JM6.0i pa.txt -O -K >> $ps) to visualize the location of Panama within the global context by the Digital Chart of the World (DCW), and to highlight the study area overlaid upon the semi-translucent image of the neighboring countries. For plotting coastal lines and borders of the countries, the 'pscoast' module was used as follows: 'gmt pscoast -R -J -Ia/thinner,blue -Na -N1/ thicker,darkred -W0.1p -Df -O -K >> $ps'. To account for complex variables (free-air gravity, geodetic data), we present an integration method of plotting data in the identical cartographic projection for comparability, and to improve the efficiency of the visualization of data with different resolution and sources.
By presenting the geographical and geophysical setting on each map, we derive geospatial information, e.g., the level of seismicity to visualize the most vulnerable areas and the variations in the gravity grids for comparison. We show through the analysis of the presented maps that the script-based method performed well for generating a series of the thematic maps covering. The 'psxy' module of the GMT was used for visualization of the earthquake events to evaluate the seismicity. The procedure was performed by reading the (x,y,z) triplets derived from the file in the IRIS .ngdc format, converted and visualized by the code that plots symbols at the coordinate locations, as follows: 'gmt psxy -R -J quakes_ PA_s.ngdc -Wfaint -i4,3,6,6s0.05 -h3 -Scc -Csteps.cpt -O -K >> $ps'. The earthquakes were plotted using the '-Sc' flag that visualizes circles where the size of the diameter signifies the earthquake's magnitude. The magnitude and spatial variability of the earthquakes were indicated in the '-i4,3,6,6s0.05' flag that signifies the numbers of columns in the table with the following names: Year, Month, Day, Lat, Lon, Depth, Mag.

Results
The results include previously unpublished new maps visualizing the complex geophysical parameters of Panama through 2D and 3D modelling. The comparative analysis of maps shows the links between the topographic and geophysical setting, e.g., the correlations between the geoid, elevations, gravity and seismicity. The color distribution palettes were adjusted to each of the maps and implemented using the CPTs; thus, the colors were produced and selected to best indicate the variability in the parameters. The topographic map and the geoid model were compared to reveal the exceptionally high values of geoid (Fig. 3) notable in the high mountainous regions of Panama (17-20 m around Cordilliera de Talamanca, southern coast of Peninsula de Azuero and eastern Panama, 77.5-78.5°W). This corresponds to the topographic roughness in the highlands.
In contrast, negative values of geoid can be seen over the Caribbean Sea basin in the north (from −4 to 0 m). High resolution mapping of the seismicity (Fig. 4) demonstrated 1740 earthquake events in Panama at various depths with different magnitudes and spatial location. A notable concentration of the earthquakes is seen in the western region of the Panama's shelf waters (~82-83.5°W), and on the border with Colombia (~77-78.5°W). The earthquakes vary by magnitudes from 2.9 to 7.8 at the depths from 0 (shallow earthquakes) to 225 m (near the west coast of Colombia), that are traceable from the IRIS seismicity data (Fig. 4).
The hallmark of the high values in free-air gravity anomalies (Figs. 5 and 6) and vertical gradient is a concentration of the land masses in the seismicity with a detected correlation in the zone of the tectonic subduction and earthquakes. This supports the mountain areas of Panama. The higher geophysical anomalies can be noted by the comparison of Fig. 5 with Fig. 6, and Fig. 1 with Fig. 2. The lowest values correlate with the depression in the seafloor along the western coasts of Panama, which continues the northern extent of the Middle America Trench. The topographic and geophysical patterns are compared to previous studies regarding the associations between the seismicity and topographic roughness of the seafloor (Bilek et al., 2003). The appearance of the high values in gravity anomalies (over 220 mGal) in the mountainous regions matches the geodynamic processes associated with the Earth structure, inner interior, and crustal surface, as well as the geodetic and geophysical effects. This considers the shape and rotation of the Earth, and land mass distribution, which results in the contrasting topographic patterns and distribution of the topographic elevations: depressions, deep oceanic trenches and mountain ranges. The results prove the reports on seismic activity and tectonic instability in the Panama region. The presented maps of seismically active areas in Panama correspond to prior studies that indicated that Costa Rica Ridge and the Panama Fracture Zone are significant tectonic features related to the thermal anomalies in Panama (Vargas et al., 2018).
The recordings of earthquakes on the oceanic bottom seismographs showed magmatic and hydrothermal activity of the upper lithosphere in Panama region that are a source of seismic hazards and risks. Using the comparison of gravity, geoid and earthquake data, the relationship between the seismic activity and topography of Panama was explored. The presented series of new maps are applicable as indicator and inventory data sources to analyze hazard vulnerability and perform risk assessment in the regions prone to seismic activity along the coasts of the Caribbean Sea and the Pacific Ocean.
Earthquakes, their spatial distribution, magnitude and frequency are an important source indicating hazard risks in the region. At the same time, earthquakes are influenced by regional geologic and tectonic factors. These include, among others, lithosphere plate subduction, run-off, topography, distribution of geologic fault lines and regional geologic context, soil type (clayey, sandy). In this context, highlighting the correlations between these parameters allows taking a deeper insight into the factors controlling the seismicity in Central America. Thus, the distribution of the continuous fields on gravity is correlated with topography. At the same time, gravity is poorly correlated with seismicity and active tectonic zones in western Panama, 81-82.5°W), although there is a certain impact of the geodynamic phenomena on geological processes. The results indicate that the most notable series of earthquakes seismic activity are concentrated in the following regions of Panama: • south-western region (square of 82-83°W, 7-10°N), owing to the tectonic faulting; • south-eastern region in the border with Colombia (square 77-78.8°W, 6-8.5°N) due to the tectonically active margins of the subducting Nazca Plate along the South American continent; • sporadic, less prominent earthquake activity in the southern coasts of the Caribbean Sea (78-80°W, 9.3-10°N), which notes the movements of the Caribbean plate colliding with the South American and Cocos plates in a more complex tectonic movement.
The tectonics of the Caribbean Sea and the marginal eastern Pacific is known to largely influence regional geologic variability and seismicity. Another factor is the conjunction with a greater topographic heterogeneity in the terrestrial areas and seafloor (Rodríguez-Zurrunero et al., 2020;Idárraga-García et al., 2021;Naranjo-Vesga et al., 2020;Duncan et al., 1999). This interplay between the geophysics, geology and tectonics explains the correlation between the seismicity and topography with geology and geophysics (Lemenkova and Debeir, 2022a). At the same time, when this correlation is analyzed within regional zones of Panama (SW, northern coasts of the Caribbean Sea, SE at the Colombia affected by the Andean geology), stronger correlations that are observed between the magnitude, distribution and frequency of earthquakes are notable.

Discussion
We presented a previously unpublished new series of maps showing geophysical setting of Panama made using GMT scripts. To the best of our knowledge, no existing maps on Panama using the GMT scripts exist in the available literature. The functionality of GMT syntax as an advanced cartographic tool showed effective results in terms of mapping quality and cartographic workflow. The maps illustrate a high level of seismicity in the Gulf of Panama, caused by active tectonic movements, which result from the colliding lithosphere plates. High seismicity in this region presents risks for the economic functionality of the Gulf of Panama, which is a strategic maritime port connecting the Atlantic and Pacific oceans, and for the environmental sustainability. The effectiveness of GMT scripts for geophysical and topographic mapping of Panama was showed both in terms of illustrative quality of the cartographic materials and technical soundness of the programming approach.
The Gulf of Panama is a strategic maritime hub of high logistic importance for transport, which connects the Atlantic and the Pacific oceans. Therefore, mapping seismicity by advanced cartographic methods in such a region is an important contribution to data analysis for mitigating possible risks and predicting earthquakes. Earthquakes are recognized as hazardous geologic events accompanied by devastating processes: landslides, tsunami, rock falls, damage of buildings and losses of human lives. Mapping earthquakes provides an important source of information for risk assessment in the regions of high seismicity. However, such studies for Panama are limited, compared with other Central American countries with a long tradition of geological studies. As a response to these needs, this paper presented a series of thematic geophysical and seismic maps of Panama plotted using the practical workflow of the advanced methods of GMT programming scripts.
Our main contribution to the cartographic studies on Latin America consists in a novel series of maps covering Panama, and the analysis of the variations in topographic and geophysical parameters. Our study proves the associations between the geological-geophysical parameters and environmental processes in Panama.
The geophysical monitoring of this region is maintained by the high quality maps aimed to indicate seismic risks. Accurate visualization of the seismically active areas serves to the prognosis of earthquake hazards.

Conclusions
Mapping the marine regions benefits from cartographic scripts due to the increased speed of plotting and automated workflow that ensures high quality of maps. This paper presented the GMT-based mapping of the region of Panama based on the high-resolution data including IRIS seismographic database . The presented methodology of this study included the examples of the explained and commented scripts made using console-based techniques of GMT. The link to the GitHub repository of the authors with open scripts is provided for repeatability of scripts. The specific focus of this article is on seismic and geophysical mapping made using the datasets on gravity and earthquakes. The presented maps enable the analysis of correlation between the earthquake distribution, seismicity, topographic and geophysical setting in Panama.
This study demonstrated the integrated analysis of the seismicity in the Gulf of Panama in the context of its topographic, geophysical and environmental setting. The importance of the Gulf of Panama as a study area is explained as the combination of the following factors: (i) a unique location connecting the two oceans, the Atlantic Ocean (Caribbean Sea) and the Pacific Ocean; (ii) the bathymetry in both oceans differs with clearly visible depressions in the seafloor; (iii) tectonic activity of the subducting plates, which largely affects the seismicity; (iv) mountainous geomorphology in the central regions of Panama.
The GMT scripting cartographic methods employed in this study offer the improvement over the GIS-based geological and geophysical mapping of Panama in several technical and methodological aspects: 1 The data-driven mapping requires rapid and automated data processing that can only be performed using the advanced methods such as GMT, the advanced cartographic toolset for spatial data processing, modelling and visualization.
2 This study relies upon the most reputable seismic database from IRIS, which is updated regularly using the data from multiple seismic stations worldwide. IRIS uses generated reports with a standardized model of the data capture and comprehensive metadata: event earthquake frequency, depth, magnitude, precise coordinates and description of locations. These data were downscaled to the regional level of Panama.