Letterdoi:10.1038/nature26144

Timing of oceans on Mars from shoreline
deformation
Robert I. Citron1,2, Michael Manga1,2 & Douglas J. Hemingway1,2

Widespread evidence points to the existence of an ancient Martian
ocean1–8. Most compelling are the putative ancient shorelines in
the northern plains2,7. However, these shorelines fail to follow
an equipotential surface, and this has been used to challenge thenotion that they formed via an early ocean9 and hence to question
the existence of such an ocean. The shorelines’ deviation from a
constant elevation can be explained by true polar wander occurring
after the formation of Tharsis10, a volcanic province that dominates
the gravity and topography of Mars. However, surface loading fromthe oceans can drive polar wander only if Tharsis formed far fromthe equator10, and most evidence indicates that Tharsis formed near
the equator11–15, meaning that there is no current explanation forthe shorelines’ deviation from an equipotential that is consistent
with our geophysical understanding of Mars. Here we show thatvariations in shoreline topography can be explained by deformation
caused by the emplacement of Tharsis. We find that the shorelines
must have formed before and during the emplacement of Tharsis,
instead of afterwards, as previously assumed. Our results imply that
oceans on Mars formed early, concurrent with the valley networks15,
and point to a close relationship between the evolution of oceans
on Mars and the initiation and decline of Tharsis volcanism, withbroad implications for the geology, hydrological cycle and climateof early Mars.
Distinct geological boundaries (contacts) lining the northern
plains of Mars for thousands of kilometres have been interpreted
as ­palaeo-shorelines and evidence of an early ocean2–4,6,7. However,observed long-wavelength deviations (by up to several kilometres) in
shoreline elevation from an equipotential have been used as an argument against the emplacement of the contacts by a body of liquid water,
the interpretation of the features as shorelines, and the existence of a
Martian ocean9. Perron et al.10 showed that the elevation changes of two
extensive contacts, Arabia (contact 1) and Deuteronilus (contact 2), canbe explained by deformation due to 30°–60° and 5°–25° of post-Tharsis
true polar wander (TPW), respectively, because a varying rotation polealso changes the orientation of a planet’s equatorial bulge, or polar flattening, altering equipotential surfaces (such as sea levels) globally. Suchlarge magnitudes of TPW can be driven by ocean loading/unloading,but only if Tharsis formed far from the equator10. If Tharsis formednear the equator, then the remnant fossil bulge would have prevented
ocean loading from causing large amounts of post-Tharsis TPW (seeExtended Data Fig. 1).
Most evidence points to the formation of Tharsis near the­equator11–15. Mars’ remnant rotational figure (fossil bulge) is close to
17.5
the equator, indicating a palaeopole of (259.5 ±​  49.5° E, 71.1−+14.4 ° N),14
the likely pre-Tharsis orientation of Mars . The pre-Tharsis palaeopole
also matches the likely orientation of Mars during valley network­formation15. Formation of Tharsis probably drove only limited(approxi­mately 20°) TPW to reach Mars’ current configuration, which
precludes the possibility that surface loads drove sufficient TPW to
deform the shorelines10,16.

We propose that the Arabia shoreline instead formed before or
during the early stages of Tharsis emplacement, which initiated >​  3.7billion years (Gyr) ago17 when the rotation pole of Mars was at thepalaeopole (259.5° E, 71.1° N) corresponding to the fossil bulge14.
The Arabia shoreline, potentially emplaced at least 4 Gyr ago6, would
have been modified by both topographic changes from Tharsis (which­dominates Mars’ topography and gravity on a global scale; see Extended
Data Fig. 2), and the approximately 20° of Tharsis-induced TPW. The
Deuteronilus shoreline, which differs less from a present-day equipotential than the older Arabia shoreline, is dated to about 3.6 Gyrago18, after most of Tharsis was emplaced. However, Tharsis had­complex and multi-stage growth that extended into the Hesperian andAmazonian17,19, meaning that the Deuteronilus shoreline would havebeen deformed by the late stages of Tharsis’ emplacement. We examine
a chronology in which shoreline deformation is due mainly to Tharsis(Table 1), and compare expected deformation due to Tharsis with theelevation profiles of the Arabia and Deuteronilus contacts.
Assuming the Arabia shoreline formed before Tharsis, and theDeuteronilus shoreline formed after most of Tharsis was emplaced,
we compare the best fits for the deformation expected from Tharsis
to the current topography of the shorelines, including an offset factor
Z to represent sea level at the time of shoreline formation. We also
examine the Isidis shoreline, which formed 100 million years (Myr)after Deuteronilus18. For the Arabia shoreline emplaced before Tharsis,
deformation is expressed as the contribution of Tharsis to Mars’ topo­
graphy along the shoreline, and the change in topography from l­ imited
Tharsis-induced TPW. For the Deuteronilus and Isidis shorelines
emplaced during the late stages of Tharsis growth, deformation is taken
as the percentage of Tharsis’ contribution to topography occurring after
the shorelines formed, and no contribution from TPW (because reorientation should occur within tens of thousands of years to a few million
years after the Tharsis plume reaches the surface20, much less than the
100 Myr or more that lies between Tharsis initiation and Deuteronilus
formation). See Methods for more details.
We show that the Arabia shoreline’s deviations from an equipotential can be explained almost entirely by deformation due to Tharsis
emplacement (Fig. 1). Our best fit (equation (3) with Z =​  −​2.3 km)
yields a root-mean square misfit σrms of 0.615 km, comparable to theerror values from Perron et al.10, and follows the slope of the shoreline
data better from 1,000 km to 6,600 km. The limited Tharsis-induced
TPW has a negligible effect. A slightly lower σrms is obtained if only 80%
of Tharsis topography was emplaced after the Arabia shoreline formed(Extended Data Fig. 3). However, the difference between the fits using
80% or 100% of Tharsis’ topography is negligible considering the s­ catterin the shoreline data. Our model therefore suggests that the Arabia
shoreline formed before or during the early stages of Tharsis’ growth.
The Deuteronilus shoreline’s deviations from an equipotential canbe explained by deformation due to the emplacement of about 17%of Tharsis topography (Fig. 2), indicating that the shoreline formed
during the late stages of Tharsis’ growth. Our best fit (equation (4) with1
Department of Earth and Planetary Science, University of California, Berkeley, Berkeley, California, USA. 2Center for Integrative Planetary Science, University of California, Berkeley, Berkeley, California,
USA.

0 0 M o n t h 2 0 1 8 | VO L 0 0 0 | NAT U R E | 1© 2018 Macmillan Publishers Limited, part of Springer Nature. All rights reserved.

RESEARCH Letter
Table 1 | Possible evolution of Martian shorelines
Time (epoch)

Event

At least 4 Gyr ago (Early Noachian)>​3.6 Gyr ago (Late Noachian/Early Hesperian)Arabia shoreline forms
Majority of Tharsis forms

Effect

3.6 Gyr ago (Early Hesperian)
3.5 Gyr ago (Late Hesperian)
3.6 Gyr ago to present (Early Hesperian−Late Amazonian)Deuteronilus shoreline forms
Isidis shoreline forms
Remainder of Tharsis forms

3.5–3.0 Gyr ago(Late Hesperian/Early Amazonian)

Loading of Utopia basin

Both Tharsis emplacement/loading and Tharsis-induced
TPW (20°) deform the Arabia shoreline

Remaining Tharsis growth deforms the Arabia,
Deuteronilus and Isidis shorelines
Loading of the Utopia basin tilts the Isidis basin anddeforms sections of the Deuteronilus shoreline

The table outlines the sequence of events required to produce a pre- or early-Tharsis Arabia shoreline, and the formation of the Deuteronilus shoreline concurrent with the later stages of Tharsis’growth. Explaining the Isidis shoreline’s topography requires subsequent loading of the Utopia basin (see Extended Data Fig. 4 and Methods).

C =​  0.17 and Z =​  −​3.68 km) yields σrms =​ 0.110 km. Our fit successfully
recovers the low elevation of the Phlegra contact, and also captures
the decrease in elevation across Utopia and Elysium West. Neither ourmodel nor the Perron et al.10 model captures the full elevation increaseof the Tantalus segment, which may result from the topographic bulge
from nearby Alba Patera18. For the Isidis shoreline, subsequent loadingof the Utopia basin is also required to explain the shoreline’s topography(see Extended Data Fig. 4 and Methods).
The relation between the shorelines and global deformation due
to Tharsis and its associated TPW is illustrated in Fig. 3a–c (also see
Extended Data Fig. 2). We estimate the volume of water that filled thenorthern plains to the Deuteronilus and Arabia shorelines by subtracting
the relevant Tharsis and TPW contributions from Mars’ topography(0.25° per pixel gridded MOLA data21) and filling the lowlands toshoreline elevation (Fig. 3d–f). We estimate a Deuteronilus ocean
volume of about 1.2 ×​  107 km3, and an Arabia ocean volume of about
4.1 ×​  107 km3. These are lower limits because we do not remove excess
terrain, such as Elysium, polar deposits, lava/sediment basin deposits,
and short-wavelength Tharsis topography (that is, variations in Tharsis
topography that occur over short length scales). For the Arabia ocean,use of a map of Mars with excess terrain removed15 yields an oceanavolume of about 5.5 ×​  107 km3. The ocean volumes we compute areslightly lower than previous estimates22 because the Tharsis topography
we subtract is negative in much of the area enclosed by the northernocean basin.
Short-wavelength deviations in shoreline elevation from our model
may be due to our assumptions that both lithospheric thickness and therate of Tharsis emplacement are spatially uniform. Spatial variations
in lithospheric thickness23 would allow for non-uniform responses topheno­mena such as TPW10, plate flexure and dynamic topography
from the Tharsis plume24. Spatially variable Tharsis emplacement could
also affect shoreline modification. Another consideration is ocean
loading, but the computed effect on shoreline elevations is small (seeExtended Data Fig. 5 and Methods).
Several other short- and long-wavelength processes could have
deformed the shorelines in the 3.5 Gyr or more since their emplacement, including dynamic topography from mantle plumes 24,­lithospheric deformation22,25, glacial erosion26 and plate flexure fromloading/unloading. For example, loading of the Utopia basin may havetilted Isidis (see Methods) and deformed sections of the Deuteronilus
shoreline18. Other loads that post-date shoreline formation include
Elysium, the polar deposits, and sections of Tharsis. Such loads

1.0
ΔTTharsis (km)ΔTTPW (km)0.5

ΔT (km)

0.0
–0.5
–1.0
–1.5
–2.0

b

0
Current shoreline topography
Perron et al. (2007) model (ref. 10)
Present paper

Elevation (km)

–1
–2
–3
–4
0

1,000

2,000

3,000

4,000

5,000

6,000

7,000

Distance along shoreline (km)Figure 1 | Comparison of Arabia shoreline topography to shoreline
deformation models. a, Change in topography ∆T caused by TPW of 20°(equation (1)) and Tharsis uplift (equation (2)), illustrating that the latteris much more important. b, Current topography of the Arabia shorelinefrom Perron et al.10 (data originally from ref. 7) compared to the Perron

et al.10 model of deformation due to post-Tharsis TPW (with Te =​ 200 km)
and our model of deformation due to Tharsis emplacement and induced
TPW (∆T Tharsis +​ ∆T TPW −​2.3 km). The starting point for the shorelineis (24.91° W, 13.48° N).

2 | NAT U R E | VO L 0 0 0 | 0 0 m o n t h 2 0 1 8© 2018 Macmillan Publishers Limited, part of Springer Nature. All rights reserved.

Letter RESEARCH–3.0
Tantalus–3.2
Elevation (km)Current shoreline topography
Perron et al. (2007) model (ref. 10)
Present paper–3.4
Cydonia–
Deuteronilus

Tempe–Chryse–Acidalia

–3.6
–3.8

PhlegraPyramus–Astapus–
Utopia–Elysium W

–4.0
–4.2

0

5,000

10,00015,00020,00025,00030,000Distance along shoreline (km)Figure 2 | Comparison of Deuteronilus shoreline topography toshoreline deformation models. Current Deuteronilus topography (dataand contact names from ref. 18) compared to the Perron et al.10 modeland our model of deformation due to partial Tharsis emplacement(0.17∆T Tharsis −​3.68 km). The starting point for the shoreline is(96.40° W, 63.69° N).

could also induce a small amount (<​2°) of post-Tharsis TPW16. Plate­flexure associated with impact basins could also deform shorelines.
While basins of over 1,000 km in diameter pre-date the Deuteronilus
shoreline, some basins may have been coincident with or post-date theArabia shoreline. Short-wavelength deformation may also be a consequence of the difficulty in identifying the shorelines themselves18.
Increased accuracy in shoreline identification and dating18 can help
to reconstruct the history of shoreline formation and modification.
Several potential shorelines6, such as the Ismenius, Acidalia, Elysium,
Aeolis and Meridiani contacts, have been relatively unexaminedowing to their high degree of discontinuity7. Shorelines may also be
mismapped; for example, portions of the Meridiani shoreline may
be part of the Arabia shoreline22. A re-evaluation of shorelines withfull consideration of the various deformation processes may enable
the development of a chronology of oceans on Mars. In particular, theMeridiani shoreline6,27 may pre-date the Arabia shoreline and havecontained a larger volume of water22.
Accurate dating of the Arabia shoreline is necessary to determine
whether the shoreline formed before or during the early stages ofTharsis’ growth. Formation of the Arabia shoreline after some l­ imited

early Tharsis growth is suggested by Arabia segments that border
Acheron Fossae and Tempe Terra6, two of the oldest Tharsis units, which
are located well north of the expected pre-Tharsis crustal dichotomy
boundary (the stark difference in elevation and crustal thickness
between the northern lowlands and southern highlands). However,
it is possible that the Archeron Fossae and Tempe Terra contacts weremisidentified as belonging to the Arabia shoreline, or that the Arabia
shoreline initially followed the pre-Tharsis dichotomy boundary, andformed the Tempe Terra and Archaeon Fossae contacts only after earlyTharsis uplift and deposition.
The decline in ocean volume from the pre- or early-Tharsis Arabia
shoreline to the late-Tharsis Deuteronilus shoreline suggests thatTharsis volcanism may have played a critical part in the evolution of a
Martian ocean. After Tharsis was mostly emplaced, by about 3.6 Gyrago, only short-lived lakes may have been stable28, although a Late
Hesperian/Early Amazonian ocean has also been suggested on thebasis of tsunami evidence8,29. Outgassing from Tharsis could havecontributed to either heating30 or cooling31 the planet, both of which
could produce the decrease in ocean volume from the Arabia shoreline to the Deuteronilus shoreline as Tharsis activity declined. Either aab

MOLA topographykm

Tharsis contribution
to topography

c

12
9
6
3
0
−3
−6

d

Deuteronilus Ocean

km
0.24

4.5

e
km

Deformation from TPW

km

Arabia Ocean

3.0

0.12

1.5

0.00

0.0

−0.12

−1.5

−0.24

f
km

Arabia Ocean(excess topography removed)

km

12

6

6

8

4

4

2

2

0

0

−2

−2

4
0
−4
−8

Figure 3 | Shoreline locations relative to current topography,
deformation due to Tharsis/TPW, and computed ocean extents.
a, MOLA topography. b, Tharsis contribution to topography (equation(2)). c, Deformation due to approximately 20° TPW from the palaeopole
corresponding to the fossil bulge14 with Te =​ 58 km (equation (1)).
d, Ocean basin filled to the Deuteronilus shoreline, with the topographyof Mars at the time of the Deuteronilus shoreline’s formation (MOLA
topography minus 17% of Tharsis topography). e, Ocean basin filled to the
Arabia shoreline, with the topography of Mars at the time of the Arabia

shoreline’s formation (MOLA topography minus Tharsis topography anddeformation due to TPW). Short-wavelength remnants of Tharsis arevisible because Tharsis topography is only modelled up to degree-5.
f, The ‘Mars without Tharsis’ map15, which is similar to e, but with shortwavelength Tharsis, Elysium, and polar topography removed, filled to the
Arabia shoreline. Shorelines are plotted for Deuteronilus (white), Arabia(magenta) and Isidis (cyan), and colour scales denote elevation (or changes
in elevation) in kilometres.
0 0 M o n t h 2 0 1 8 | VO L 0 0 0 | NAT U R E | 3© 2018 Macmillan Publishers Limited, part of Springer Nature. All rights reserved.

RESEARCH Letter
large ocean was in place before Tharsis volcanism initiated, and shrankas Tharsis volcanism cooled the planet, or an ocean arose as a resultof heating caused by Tharsis outgassing and decreased in volume asTharsis volcanism declined. It is also possible that each shoreline represents the transient warming of an otherwise frozen ocean or g­ lacialstate30, producing a liquid ocean in periods of heightened Tharsis­activity (which, owing to enhanced surface heat flux, may also haveresulted in catastrophic circum-Tharsis groundwater discharge1). If
episodic warming was sufficient to melt most of Mars’ glaciers, thedecrease in ocean volume may record a declining surface water budget.
Although geochemical evidence for a northern ocean is ambiguous32,
an ocean supported by the degassing of sulfur from Tharsis could
explain the lack of widespread carbonate deposits observed in thenorthern plains33.
The evolution of water on Mars is critical to understanding the past
climate and habitability of the planet. Although shorelines on Mars haveprovided compelling evidence for a Martian ocean, to explain theirdeviations from an equipotential has been a challenge. We show thatthe topography of Martian shorelines can be quantitatively explainedby deformation due to the emplacement of Tharsis and resulting TPW(in the case of the Arabia shoreline) or by the latter stages of Tharsis
emplacement (in the case of the Deuteronilus shoreline). Formationof the Arabia shoreline before (or during the early stages of) Tharsis
emplacement suggests that the Arabia ocean was concurrent with valley
network incision15, which probably occurred as part of a globally activehydrosphere capable of supporting such an ocean5. The consistency
between the topography of the Martian shorelines, their ages, and thechronology of topographic changes due to Tharsis emplacement andassociated TPW, suggests that the Arabia and Deuteronilus contacts
are evidence that Martian oceans existed, and may have been linked
to Tharsis volcanism.
Online Content Methods, along with any additional Extended Data display items andSource Data, are available in the online version of the paper; references unique tothese sections appear only in the online paper.
received 1 October 2017; accepted 25 January 2018.
Published online 19 March 2018.
1.

Baker, V. R. et al. Ancient oceans, ice sheets and the hydrological cycle on Mars.
Nature 352, 589–594 (1991).
2. Parker, T. J., Saunders, R. S. & Schneeberger, D. M. Transitional morphology inWest Deuteronilus Mensae, Mars: implications for modification of the lowland/
upland boundary. Icarus 82, 111–145 (1989).
3. Parker, T. J., Gorsline, D. S., Saunders, R. S., Pieri, D. C. & Schneeberger, D. M.
Coastal geomorphology of the Martian northern plains. J. Geophys. Res. Planets
98, 11061–11078 (1993).
4. Head, J. W. et al. Possible ancient oceans on Mars: evidence from Mars Orbiter
Laser Altimeter data. Science 286, 2134–2137 (1999).
5. Di Achille, G. & Hynek, B. M. Ancient ocean on Mars supported by global
distribution of deltas and valleys. Nat. Geosci. 3, 459–463 (2010).
6. Clifford, S. M. & Parker, T. J. The evolution of the Martian hydrosphere:
implications for the fate of a primordial ocean and the current state of thenorthern plains. Icarus 154, 40–79 (2001).
7. Carr, M. H. & & Head, J. W. Oceans on Mars: an assessment of the observational
evidence and possible fate. J. Geophys. Res. Planets 108, 5042 (2003).
8. Rodriguez, J. A. P. et al. Tsunami waves extensively resurfaced the shorelines ofan early Martian ocean. Sci. Rep. 6, 25106 (2016).
9. Malin, M. C. & Edgett, K. S. Oceans or seas in the Martian northern lowlands:
high resolution imaging tests of proposed coastlines. Geophys. Res. Lett. 26,
3049–3052 (1999).
10. Perron, J. T., Mitrovica, J. X., Manga, M., Matsuyama, I. & Richards, M. A.
Evidence for an ancient martian ocean in the topography of deformed
shorelines. Nature 447, 840–843 (2007).
11. Willemann, R. J. Reorientation of planets with elastic lithospheres. Icarus 60,
701–709 (1984).
12. Roberts, J. H. & Zhong, S. The cause for the north-south orientation of thecrustal dichotomy and the equatorial location of Tharsis on Mars. Icarus 190,
24–31 (2007).

13. Daradich, A. et al. Equilibrium rotational stability and figure of Mars. Icarus
194, 463–475 (2008).
14. Matsuyama, I. & Manga, M. Mars without the equilibrium rotational figure,Tharsis, and the remnant rotational figure. J. Geophys. Res. Planets 115,
E12020 (2010).
15. Bouley, S. et al. Late Tharsis formation and implications for early Mars. Nature
531, 344–347 (2016).
16. Kite, E. S., Matsuyama, I., Manga, M., Perron, J. T. & Mitrovica, J. X. True Polar
wander driven by late-stage volcanism and the distribution of paleopolardeposits on Mars. Earth Planet. Sci. Lett. 280, 254–267 (2009).
17. Anderson, R. C. et al. Primary centers and secondary concentrations of tectonic
activity through time in the western hemisphere of Mars. J. Geophys. Res.
Planets 106, 20563–20585 (2001).
18. Ivanov, M. A., Erkeling, G., Hiesinger, H., Bernhardt, H. & Reiss, D. Topography of
the Deuteronilus contact on Mars: evidence for an ancient water/mud ocean
and long-wavelength topographic readjustments. Planet. Space Sci. 144,
49–70 (2017).
19. Dohm, J. M., Baker, V. R., Maruyama, S. & Anderson, R. C. In Superplumes:
Beyond Plate Tectonics 523–536 (Springer, 2007).
20. Rouby, H., Greff-Lefftz, M. & Besse, J. Rotational bulge and one plumeconvection pattern: influence on Martian true polar wander. Earth Planet. Sci.
Lett. 272, 212–220 (2008).
21. Smith, D. E., Zuber, M. T., Neumann, G. A., Guinness, E. A. & Slaveny, S. Mars
Global Surveyor Laser Altimeter mission experiment gridded data record.
MGS-M-MOLA-5-MEGDR-L3–V1.0, http://pds-geosciences.wustl.edu/missions/
mgs/megdr.html (NASA Planetary Data System, 2003).
22. Ruiz, J., Tejero, R., Gómez-Ortiz, D. & López, V. in Space Science: New Research
141–164 (Nova Science, 2006).
23. Grott, M. & Breuer, D. On the spatial variability of the Martian elastic
lithosphere thickness: evidence for mantle plumes? J. Geophys. Res. Planets
115, E03005 (2010).
24. Roberts, J. H. & Zhong, S. Plume-induced topography and geoid anomalies
and their implications for the Tharsis rise on Mars. J. Geophys. Res. Planets
109, E03009 (2004).
25. Ruiz, J., Fairén, A. G., Dohm, J. M. & Tejero, R. Thermal isostasy anddeformation of possible paleoshorelines on Mars. Planet. Space Sci. 52,
1297–1301 (2004).
26. Davila, A. F. et al. Evidence for Hesperian glaciation along the Martiandichotomy boundary. Geology 41, 755–758 (2013).
27. Edgett, K. S. & Parker, T. J. Water on early Mars: possible subaqueous
sedimentary deposits covering ancient cratered terrain in western Arabia andSinus Meridiani. Geophys. Res. Lett. 24, 2897 (1997).
28. Kite, E. S. et al. Methane bursts as a trigger for intermittent lake-forming
climates on post-Noachian Mars. Nat. Geosci. 10, 737–740 (2017).
29. Costard, F. et al. Modeling tsunami propagation and the emplacement of
thumbprint terrain in an early Mars ocean. J. Geophys. Res. Planets 122,
633–649 (2017).
30. Halevy, I. & Head, J. W. Episodic warming of early Mars by punctuatedvolcanism. Nat. Geosci. 7, 865–868 (2014).
31. Tian, F. et al. Photochemical and climate consequences of sulfur outgassing on
early Mars. Earth Planet. Sci. Lett. 295, 412–418 (2010).
32. Pan, L., Ehlmann, B. L., Carter, J. & Ernst, C. M. The stratigraphy and history ofMars’ northern lowlands through mineralogy of impact craters: acomprehensive survey. J. Geophys. Res. Planets 122, 1824–1854 (2017).
33. Halevy, I., Zuber, M. T. & Schrag, D. P. A sulfur dioxide climate feedback on earlyMars. Science 318, 1903–1907 (2007).
Supplementary Information is available in the online version of the paper.
Acknowledgements We thank J. T. Perron for providing the data for the Arabia
shoreline (originally from ref. 7), and M. A. Ivanov for providing the data for
the Deuteronilus and Isidis shorelines. We thank I. Matsuyama for discussions
regarding this research. R.I.C. and M.M. are supported by NSF EAR-1135382.
D.J.H. is supported by the Miller Institute for Basic Research in Science.
Author Contributions All authors discussed the research idea, methods andinterpretation of results. R.I.C. and M.M. developed the hypothesis with input
from D.J.H. R.I.C. performed the calculations and wrote the manuscript, withguidance, comments, and revisions from M.M. and D.J.H.
Author Information Reprints and permissions information is available atwww.nature.com/reprints. The authors declare no competing interests.
Readers are welcome to comment on the online version of the paper.
Publisher’s note: Springer Nature remains neutral with regard to jurisdictional
claims in published maps and institutional affiliations. Correspondence andrequests for materials should be addressed to R.I.C. (ricitron@berkeley.edu).
Reviewer Information Nature thanks S. Bouley and M. Zuber for theircontribution to the peer review of this work.

4 | NAT U R E | VO L 0 0 0 | 0 0 m o n t h 2 0 1 8© 2018 Macmillan Publishers Limited, part of Springer Nature. All rights reserved.

Letter RESEARCH
Methods

Arabia shoreline (pre- or early-Tharsis formation). We assume deformation of
the Arabia shoreline since its formation is due to global changes in topography
resulting from Tharsis’ formation (emplacement and loading) and the approximately 20° of Tharsis-induced TPW (Table 1).
The topographic response to TPW is given by the change in the flattening of theplanet caused by the difference between the centrifugal potential at the initial andfinal rotation poles10. For a shoreline in place before TPW occurs, the deformationof the shoreline topography due to TPW10 is:∆TTPW(θ , ϕ) =

ω 2a 2[P 2,0(cosγ) − P 2,0(cosθ)][h2 − (1 + k2)]
3g

(1)

where a is the mean planetary radius, ω is the rotation rate, g is the surface gravity,
γ is the angular distance between a given current colatitude and longitude (θ, ϕ)and the palaeopole and h2 and k2 are the secular (fluid-limit) degree-2 Love numbers that depend on the density and elastic structure of Mars. The unnormalized
degree-2 Legendre polynomial P2,0(cosη) =​  1 (3cos2η −​ 1).
2
The change in topography due to the emplacement of Tharsis and its associatedloading is:∆TTharsis(θ , ϕ) = S Tharsis(θ , ϕ) − NTharsis(θ , ϕ)(2)where STharsis and NTharsis are Tharsis’ contribution to the shape and geoid of Mars,
respectively. We use gravity and shape coefficients for Tharsis up to degree-5 from
Matsuyama and Manga14.
The current topography of the Arabia shoreline should therefore follow the
deformation profile ∆TA, given by:∆TA(θ , ϕ) = ∆TTPW(θ , ϕ) + ∆TTharsis(θ , ϕ) + Z(3)

where Z is a constant to adjust for the sea level at the time of the shoreline’s
emplacement. We minimize the least-squares misfit (σrms) between equation (3)and the shoreline elevation data for the Arabia contact examined in ref. 10 (dataoriginally from ref. 7). We assume a fixed palaeopole (259.5° E, 71.1° N), corresponding to the fossil bulge14. We use an elastic lithosphere thickness Te =​ 58 km,
the expected value at the time of Tharsis’ emplacement14, corresponding to h2 =​  2.0and k2 =​  1.115.
We also test whether the Arabia shoreline can be explained by deformation
due to only a certain percentage of Tharsis’ emplacement and associated ­loading,by multiplying ∆T Tharsis in equation (3) by a factor C, corresponding to thepercentage of Tharsis topography emplaced after shoreline formation (seeExtended Data Fig. 3).
Deuteronilus shoreline (late-stage Tharsis formation). The Deuteronilus shoreline post-dates the initiation of Tharsis by >​100 Myr, and therefore probably
formed after a large portion of Tharsis was emplaced. The shoreline also probably
post-dates most Tharsis-induced TPW, which should have occurred within a few
million years of load emplacement20. Estimates of load-driven TPW on Mars suggest timescales less than 10 Myr34–36, well within the required pre-Deuteronilus
timescale. Although a fraction of the 20° of Tharsis-induced TPW may be due torelaxation of the lithosphere and occur on longer timescales37, this should havea negligible effect given the small influence of TPW on shoreline deformation(Fig. 1a). Accordingly, we assume that deformation to the Deuteronilus shoreline since its formation is due to the topographic response of Mars to only the latestages of Tharsis’ emplacement and associated loading. The current topography of theDeuteronilus shoreline should therefore follow the deformation profile ∆TD, given by:∆TD(θ, ϕ) = C ∆TTharsis(θ , ϕ) + Z(4)where ∆T Tharsis is deformation due to Tharsis (equation (2)), C is a constantrepresenting how much of Tharsis formed after Deuteronilus was formed, andZ is a constant to adjust for sea level at the time of the shoreline’s formation. We
minimize the misfit between equation (4) and the Deuteronilus shoreline elevation
data18, to determine the optimal amount of Tharsis topography that should postdate the shoreline’s formation.
Isidis shoreline (late-stage Tharsis formation). Because the Isidis shoreline is
100 Myr younger than Deuteronilus, we assume a similar deformation profile andcompare equation (4) to the Isidis shoreline data18. We use the same value of C optimized for the Deuteronilus shoreline, but allow Z to vary, reflecting that sea level
could change in the 100 Myr between the formation of the Deuteronilus contact
and the Isidis contact, but deformation from Tharsis topography should not change
substantially. Although C should be slightly less for Isidis, optimizing for C would
result in an unrealistic C =​ 0 because deformation due to Tharsis along the Isidis
shoreline results in a tilt that is opposite to the present tilt (Extended Data Fig. 4).

For the Isidis shoreline, our model predicts that deformation due to Tharsis
would have tilted Isidis opposite to its present tilt (Extended Data Fig. 4). While
this appears contradictory, the mismatch is possible if Isidis was tilted to its present
orientation by loading of the Utopia basin. The Utopia basin has a large positive
gravity anomaly38,39, indicating about 18 km of excess fill40. Such a load would havecaused elastic plate flexure and a peripheral bulge, which could have tilted the Isidis
basin. Using a plate flexure model, McGowan and McGill41 show that loading of
Utopia could have tilted Isidis to an even greater extent than currently observed.
Therefore, some amount of reverse tilting (as our model predicts) is possible. The
timing of Utopia loading relative to the subsequent Tharsis deformation is irrelevant provided that Utopia loading also occurred after the Isidis shoreline formed.
We expect loading of Utopia to occur after Isidis shoreline formation because ashrinking Martian ocean would evaporate from the Utopia basin last, depositing
the non-volatile component of the ocean there. Additionally, if the ocean became
cold and glacial during its decline42, then receding glaciers may also have loaded
Utopia with excess sediment. The deposits in the base of Utopia basin date to
the early Amazonian (<​3–3.46 Gyr ago)43,44, after the emplacement of the Isidis
shoreline. The eastern portion of Utopia also contains volcanic deposits from
Elysium that date to the Amazonian44, which could also contribute to loading.
While loading from the ocean itself is expected to produce some plate flexure, it
is not sufficient to explain the tilt of the Isidis basin (see Extended Data Fig. 5c),and water loading/unloading of Utopia is also insufficient to explain Isidis’ tilt41.
Therefore, deposition of material from a receding liquid, muddy or frozen ocean
may explain the tilt of the Isidis basin, even if some reverse tilting is caused bydeformation due to Tharsis.
Effect of elastic lithosphere thickness. The gravity and shape coefficients we use
to subtract Tharsis topography are based on an assumed Te =​ 58 km, the expectedvalue at the time of Tharsis loading14. However, the estimate of Te yields a 90%
confidence interval with a minimum and maximum of 26 km and 92 km,­respectively14. A thinner or thicker Te would alter the deformation due to Tharsis(and TPW) because the Love numbers used to compute Mars’ deformation wouldchange. To estimate the effect of Te =​ 26 km or 92 km on deformation due toTharsis, we recompute Tharsis’ gravity and shape coefficients following the methodof ref. 14. Using a fixed Tharsis centre location (258.6° E, 9.8° N), Matsuyama andManga14 compute the degree-2 gravity coefficients of Tharsis using a minimization
technique with four unconstrained model parameters (Te, non-dimensional Tharsis
load Q, palaeopole colatitude θR and palaeopole longitude ϕR), where the palaeopole corresponds to the axis of rotation when the fossil (remnant) bulge wasformed. This results in probability density functions for each unconstrained
parameter, with the weighted averages (expected values) used to compute the­gravity and shape coefficients. We redo this analysis, as described in section 5 ofref. 14, but with Te treated as a constrained parameter. This allows us to estimate
the expected values of Q, θR and ϕR for a given value of Te. We find that for
Te =​ 26 km, Q =​  3.95, θ R =​  17.9°, and ϕ R =​  259.1°. For Te =​ 92 km, Q =​  1.57,
θ R =​  14.2°, and ϕ R =​ 259.3°. Tharsis’ degree-2 gravity coefficients are recomputed
using these values. The degree-3 to -5 gravity coefficients of Tharsis are computedfrom minimization against the observed degree-3 to -5 gravity coefficients, andare therefore not dependent on Te. Shape coefficients for Tharsis are computed up
to degree 5 following section 7 of ref. 14. We compute the load Love numbers using
the ALMA code45, with a five-layer model as described in ref. 14.
We construct new best-fit deformation profiles for Te =​ 26 km and 92 km,
but with the corresponding Tharsis gravity and shape coefficients that we computed for each Te. The best-fit profiles for Te =​ 26 km and 92 km are comparedwith the ­nominal Te =​ 58 km profiles in Extended Data Fig. 3. All best-fit profiles
are ­relatively similar, showing that changes in Te do not have much effect on ourconclusions.
Effect of plate flexure. Although Perron et al.10 found that plate flexure due toloading of the ocean basin should not substantially affect the shoreline elevations,
their analysis was for Te =​ 200 km, whereas we use Te =​ 58 km. The ocean basin
resulting from our analysis also has less volume and a different shape, because
we subtract Tharsis topography, which has a negative component in much of the
Borealis basin. To compute plate flexure due to ocean loading, we expand thesurface density of the ocean load in spherical harmonics and compute the associated displacement using the method described in ref. 46. For the Arabia ocean,the ocean load is computed by subtracting the pre-Tharsis topography of Mars
from the best-fit Arabia ocean elevation (Z =​  −​2.3 km). The pre-Tharsis Martian
topography is computed by subtracting the deformation due to Tharsis and TPW,
equations (2) and (1), from Mars’ current topography (0.25° per pixel gridded
MOLA data21). For the ocean level corresponding to the Deuteronilus and Isidis
shorelines, only 17% of deformation due to Tharsis was subtracted from Mars’current topography, and the ocean elevation Z was set to −​3.68 km and −​3.95 km,
respectively. We use a Young’s modulus of 70 GPa, Poisson ratio of 0.25, and anassumed value of Te =​ 58 km. The loaded shoreline topography is compared to the© 2018 Macmillan Publishers Limited, part of Springer Nature. All rights reserved.

RESEARCH Letter
unloaded topography in Extended Data Fig. 5. We compute a maximum magnitude of deflection of 134 m, 84 m and 57 m, for the Arabia, Deuteronilus and Isidis
shorelines, respectively. The mean magnitude of deflection is 35 m for the Arabia
shoreline and 17 m for the Deuteronilus and Isidis shorelines. Deformation of theshorelines due to unloading of the ocean basin is negligible.
Data availability. The data that supports the findings of this study are available on
request from the corresponding author. Gravity and shape coefficients for Tharsis
are included in the Supplementary Information. Shoreline data should be requested
from the respective sources.
34. Harada, Y. Long-term polar motion on a quasi-fluid planetary body with anelastic lithosphere: semi-analytic solutions of the time-dependent equation.
Icarus 220, 449–465 (2012).
35. Chan, N.-H. et al. Time-dependent rotational stability of dynamic planets withelastic lithospheres. J. Geophys. Res. Planets 119, 169–188 (2014).
36. Hu, H., van der Wal, W. & Vermeersen, L. L. A. A full-Maxwell approach for large
angle polar wander of viscoelastic bodies. J. Geophys. Res. Planets 122,
2745–2764 (2017).
37. Moore, K. M., Chan, N.-H., Daradich, A. & Mitrovica, J. X. Time-dependent
rotational stability of dynamic planets with viscoelastic lithospheres. Icarus
289, 34–41 (2017).

38. Sjogren, W. L. Mars gravity: high-resolution results from Viking Orbiter 2.
Science 203, 1006–1010 (1979).
39. Zuber, M. T. et al. Internal structure and early thermal evolution of Mars from
Mars Global Surveyor topography and gravity. Science 287, 1788–1793(2000).
40. Searls, M. L., Banerdt, W. B. & Phillips, R. J. Utopia and Hellas basins,
Mars: twins separated at birth. J. Geophys. Res. Planets 111, E08005(2006).
41. McGowan, E. M. & McGill, G. E. Anomalous tilt of Isidis Planitia, Mars. Geophys.
Res. Lett. 33, L08S06 (2006).
42. Fairén, A. G. A cold and wet Mars. Icarus 208, 165–175 (2010).
43. Werner, S. C. & Tanaka, K. L. Redefinition of the crater-density and absolute-age
boundaries for the chronostratigraphic system of Mars. Icarus 215, 603–607(2011).
44. Tanaka, K. L., Robbins, S. J., Fortezzo, C. M., Skinner, J. A. & Hare, T. M. The
digital global geologic map of Mars: chronostratigraphic ages, topographic andcrater morphologic characteristics, and updated resurfacing history. Planet.
Space Sci. 95, 11–24 (2014).
45. Spada, G. ALMA, a Fortran program for computing the viscoelastic Love
numbers of a spherically symmetric planet. Comput. Geosci. 34, 667–687(2008).
46. Arkani-Hamed, J. The lunar mascons revisited. J. Geophys. Res. Planets 103,
3709–3739 (1998).

© 2018 Macmillan Publishers Limited, part of Springer Nature. All rights reserved.

Letter RESEARCH

Extended Data Figure 1 | Illustration of the feasibility of post-Tharsis
TPW depending on the location of Tharsis’ formation. a–c, Tharsis(orange zone) forms far from the palaeo-equator (a), causing large-scale
TPW as the planet reorients so Tharsis is at the equator (b). After thereorientation, the fossil bulge is far from the current equatorial bulge,
making the rotation pole sufficiently unstable to allow for ocean loading(blue zone) to cause subsequent TPW along an arc 90° from Tharsis(c, blue dashed line)10. d–f, Alternatively, Tharsis forms near the palaeoequator (d), causing limited (approximately 20°) TPW (e). The position of
the fossil bulge near the equator stabilizes the planet against subsequent
TPW caused by oceans and other surface loads (f, see supplementary
figure 1 of ref. 10). Formation of Tharsis near the equator is supported byrefs 14 and 15.

© 2018 Macmillan Publishers Limited, part of Springer Nature. All rights reserved.

RESEARCH Letter

Extended Data Figure 2 | Map of shoreline locations, MOLA
topography, and Tharsis deformation. Arabia (magenta) shoreline
data are from ref. 10 (data originally from ref. 7). Deuteronilus (white)

and Isidis (cyan) shoreline data and regional names are from ref. 18. The
contribution of Tharsis to Mars’ topography up to degree-5 (equation (2))
is displayed as 1-km dark grey contours (dashed contours are negative).

© 2018 Macmillan Publishers Limited, part of Springer Nature. All rights reserved.

Letter RESEARCH

Extended Data Figure 3 | Effect of elastic lithosphere thickness on
deformation due to Tharsis. a, Current Arabia shoreline topography
compared to displacement due to TPW and Tharsis deformation(equation (3)). Tharsis gravity and shape coefficients are computed
separately for Te =​ 26 km, 58 km and 92 km (see Methods), which eachyield a corresponding best-fit offset Z and error σrms. Dashed lines

show the best fit when the percentage of Tharsis topography added
after shoreline formation was allowed to vary by a factor C. Solid lines
assume 100% percent of Tharsis topography was emplaced after shoreline
formation (C =​  1). b, Deuteronilus shoreline topography compared to thebest-fit displacement due to Tharsis loading (equation (4)) for Te =​ 26 km,
58 km and 92 km.

© 2018 Macmillan Publishers Limited, part of Springer Nature. All rights reserved.

RESEARCH Letter

Extended Data Figure 4 | Comparison of Isidis shoreline topography
to shoreline deformation models. Current Isidis shoreline topography(elevation data from ref. 18) compared to the Perron et al.10 model
for Te =​ 200 km and our model of deformation due to partial Tharsis
emplacement (0.17∆T Tharsis −​3.95 km). The topography of the Isidis

shoreline can be explained by subsequent loading of the Utopia basin(see Methods). The starting point for the shoreline is (82.32° E, 7.36° N),near the southwest (SW) rim, with shoreline data proceeding clockwisethrough the northeast (NE) rim.

© 2018 Macmillan Publishers Limited, part of Springer Nature. All rights reserved.

Letter RESEARCH

Extended Data Figure 5 | Effect of plate flexure due to ocean loading on shoreline topography. Current shoreline elevations are plotted against
displaced elevations for the Arabia shoreline (a), the Deuteronilus shoreline (b) and the Isidis shoreline (c).

© 2018 Macmillan Publishers Limited, part of Springer Nature. All rights reserved.