Our work relies on three datasets: (1) global digital elevation models and bathymetric data for both Earth and Mars; (2) maps of fluvial features on Earth (major global rivers and deltas) and Mars (valley networks, fluvial ridge systems, outlet canyons and interpreted deltas), along with maps of oceanic features on Earth (continental shelf, shelf break and ocean floor) and interpreted shorelines on Mars; and (3) results of elevation, slope, curvature and landscape classifications for each cell on both Earth and Mars. We describe the data and outline our methods for each dataset below and briefly explain the flowchart in Supplementary Fig. 1.
Global digital elevation data and bathymetric data
Earth
Three global digital elevation models integrate both land and ocean surfaces at different resolutions and levels of consistency: (1) the ETOPO Global Relief Model with a general average resolution of about 1.85 km per pixel63; (2) the SRTM30_PLUS Estimated Topography with a resolution of around 1 km per pixel64; and (3) the General Bathymetric Chart of the Oceans (GEBCO) with a resolution of approximately 500  m (ref. 65). For our global-scale topographic analysis, we used the ETOPO1 Global Relief Model because of its consistent pixel resolution across both terrestrial and oceanic regions. ETOPO1 provides a uniform 1 arcmin resolution (about 1.85 km per pixel), integrating satellite altimetry, shipboard echo-sounding and terrestrial measurements into a cohesive dataset63. This consistency is important for studies requiring seamless data across different terrains, ensuring that both land and ocean topography are represented with the same level of detail. By contrast, the SRTM30_PLUS Global Bathymetry and Topography dataset, while offering higher resolution for land areas (30 arcsec), lacks uniformity as it focuses primarily on terrestrial regions and provides less detailed coverage for the oceans64. Moreover, the GEBCO dataset, although detailed for ocean bathymetry, does not offer the same consistent pixel resolution for land topography, leading to potential discrepancies when integrating land and ocean data65. Therefore, ETOPO1 was selected to ensure uniform resolution and comprehensive coverage across both terrestrial and oceanic environments, addressing the need for consistent pixel data in our analysis.
Mars
We used the global Mars Orbiter Laser Altimeter (MOLA) gridded topography, which offers a pixel resolution of 463 m per pixel66. This dataset is derived from more than 600 million measurements covering the entire Martian surface. These measurements were meticulously adjusted to ensure consistency, providing a uniform pixel resolution across the entire terrain of Mars67,68.
Data resampling
We resampled both digital elevation models to multiple resolutions—2.5 km, 5 km and 10 km—for several key reasons: (1) Resampling the topographic data of Earth and Mars to a uniform resolution is essential to apply uniform analytical methods and enable direct comparison of topographic features across both planetary surfaces. (2) These specific resolutions were selected to intentionally exclude fine-scale landforms on Mars, as these are generally younger in age69. Our study, however, focuses on older, broader-scale topographic features that provide insights into ancient surface processes. The resampling was conducted using the ‘Resample’ tool in ArcGIS70,71, in which we applied the ‘nearest neighbour’ option to preserve the exact elevation values, minimizing significant interpolation and smoothing, and thereby ensuring the integrity of the original data72,73.
Maps of fluvial and oceanic features on both Earth and Mars
To identify a rough search zone on Mars for the transition from landscape to seascape, we used maps of the major world rivers and deltas on Earth40,41. These typically indicate where the terrestrial landscape ends and the oceanic zone begins. However, this assumption holds only if the rivers and deltas were active simultaneously. Changes in sea levels could alter this relationship, but the maps still provide a useful approximation of the extent of the transition zone. Moreover, we used an extensive dataset mapping seafloor geomorphic feature42, which not only helps to define the zone but also offers insights into how oceanic geomorphic features evolve spatially. We focused on the key geomorphic features that define the transition: the continental shelf, the shelf-break slope, the continental slope, the continental rise and the key ocean floor landforms (abyssal and hadal zones).
We relied on global mapping of valley networks44,74,75, outlet canyons44, depositional rivers (fluvial ridges)45 and interpreted deltas4,10,14,15,16,17,76. Furthermore, we used maps of topographic contacts, previously interpreted as shorelines10,22,77, to assess how the water-formed landscape functioned.
Given the debate over whether Martian deltas formed in open or closed basin systems, we chose to compile the available delta datasets and then apply specific filtering criteria4,10,14,15,16,17,76,78,79. We selected deltas that (1) are open to downstream flow and located along the dichotomy boundary and/or (2) exhibit complex stacking patterns interpreted as evidence of formation within either regressive or transgressive depositional environments. This filtering resulted in a set of 48 deltas (Supplementary Table 1 and Supplementary Fig. 4). We further examined these deltas and classified them into two categories based on their morphology: single-lobate deltas and stacked deltaic systems. To further cross-validate our compilation, we calculated the elevations of all channels and lobes—including those preserved as ridges and interpreted as erosional remnants of ancient deltas—within the largest deltaic systems in Aeolis Dorsa and Hypanis14,15,16,17, to capture elevation changes potentially associated with past sea-level fluctuations (Fig. 5e,f and Extended Data Fig. 7b,c).
Setting a search zone for landscape-to-seascape transitions
We converted the shapefiles of terrestrial rivers, deltas and Martian rivers, deltas and proposed shorelines into points in ArcGIS Pro70. We then calculated the elevation of each point and plotted the results to examine where the continental zone transitions to the oceanic zone. For oceanic geomorphic features (polygons), we used the zonal statistics tool in ArcGIS Pro to calculate the number of pixels within each polygon, total area and the 10th–90th percentile elevation values for each zone.
We extracted elevation data for the major global terrestrial rivers (195,022 data points), global deltas (10,848 data points), the continental shelf (14,820,634 data points), the continental slope (7,606,463 data points) and the key ocean floor landforms, including the continental rise (12,144,045 data points), abyssal plain (116,749,407 data points) and hadal zone (1,238,491 data points) on Earth to identify the upper and lower bounds of the continental shelf. On Mars, the analysis included valley networks (3,294,322 data points), depositional rivers (16,515 data points), outlet canyons (248,865 data points), deltas (48 data points), and the Arabia (10,192 data points) and Deuteronilus (42,900 data points) shorelines, to establish the upper bound of the potential Martian shelf. On Earth, river deltas prograde across and rest atop continental shelves, and the transition from these deltaic deposits to the deep ocean typically takes place within the upper 2.5 km below sea level. We, therefore, use this depth interval to define the search window for a potential shelf on Mars.
Raster to points of elevation, slope and curvature
To obtain elevation at each point of the resampled raster files, we used the ArcGIS ‘Add Surface Information’ tool to sample elevation values from grids at resolutions of 2.5 km, 5 km and 10 km (refs. 70,71; Supplementary Figs. 2 and 3). For each point, the z-value is derived from its x–y coordinates on the underlying surface.
To calculate the steepness of each cell on both terrestrial and Martian surfaces, represented as raster grids, we used the ‘Slope’ tool in ArcMap71. Slope (degrees) was calculated in ArcMap (Spatial Analyst) using the Slope tool (3 × 3 neighbourhood), which estimates ∂z/∂x and ∂z/∂y using a finite-difference gradient and computes \({S}^{^\circ }=\arctan (\sqrt{({\partial z/\partial x)}^{2}+{(\partial z/\partial y)}^{2})})\times 57.29578\).
Curvature shows the shape of the slope, indicating whether it is convex (that is, ridges and plateaus) or concave-up surfaces (that is, valleys and depressions), which is particularly useful for identifying transitions between landscape and seascape. We calculated curvature using the ‘Curvature’ function in ArcMap71, which fits a plane to the nine surrounding cells in a 3 × 3 window to determine surface curvature. The primary output provides cell-by-cell curvature values (second derivative of elevation). The positive values indicate upwardly convex surfaces, negative values indicate upwardly concave surfaces, and values near zero represent flat or nearly planar areas. In ArcGIS, curvature values are reported in units of one hundredth (1/100) of the DEM z-unit (here, metres; z-factor = 1). To treat concave and convex areas equally in terms of magnitude, negative values of both planets are multiplied by −1, allowing them to be plotted alongside positive values with different colours (Fig. 2).
Landscape classification
To map the transition zone between continental and oceanic landscapes on both Earth and Mars, we used the ‘Geomorphons’ tool53. This algorithm uses the concept of Local Ternary Patterns (LTP) to analyse terrain features by comparing the elevation of each pixel with its neighbouring pixels54. Instead of a simple binary comparison, LTP classifies differences into three categories: values that are (1) similar to the centre pixel, (2) significantly higher or (3) significantly lower. This approach reduces the impact of noise and provides a more nuanced representation of terrain, capturing subtle variations in pixel elevation. The Geomorphons tool classifies each cell of an input raster into common landforms, including flat areas, ridges, shoulders, spurs, slopes, pits, footslopes, hollows and peaks.
On Earth, the transition typically occurs on a relatively flat surface, known as the continental shelf29. To detect these flat surfaces, we used the Geomorphons tool with a specific ‘flat terrain angle threshold’. The tectonic system of Earth, driven by active plate tectonics, differs significantly from Mars, which lacks substantial tectonic activity. This absence of tectonism on Mars results in longer topographic wavelengths compared with Earth52. As a result, we applied different flat-terrain angle thresholds for the two planets.
For Earth, we conducted 40 Geomorphons experiments (Supplementary Table 2) to (1) classify the surface of Earth into shelf cells and non-shelf cells by comparing the detected flat cells to the mapped continental shelf; (2) determine the flat angle threshold that fully detects the terrestrial continental shelf; and (3) establish a range of shelf area detection at each angle threshold (Extended Data Figs. 2–4 and Supplementary Table 2), which will be used to assess the percentage and accuracy of shelf detection on Mars. In each experiment, we applied a different flat-terrain angle threshold and found that the continental shelf was fully detected at an angle of 1.22° (Supplementary Table 2). However, increasing the flat angle threshold led to the detection of both shelf and non-shelf areas. For instance, we detected 100% of the shelf (32,308,476 km2), but it also identified 238,719,226 km2 of non-shelf areas, resulting in a precision of 12%. We, therefore, used the experiments to set an accuracy range for each flat angle threshold, which would be used to detect the Martian shelf (Extended Data Fig. 3 and Supplementary Table 2).
For Mars, no definitive maps of oceanic features exist, apart from two long-debated proposed shorelines10,22,77. To map potential shelf-oceanic zones, we focused on a region in the northern lowlands (Fig. 2c) characterized by a distinct flat zone compared with its surroundings. This region also preserves 48 deltaic systems, some of which are connected to interpreted submarine-channel belts that are thought to have formed along an ancient oceanic margin4,10,14,15,16,17 (Supplementary Table 1). Moreover, it contains numerous valley network termini and fluvial depositional ridges, both of which are believed to represent the endpoints of fluvial systems. The region also preserves the two debated shorelines. This zone is well-defined by elevation, ranging from −1,800 m to −3,800 m, with a distinct median slope of 0.31° at a grid resolution of 5 km. We ran the tool at this flat angle threshold and found that nearly the entire northern lowland was marked as a relatively flat surface. To refine the results, we combined our topographic analysis (elevation, slope and curvature) with key morphological features (valley network termini, fluvial ridges, deltas and the two proposed shorelines). This allowed us to spatially constrain the flat surface zone between −1,800 m and −3,800 m to areas coinciding with geomorphic indicators of a landscape-to-seascape transition, mostly located between about 30 °S and about 70 °N. On Earth, a flat angle threshold of 0.31° would detect nearly 69–71% of the continental shelf area, giving us confidence that the detected surface on Mars corresponds to a similar transition (Supplementary Table 2).
Statistical analysis
To test whether surface steepness differs between elevation zones, we first computed median slope and curvature values at 200 m elevation intervals on both Earth and Mars. The resulting Martian profiles show an intermediate-elevation, low-slope, low-curvature interval between −1,800 m and −3,800 m, bounded by higher slopes and curvature at elevations >−1,800 m and lower slopes and curvature at elevations <−3,800 m (Fig. 2 and Extended Data Fig. 1). On this basis, we defined three elevation bands: >−1,800 m, −1,800 m to −3,800 m, and <−3,800 m. We then applied a Kruskal–Wallis H test80 to these three elevation bands to quantify whether their slope distributions differ, without assuming normality. The test showed a highly significant difference in median slopes (H = 27.50, P = 1.07 × 10−6; Extended Data Fig. 1), indicating that the slope populations of the three elevation zones are statistically distinct. The Kruskal–Wallis test is used here to assess only the distinctness of elevation zones defined from the slope–curvature–elevation relationship, not to locate the breaks themselves.
Limitations
Our results are subject to several limitations. One is the potential alteration of topography due to true polar wander and the emplacement of the Tharsis volcanic province, which probably caused uplift near Tharsis and subsidence farther away11,12. However, analysing the surface as geomorphic domains helps mitigate this, because this deformation would affect broad regions rather than the specific elevation ranges of landform mosaics. A second source of uncertainty is the isostatic response to ocean unloading, which on Earth can modify elevations by several hundred metres following ocean retreat81. However, recent estimates for Mars suggest that isostatic rebound probably ranged from several tens to just more than 100 m (ref. 12). Yet, the approximately 2 km elevation span of our detected shelf-like zone exceeds expected rebound estimates and remains consistent with depositional features. A third source of uncertainty is long-term burial, exhumation and erosion46. Although these processes may have introduced regional variability, they are unlikely to alter the broader topographic patterns we identify at the global scale. A fourth source of limitation is the erosion and sediment redistribution along the dichotomy boundary by Hesperian-aged outflow floods13,48, which probably deposited substantial volumes of sediment along the northern dichotomy, particularly in Chryse Planitia. These outflow events probably contributed to locally flattening the surface there. However, similarly flat, low-slope surfaces are also present at other key sites, such as Aeolis Dorsa—rich in stacked deposits of different origins and interpretations, including fluvial and deltaic deposits, and possibly even submarine deposits15,16,17,82,83,84,85,86—and along the remaining segments of the proposed shelf. This broader distribution, together with independent evidence for sea-level changes recorded by deltaic deposits at Hypanis14, suggests that although Hesperian-aged outflow floods helped flatten the surface in Chryse Planitia, they were probably not the primary cause of surface flattening across the northern lowlands.
