The hypothesis of the existence of a Martian ocean in the northern lowlands remains an intriguing open question about the early stages of the planet’s evolution^1,2,3,4. The existence of the ocean may have significantly influenced the climate and atmosphere of early Mars (e.g.,⁵) and left geological records of its presence^2,3.

China’s Mars rover Zhurong, onboard the Tianwen-1 probe, successfully landed in southern Utopia Planitia on Mars in May 2021^6,7. This region has long been hypothesized to be a part of an ancient ocean that once covered the northern lowlands^2,8,9,10. The Zhurong landing site represents the tenth in situ investigation on the Martian surface and the first in southern Utopia Planitia. Using its onboard payloads, such as the Navigation and Terrain Camera (NaTeCam), Mars Surface Composition Detector (MarSCoDe), Multispectral Camera (MSCam), and Mars Rover Penetrating Radar (RoPeR)¹¹, the Zhurong rover has investigated the topography, regolith structure and geology of the roving area, with the scientific objectives of analyzing the elements, minerals, and rock types to search for signatures of water or ice in the landing area¹².

Figure 1 shows the Zhurong landing site (109.925° E, 25.066° N) in southern Utopia Planitia on a geologic map¹³. The high plateau, Nepenthes Planum, which defines the boundary between the plains in the north and the highlands in the southern hemisphere, is approximately 400 km southward of the landing site. The volcano Elysium Mons is approximately 2000 km to its east. Zhurong landed in a Late Hesperian lowland (lHl) unit¹³. Southern to the lHl unit is a Late Hesperian transition (lHt) unit. An Amazonian and Hesperian volcanic (AHv) unit and an Amazonian volcanic (Av) unit are in the east, close to Elysium Mons.

The Zhurong landing site (red cross) shown on the geologic map by Tanaka et al.¹³. Regions of different colors indicate different geologic units. The red box shows the study area, covering an area of ~ 1800 km × 800 km including the lHl, lHt, Av, and AHv units. The thick black line is the Deuteronilus contact¹⁷ indicating a possible shoreline across the region. The background map is a shaded relief of the MOLA topography. The software ArcGIS 10.1 (https://www.esri.com/) was used to create the map in the figure.

Full size image

The lHl is a Hesperian-aged sedimentary geologic unit mapped from orbital data that extends across most of the Martian northern plains¹³. It is mainly composed of Vastitas Borealis Formation (VBF) materials^9,14,15. The VBF deposits might have originated from fluvial, lacustrine, or marine sediments^16,17,18, could be the result of sedimentation from highlands debris flows^19,20, could have formed due to groundwater discharge/mud volcanism^4,21,22,23, or might have been formed during volcanic activity²⁴. The volcanic origin of the deposits is dismissed due to the material’s physical properties, such as fine grain size¹⁴ and dielectric constant value²⁴. Thus, in our paper, we focus on the vital role of water/mud in deposit formation in the Utopia region.

Previous studies of the northern lowlands on Mars (e.g.,²⁵) have revealed curvilinear geomorphic features (such as the Deuteronilus contact shown in Fig. 1) in the Pyramus-Astapus, Utopia, and Western Elysium regions that extend for over 7700 km^1,2,17. The Deuteronilus contact has a global altitude variation of about 500 m, and in southern Utopia Planitia, it has a mean elevation of approximately − 3580 ± 270 m¹⁷. The Deuteronilus contact delineates the border of a depositional unit that overlaps the underlying surface. It has been interpreted as the shoreline of an ancient ocean formed in Late Noachian to Early Hesperian^3,9,17,26,27. In southern Utopia Planitia, the Deuteronilus contact closely coincides with the boundary between the lHl and lHt units (Fig. 1). The mean elevation of the Deuteronilus contact is close to the mean southern boundary of the VBF deposits, at −3658 ± 282 m²⁸, and is comparable with the elevation of outflow deltas in Chryse Planitia, which is −3724 ± 153 m²⁹. The channels in Chryse Planitia carried a significant amount of water from the highlands and might have filled the northern lowlands in the Hesperian period^2,3,5,27. The brightness of the VBF unit on Thermal Emission Imaging System (THEMIS) infrared (IR) images is higher in the daytime images and lower in nighttime images¹⁷. This variation in brightness corresponds to a significant presence of fine-grained materials in these deposits²³ and is consistent with the sedimentary nature of the unit. A number of coastal landforms^30,31,32 along the possible shoreline^2,3 and the depositional behavior of the VBF unit^14,17 support a sedimentation mechanism of origin involving heavily silted water. This may have occurred as debris flow sedimentation^19,20 and/or from a single standing water/mud reservoir^2,3,5,27. However, it should be noted that VBF deposits have a homogeneous distribution and lack well-defined flows²⁰, which are typical of debris flows³³. Apparently, the outflow channels of the Chryse Planitia region are a main source of water/mud that spread around the northern lowlands, forming an extensive single-standing reservoir² or isolated bodies that may have formed independently^4,19. As will be discussed below, widely spread water/ice-related landforms in southern Utopia indicate a significant amount of water in this region. The combination of these factors allows us to assume the presence of an extensive reservoir of water/mud in this region.

In this study, we present a comprehensive analysis of the Zhurong landing area using remote sensing data and in situ observations. For the systematic investigation using remote sensing data, we have selected an area of approximately 1800 km × 800 km covering the Zhurong landing site and its neighborhood of different geologic units¹³, as shown by the red box in Fig. 1. We particularly examined the lHl unit in the study area and found that it was not uniform in terms of distributions of geomorphological features such as pitted-cones^34,35, polygonal troughs^36,37, etched flows¹⁴, and impact craters as shown in Figs. 2, S1, and S2. Pitted-cones are small cone-like structures with summit pits and are common on the vast plains in the northern lowlands of Mars³⁸. Different mechanisms have been suggested for the formation of pitted cones, including mud volcanism^22,39, rootless cones⁴⁰, cinder/scoria cones⁴¹, tuff cones/rings⁴², and pingos⁴³. Some pitted cones in the Hephaestus Fossae region, located near Elysium Mons in the eastern part of the study area, exhibit morphometric and spectral similarities to rootless cones⁴⁴. However, most studied pitted-cones in the southern Utopia Planitia suggest that mud volcanism is the primary mechanism of formation based on morphometric parameters, which may be related to explosive eruptions driven by volatiles^45,46. Some previous studies have associated mud volcanic pitted cones in Isidis Planitia with a shallow marine environment and/or the formation of ice sheets^47,48,49,50. Polygonal troughs are curvilinear extensional features cutting the surface into polygonal terrains. Some troughs have concentric shapes related to covered ghost craters^51,52 or stealth craters⁵³. However, in the study area, concentric troughs likely formed over the covered ghost craters^17,54,55. One of the possible interpretations of the formation of the troughs in Utopia Planitia is the surface subsidence after the elimination of the ice sheet load and/or fluid expulsion through the cracks in the ice shield^36,51,56. Etched flows are spatially associated with troughs. The edges of etched flows are disintegrated into numerous small mesas and chaotically oriented short ridges. The absolute model age (~ 3.6–3.5 Ga)¹⁴ and structure of the material suggest the origin of etched flows as the results of eruptions of liquefied (mud-like) materials from either partly or wholly liquefied sources^14,57. Troughs, and the associated etched flows, were presumably formed in an environment with a significant amount of water/ice (e.g.³⁶).

Morphologically different terrains and the new division of geological units in the study area. (a) Geological units defined in our study and their absolute model ages (see Fig. 3). Unit boundaries of Av, AHv, lHt, and lHl are from Tanaka et al.¹³, with lHl unit further divided into HNl_sm and HNl_dm and lHt renamed as lNt based on AMA analysis and distribution of geomorphological features; (b) map of pitted-cone density, polygonal troughs, and etched flows; (c) MOLA elevation map⁵⁸, (d) map of crater density, pancake-like and rampart craters; (e) THEMIS-IR mosaic of night-time images (100 m/pixel)⁶¹. The software ArcGIS 10.1 (https://www.esri.com/) was used to create the maps in the figure.

Full size image

The remote sensing observations were supplemented with in situ analysis discussed in previous research and continued in this study to unveil features consistent with the existence of a nearshore zone in southern Utopia. Based on these findings, this study further suggests an evolution scenario for the nearshore zone in southern Utopia in the context of the Martian ocean theory. Details are presented in the following sections.

Remote sensing and in situ data used for analysis

We used the following remote sensing data to determine the morphological features and provide a geological analysis of the surface. Elevation data were derived from the Mars Orbiter Laser Altimeter (MOLA)⁵⁸ with a resolution of ~ 460 m/pixel. The Context Camera (CTX) image mosaic (~ 5 m/pixel)⁵⁹ and High Resolution Imaging Science Experiment (HiRISE, ~ 25–30 cm/pixel)⁶⁰ images were used for detailed geomorphologic and geologic analysis. The Thermal Emission Imaging System (THEMIS) data with a resolution of 100 m/pixel were used for the analysis of thermal properties^61,62.

For the analysis of sedimentary rocks we used the data collected by the Navigation and Terrain Camera (NaTeCam⁶³), Mars Surface Composition Detector (MarSCoDe⁶⁴), and the Multispectral Camera (MSCam⁶⁵) on board the Zhurong rover. The NaTeCam has two identical cameras separated by a rigid base with a length of 270 mm. Each camera has a pixel resolution of 0.396 mrad and a focal length of 13.17 mm. The NaTeCam cameras are capable of acquiring stereo panoramic images around the rover, which allows the generation of digital elevation models (DEMs) of favorable resolutions (3 mm/pixel at a distance of ~ 5 m from the rover) for geometric measurement⁷. The MarSCoDe includes a Laser-Induced Breakdown Spectrometer (LIBS) (240–850 nm wavelengths range) and a short-wave infrared (SWIR) spectrometer (850–2400 nm)¹¹. The LIBS can remove the dust or coating on rocks to determine the original rock composition. It can then allow the SWIR imager for detailed analysis of mineralogy and classification of the rocks. The MSCam has eight bands with a spectral range from 460 to 1050 nm. It has a higher pixel resolution of 0.146 mrad compared with the NaTeCam.

Crater size-frequency distribution (CSFD) analysis for surface age estimation

CSFD analysis was conducted to estimate the absolute model ages (AMA) of each geological unit. The CTX image mosaic (~ 5 m/pixel) covering the study area was used for the automatic detection of craters (≥ 200 m) based on an active machine-learning approach⁶⁶. The approach has a detection rate of approximately 85%. Extensive manual checking was further conducted to increase the detection rate. To minimize the influence of secondary craters on the reliability of the age estimation, only large craters (≥ 800 m in diameter) were used, and possible secondary craters were manually identified and excluded based on features such as the radial chain or clustering distributions and visual distinctness from primary craters⁶⁷. For example, the crater density map (Fig. 2d) shows several regions of crater clusters and radial patterns, likely representing secondary craters of distant primary craters. After removing possible secondary craters, the AMA analysis using CraterStats-II software⁶⁸ by applying the production function of Ivanov⁶⁹ and the chronology function of Hartmann and Neukum⁷⁰ were performed to estimate the ages of the units.

Spectral analysis

The MarSCoDe includes a Laser-Induced Breakdown Spectrometer (LIBS) and a short-wave infrared (SWIR) spectrometer. The SWIR spectrometer can identify the mineralogy of the surface with a standoff distance of 1.6–7 m. The SWIR data have 321 bands with 850–2400 nm wavelengths, 3–12 nm spectral resolution, and a field-of-view of 36.5 mrads. The SWIR spectra were acquired from the Lunar and Planetary Data Release System (https://moon.bao.ac.cn/) and processed into level 2B⁷¹ through dark current calibration and absolute radiometric calibration. The accuracy of the absolute radiometric calibration is better than 5%⁶⁴.

The absolute reflectance of the target was calculated and corrected by a standard white calibration board with an average reflectance of 0.99 across the full wavelength range⁷².

where Ref_tar and Ref_cal are the reflectance of the science target and calibration board, respectively. Rad_tar and Rad_cal are the radiance data of the science target and calibration board.

The SWIR reflectance curves of the targets were compared with laboratory spectra, including hydrated silica (opal TM8896, chalcedony CU91-6A), sulfate (kieserite F1CC1S, gypsum BKRUB556), phyllosilicate (allophane CIJBA) obtained from the USGS spectral library⁷³.

A dichotomy of different terrains at Zhurong landing area from remote sensing

In the study area, we have mapped 11,895 pitted-cones (with sizes from 34 m to 8.6 km in diameter), 2773 polygonal troughs, and 77 etched flows based on a CTX image mosaic (~ 5 m/pixel)⁵⁹ covering the study area. Our mapping results are comparable with previous investigations but are conducted at a larger scale compared to those by Ivanov et al.¹⁴ and Mills et al.⁷⁴, or cover a larger area compared to Zhao et al.⁷⁵ and Ye et al.⁴⁶. A density map was generated for the pitted-cones by counting their numbers within a moving circular window of 100 km². The polygonal troughs and etched flows are shown on the pitted-cone density map in Fig. 2b. Examples of the pitted-cones, polygonal troughs, and etched flows can be seen in Fig. S1. Figure 2b shows that the polygonal troughs and etched flows are spatially associated with each other and distributed mainly in the northern part of the lHl unit (lower elevation) in the study area, which corresponds to previous mapping (e.g.,^14,15); on the contrary, the pitted-cones are mainly distributed in the southern part of the lHl unit (higher elevation). Therefore, the study area could be subdivided into northern and southern sub-units. This variation of spatial distributions of the pitted-cones and the polygonal troughs (and etched flows) at different elevation levels indicates that they could have formed under different environmental regimes.

We also examined the distribution of impact craters in the study area. Over 460,000 craters with diameters ranging from 200 m to 32.8 km were mapped based on the CTX image mosaic using an active machine learning algorithm⁶⁶ followed by extensive manual checking, from which a crater density map was generated by counting the number of craters in a moving circular window of 100 km², as shown in Fig. 2d. The northern part of the lHl unit has lower crater densities as compared with the southern part. The higher-density of craters in the southern portion of the unit, as shown by the red clusters in Fig. 2d, is due to the concentration of non-randomly distributed secondary craters. We also mapped numerous craters with a pancake-like or rampart ejecta in the study area as shown in Fig. 2d. Both types of ejecta have been interpreted to indicate the presence of volatiles (water, ice, or both) in the target materials (e.g.,^76,77), with the pancake-like ejecta indicative of more abundant volatiles and rampart ejecta suggesting a smaller amount of volatiles in target materials^14,78. Figure 2d shows that the rampart craters (black lines) are distributed uniformly; however, the craters with pancake-like ejecta (purple lines) are mostly concentrated in the northern part of the area and their numbers are increasing northwards, indicating the presence of more volatile-rich materials there. Figure S2 shows the distribution of mapped craters and detailed examples of craters with pancake-like or rampart ejecta.

Based on the above observations, the lHl unit was divided into two sub-units: HNl_dm (northern part) and HNl_sm (southern part). The border between these two sub-units was created by extracting the centerline between the feature density maps of the northern cluster of features (polygonal troughs and etched flows) and the southern cluster of pitted-cones (Fig. S3). The centerline was further adjusted in local regions according to the local variations in thermal properties as retrieved from the THEMIS IR day and night images⁶². Based on mean distance calculations, the division between these two sub-units generally follows the −4200 m contour (Fig. 2c), however, it deviates from the −4200 m contour in some local regions due to differences in the THEMIS-IR image (see enlarged views in Fig. S3). A notable difference between the centerline of feature density maps and the adjusted line in the southeast part of the northern unit is likely related to the formation of rootless cones near Elysium Mons⁴⁴. Resurfacing processes and long-term evolution of the area might have also seriously influenced the unit’s boundary in local regions.

Crater size-frequency distribution analysis⁶⁸ was conducted to estimate the absolute model ages of the identified geological units shown in Fig. 2a. The specific crater size-frequency plots for each unit are shown in Fig. 3. There were 89 ghost craters (see detailed examples in Fig. S4) identified in the lHl and lHt units. They represent relics of impact craters formed before the possible flooding in Utopia Planitia and the emplacement of VBF^51,56 in this region. Including the ghost craters and prominent craters in CSFD analysis, the estimated ages of the HNl_sm and HNl_dm units are ({3.62}_{-0.02}^{+0.02}) Ga and ({3.68}_{-0.03}^{+0.02}) Ga, respectively, which likely correspond to the AMAs of the pre-VBF surfaces in Utopia Planitia. The age ({3.68}_{-0.03}^{+0.02}) Ga estimated from HNl_dm is more statistically robust and more indicative of the southern Utopia flooding, as more ghost craters were manifested under an environment with a significant amount of water/ice in HNl_dm than in HNl_sm, but still could be underestimated due to subsequent resurfacing processes. Previous studies^14,53,79 also indicated the possible flooding time in Utopia Planitia at about 3.67–3.7 Ga. The unit lNt (renamed from lHt) has a significant amount of partly covered different-sized ghost craters with noticeable rims, which makes it possible to estimate the AMA using all crates (ghost and prominent): ({3.65}_{-0.04}^{+0.03}) Ga. This age indicates the time of formation for the shoreline in Late Noachian and can be considered the upper limit for the flooding of southern Utopia. Without including the ghost craters in the CSFD analysis, the estimated ages of the HNl_sm and HNl_dm units are ({3.5}_{-0.03}^{+0.02}) Ga and ({3.42}_{-0.03}^{+0.03}) Ga, respectively, representing the AMAs after the formation of VBF in these two units in Early Hesperian.

Crater size-frequency plot for each unit and their estimated ages. In the first two plots, the lines and text marked in black represent results that include only prominent craters. The lines and text marked in red represent results that include both prominent and ghost craters, but only the large ones are used for age estimation. The lower right is a Martian geological timescale⁷⁰ with the estimated ages of the units.

Full size image

It should be noted that the estimated ages of the units from the above AMA analysis are relatively older than those indicated on the geological map by Tanaka et al.¹³, which describes them as lHl (Late Hesperian lowland unit) and lHt (Late Hesperian transition unit). However, as described in Tanaka et al.¹³, the lHl unit represents a vast area of the northern lowlands with a similar origin and composition, and it was affected by various resurfacing processes that could alter its formation age in different parts. The scale of global mapping by Tanaka et al.¹³ did not allow for a more detailed division and age analysis of the lHl unit, which is the aim of this paper. Therefore, based on the above AMA analysis and in accordance with the naming style in Tanaka et al.¹³, we have renamed the lHl unit in the study area as HNl (Hesperian and Noachian lowland unit) and its sub-units as HNl_sm and HNl_dm. The lHt unit has also been renamed as the lNt unit (Late Noachian transition unit).

Including the ghost craters in CSFD analysis, the three units HNl_sm, HNl_dm and lNt have similar ages, all older than 3.62 Ga. Our results indicate that the emplacement of the VBF (possibly due to flooding) in Utopia Planitia may have happened at about 3.65–3.68 Ga. The two volcanic units AHv and Av have the AMAs of about ({3.28}_{-0.06}^{+0.05}) Ga and ({3.24}_{-0.13}^{+0.08}) Ga, respectively, which is consistent with their formation by lava flows/pyroclastic/lahar deposits from Elysium Mons in the east⁸⁰.

In situ observations of past water activities at Zhurong landing site

Clear mineralogic signatures of the past presence of water or ice in southern Utopia Planitia are rarely seen in the remote sensing data, such as the hyperspectral images collected by the Compact Reconnaissance Imaging Spectrometer for Mars (CRISM) onboard the Mars Reconnaissance Orbiter (MRO), due to observational biases from surface dust mantling over the region⁸¹. In situ measurements by sensors onboard the Zhurong rover could offer direct observations from the surface for the possible past presence of water in the landing area. In this study, we analyzed Zhurong rover data and discuss in situ observations in conjunction with previous related works^72,82,83,84.

Since its landing on 15 May 2021, Zhurong has traveled about 2 km from its landing location toward the south (Fig. 4). Previous analyses of the surface along the traverse of the Zhurong rover revealed a duricrust rich in hydrous minerals⁸⁵, which might be formed during the vapor–frost cycling at the atmosphere–soil interface, with the source of water from the atmosphere or groundwater/ground ice^72,84,86. The small light-toned rocks scattered on the ground could be related to duricrust weathering. The dark-toned rocks (e.g., those examined on Sols 45 and 50) have angular shapes and various sizes and could be interpreted as excavated basaltic rocks⁷². On the other hand, the massive light-toned rocks (e.g., those examined on Sols 32, 43, 65, and 79 shown in Fig. 4a) could have a sedimentary origin as part of the VBF deposits and, consequently, interacted with liquid water. We found 15 massive light-toned rocks on the NaTeCam images (Sols mentioned above and Sols 94 and 103), of which four were analyzed by MarSCoDe. The results of our spectral analysis (Fig. 4b) are generally consistent with previous works^72,84,86. The reflectance curves for the rocks measured on Sols 32, 43, 65, and 79 are characterized by absorptions at 1.94 µm and ~ 2.2 μm as well as a relatively weak absorption at 1.45 μm. Sols 32 and 65 samples have a sign of asymmetric 1.93- to 1.95-μm absorption features related to the structural H₂O⁷³. The absorption peak at 1.94 µm could be related to the stretching and bending of the water molecule and some sulfate amounts. However, the absence of noticeable spectral features at ~ 1.75 μm and 2.3 μm indicates a lack of gypsum in analyzed rocks. The absorption at ~ 2.2 μm is not narrow enough to be related to Al-phyllosilicates but could be explained by the existence of amorphous materials (hydrated silica, allophane, hydrated primary glass) and some zeolites. The mixture of hydrated silica (mainly opal or imogolite) and allophane, however, is the most plausible interpretation given the width between 2.21 and 2.26 μm⁸⁷, which conforms to the findings of Liu et al.⁷². Reflectance curves of two samples measured from the dark-toned rocks (on Sols 45 and 50) are also shown in Fig. 4b for comparison. They show absorption bands at 1.45 μm and 1.94 μm but the lack of metal-OH absorption bands. These samples are interpreted as weathered basaltic rocks.

SWIR reflectance data collected by the Zhurong rover showing rocks containing hydrated silica along the rover traverse. (a) The rover traverses through Sol 276 marked on a satellite image with the inset images showing the NaTeCam images of rocks investigated by the MarSCoDe; (b) the SWIR reflectance curves of the targets measured by the MarSCoDe; (c) Laboratory spectra of hydrated silica (opal, chalcedony) obtained from the USGS spectral library⁷³ for reference.

Full size image

However, the formation process of hydrated silica is not unique to marine systems, and they could be formed as well in contact with liquid water or weathering, hydrothermal and glacial/periglacial alteration^87,88,89. Aluminosilicates in light-toned rocks could also be formed in periglacial environments with limited water–rock interaction⁹⁰. The lack of significant absorption at 1.94 µm and ~ 2.2 μm may indicate the lack of large amounts of sulfates and phyllosilicates, implying restricted water–rock interaction and potentially arising from limited/short-time water involvement. Therefore, spectral analysis shows that the analyzed rocks could be formed in a relatively recent Amazonian environment with an atmosphere or groundwater/ground ice water source^72,84, or considering the aforementioned variety of water-related mineralogy and the morphometrical parameters of the analyzed massive light-toned rocks, it is also possible to assume that they were formed in a short-lived frozen Hesperian ocean⁷².

Xiao et al.⁸³ described rocks along the Zhurong traverse (on Sols 32, 43, 79 and 103) with well-etched faces and lamination features (see examples in Fig. 5S), typical for sedimentary environments⁹¹. In Xiao et al.⁸³’s interpretation sedimentary rocks with cross-bedding features observed at the Zhurong landing site might be formed in a medium to low-energy shallow marine environment with alternating strengths and directions of currents. Examined rocks were formed on the ocean floor and overlain by VBF material, but later were redeposited on the surface during nearby impact events. The existence of tides in a short-lived frozen Hesperian ocean is a highly controversial issue, especially considering that the estimated amplitude of tides for an open ocean on Mars is about 1 cm⁹. However, in the nearshore zone, the cumulative effect of possible currents and tides could drastically affect floor deposits.

Conceptual illustration of the three-stage evolution of the nearshore zone in southern Utopia Planitia.

Full size image

An additional possible confirmation of water-related sedimentation could come from radar data. The RoPeR on board the Zhurong rover carried out radar imaging of the subsurface along the rover traverse (down to ~ 80 m depth) using its low-frequency channel (15–95 MHz), and the radargrams showed a multi-layered subsurface structure beneath the Zhurong landing area⁸². Impact gardening and subsequent long-term weathering likely formed a fine-grained regolith upper layer no thicker than 10 m⁹². According to Li et al.⁸², the second layer extends from 10 to 30 m. It starts from small boulders at the upper part and ends with cobbles at the lower, whose clast sizes increase with depth. The third layer is significantly thicker and extends from 30 to 80 m. Clast sizes have a similar trend, increasing in particle size from the small boulders to the coarser-grained blocks. The low dielectric permittivity (3–7) of the second and third layers⁸² agrees with the estimated dielectric permittivity for VBF (~ 5) by the Mars Advanced Radar for Subsurface and Ionosphere Sounding (MARSIS)²⁴. The sorting of sediments in still waters could explain the repetitive vertical structure of the two layers with finer particles in the upper part and coarser-grained blocks in the lower part. For this scenario, the area might be influenced by two episodes of flooding sedimentation with a long water precipitation period, which translated into the formation of two VBF layers. An alternative mechanism for the formation of the two VBF layers could be the concurrent VBF formation process related to ice sublimation in the upper layer and water sedimentation in the bottom layer, which will be further discussed in the next section.

The above remote sensing observations of the Zhurong landing area provide evidence consistent with the possible existence of a single-standing water/ice reservoir. The spatial analysis allowed us to divide the southern Utopia region into a unit with water-related features less dependent on the amount of water/ice (HNl_sm) and a unit with water-related features more dependent on the amount of water/ice (HNl_dm). The age of formation of these units correlates with previous age estimations and supports the regressive ocean model. The in situ observations supplement remote sensing results, although they hardly provide direct evidence of the existence of an ancient ocean. Below, we describe and discuss our model of the nearshore zone of the short-lived ocean in Late Noachian and Early Hesperian in southern Utopia.

The new geological units (HNl_sm and HNl_dm) are consistent with the geological mapping by Tanaka^13,15 in general but provide more detailed information about the local processes and deposits in southern Utopia. The border between the lNt and HNl units lies along the Deuteronilus contact with a mean elevation of about −3580 ± 270 m¹⁷. Two different types of ghost craters were found in lNt and HNl units (Fig. S4). For the lNt unit, the different-sized ghost craters were partly covered with their rims exposed; while for the HNl units, most of the ghost craters were completely covered and manifest as up by concentric troughs. The density of ghost craters in the northern unit HNl_dm (4.35 × 10^–3 km⁻²) is higher than in the southern unit HNl_sm (1.64 × 10^–3 km⁻²). Features more dependent on the amount of water/ice, such as craters with pancake-like ejecta, polygonal troughs and etched flows^14,78 are mainly located in the northern unit HNl_dm (Figs. S1 and S2). In contrast, pitted cones of mud volcanic origin^45,46 are mostly distributed in the southern unit HNl_sm (Fig. S1), indicating a lesser presence of water^47,48,49,50. In the framework of the Martian ocean theory⁵, this spatial separation of different geomorphological features by the specific topographic contour levels (i.e., the − 3580 m contour separates the lNt and HNl units, and the − 4200 m contour separates the HNl_sm and HNl_dm units) suggests that they were formed in different marine environments that existed at different water levels at different times.

A model shown in Fig. 5 illustrates three successive stages of the evolution of a putative large water/mud body in Utopia Planitia. At first (Fig. 5a), floods filled the northern lowlands and formed a large body of water in Utopia Planitia^2,3,5,27 at about 3.65–3.68 Ga. In Utopia Planitia, the water surface reached the Deuteronilus contact at about − 3580 ± 270 m elevation¹⁷. Unit lNt represents the highland-lowland transition zone¹³ in the foreshore area and was influenced by the ocean body that flooded this area and filled the craters with sediments. In this period, foreshore and shallow marine areas could be affected by the advance and retreat of the ocean. This putative ocean would have covered units HNl_sm and HNl_dm, and dividing these units depends on the abundance of water started in the next stage.

The ocean surface was likely frozen in a geologically short period from 10⁴ to 10⁵ years⁹. Materials of the VBF formed during this time by deposition/sorting of the sedimentary load from the liquid part of the water body^2,3 and by ice sublimation in its solid part⁹. Fluid expulsion through cracks¹⁴ and pitted-cone mud volcanism^4,22,23 could also transport heavily silted water to the surface, leading to sublimation and the deposition of fine particles. Therefore, two fining-upwards layers of VBF deposits could be formed: on the ocean bed and the ice shield’s surface. The high concentration of suspended solid particles and the rapid freezing of the standing water body defined the depositional relief at the VBF border. Surface VBF deposits might decrease the ice sublimation rate^28,94 and preserve liquid water under the ice sheet. After the covering of the ocean by ice and VBF deposits, the value of subsurface infiltration of water has increased. The HNl_sm unit could have a smaller amount of liquid water covered by upper layer of VBF deposits, which was insufficient to form etched flows and pancake-like ejecta. A lower amount of water/ice did not form sufficient surface subsidence due to volatile infiltration, which did not allow the formation of a significant number of polygonal and circular troughs and resulted in fewer ghost craters. However, the abundance of pitted cones with a mud volcanic origin still indicates a small amount of volatiles necessary for their formation²².

On the ocean floor, water-related minerals start forming, such as the hydrated silica identified by the MarSCoDe (Fig. 4b). The variety of minerals in the rock samples indicates the plausible presence of a short-lived frozen ocean. The sedimentary rocks with cross-bedding features described by Xiao et al.⁸³ could be formed at the first stage of water body evolution before the complete freezing. The massive light-toned rocks examined by the MarSCoDe could have been redeposited on the surface during nearby impact events and may represent portions of the VBF that were exhumed from depth⁷².

On the contrary, the HNl_dm unit could have more abundant water sources, which maintained a prolonged supply of volatile materials and allowed larger-scale features such as polygonal troughs and etched flows to form. Materials with more significant amounts of volatiles (water) also permitted the preferential formation of craters with pancake-like ejecta in this unit, whereas such craters were less frequent in unit HNl_sm. Therefore, the − 4200 m contour could separate the shallower marine environment (HNl_sm) of less abundant water/ice and the deeper marine environment (HNl_dm) of more abundant water reservoirs. Gradually, liquid water solidified and/or the ice sheet disappeared from unit HNl_sm and later from unit HNl_dm. Based on CSFD analyses, the AMAs of units HNl_sm and HNl_dm are ({3.5}_{-0.03}^{+0.02}) Ga and ({3.42}_{-0.03}^{+0.03}) Ga, respectively. Particularly, the younger age of unit HNl_dm indicates the extended duration of water-related influence and resurfacing as this unit formed. The AMA analysis of the coastline and nearshore zone of southern Utopia implies the existence of a short-lived ocean in Late Noachian and Early Hesperian in this area. This result is comparable to previous publications (e.g.,^17,57) but provides a more accurate dating of the nearshore zone. It is likely that the central part of Utopia preserved water-content deposits for a longer time, until the Late Hesperian¹⁴.

Close to the final stage of the nearshore zone evolution (Fig. 5c), the freezing of the ocean remnants under the VBF deposits could have formed either a single body of ice or a series of lens-like bodies; some of which could be responsible for the formation of features observed on the surface of Utopia Planitia, such as textured terrain or scalloped terrain (e.g.,⁹⁵). The subsurface ice lenses could persevere into the Amazonian period and perhaps contribute to the formation of surface duricrusts enriched in hydrous minerals^84,86.

The presence of VBF material in the study area plays an essential role in the ocean theory⁵. The ~ 10 m upper regolith layer represents the fine-grained material formed during Amazonian weathering. Below, two fining-upwards layers with similar vertical structures (finer particles in the upper part and coarser-grained blocks in the lower part) were found by the RoPeR radar data along the rover traverse, interpreted as VBF deposits⁸². These two subsurface layers might have been formed by sedimentation during the advance and retreat of the water body,however, in this case, it remains unclear the role of deposits accumulated on the frozen ocean surface in the stratigraphy discovered by the radar. The most probable scenario is a concurrent VBF formation process whereby the lower layer formed as the sediment load from water/mud ocean body^2,3 settled and stored, and the upper layer formed by ice sublimation⁹ during the second stage (Fig. 5b) of the ocean evolution.

The Utopia Planitia region on Mars exhibits water-related features on its surface. The AMA, geomorphological and spatial analysis at the Zhurong landing area in southern Utopia allowed us to divide the region (lHl) into two units: one with water-related features less dependent on the amount of water/ice (HNl_sm) and one with water-related features more dependent on the amount of water/ice (HNl_dm). The HNl_sm unit is located at the southern border of Utopia Planitia and is characterized by pitted cones of mud volcanic origin^45,46, indicating a lesser presence of water/ice^47,48,49,50. Northwards lies the HNl_dm unit, which has an increased concentration of craters with pancake-like ejecta, polygonal troughs, and etched flows, indicating a greater presence of water^14,78. In situ measurements by sensors onboard the Zhurong rover hardly provide direct evidence of the existence of an ancient ocean; however, they do not contradict this model and could be considered within the context of ocean theory. The water-related minerals, such as hydrated silica (mainly opal and imogolite) and allophane, identified by MarSCoDe, could have formed under conditions unrelated to the ocean theory, but they could also have formed in the short-lived frozen Hesperian ocean. Cross-bedding features interpreted by Xiao et al.⁸³ as sedimentary rocks affected by currents and tides are controversial but could have formed after the flooding and before ice shield formation. The layering structure found by RoPeR radar data analysis^82,92 could be interpreted as VBF deposits formed during ice sublimation and water sedimentation.

Geomorphological, spatial, and AMA analyses in this study suggest an evolution model for the nearshore zone of southern Utopia Planitia:

1. 1.

   The shallow (HNl_sm) and deep (HNl_dm) marine environments formed by the flooding of Utopia Planitia in ({3.65}_{-0.04}^{+0.03})–({3.68}_{-0.03}^{+0.02}) Ga in Late Noachian. The water reached the Deuteronilus contact at about −3580 ± 270 m elevation¹⁷ and formed a coastline.

2. 2.

   The ocean surface was likely frozen in a geologically short period, and VBF deposits and water-related geomorphological features started to form. The VBF material is deposited by sedimentary load from the liquid part of the water body^2,3 and by ice sublimation in its solid part⁹, forming the layering structure of deposits. The spatial separation of the water-related features by the −4200 m contour indicates different water amounts, with less water/ice in the HNl_sm unit and more water/ice in the HNl_dm unit.

3. 3.

   Gradually, liquid water solidified, and/or the ice sheet disappeared, completing the formation of HNl_sm in ({3.5}_{-0.03}^{+0.02}) Ga and HNl_dm in ({3.42}_{-0.03}^{+0.03}) Ga units. Therefore, the water body’s existence was dated from about 3.65 until 3.42 Ga (Late Noachian and Early Hesperian). The subsurface ice remnants in a single ice body or separated lens format could persevere into the Amazonian period and could be responsible for forming scalloped terrain or surface duricrusts in Utopia Planitia.

It should be noted that the Zhurong rover landed in the HNl_sm unit, approximately 20 km south of the boundary between the HNl_sm and HNl_dm units, and traveled further south after landing. This means the rover will not be able to collect in situ data in the northern HNl_dm unit to directly verify the proposed model. Nevertheless, possible future missions to land in the HNl_dm unit or orbital observations (e.g., radar remote sensing) focused on the comparison between the HNl_sm and HNl_dm units will help to further verify the proposed model.