Higher than present global mean sea level recorded by an Early Pliocene intertidal unit in Patagonia (Argentina)

Reconstructions of global mean sea level from earlier warm periods in Earth’s history can help constrain future projections of sea level rise. Here we report on the sedimentology and age of a geological unit in central Patagonia, Argentina, that we dated to the Early Pliocene (4.69–5.23 Ma, 2σ) with strontium isotope stratigraphy. The unit was interpreted as representative of an intertidal environment, and its elevation was measured with differential GPS at ca. 36 m above present-day sea level. Considering modern tidal ranges, it was possible to constrain paleo relative sea level within ±2.7 m (1σ). We use glacial isostatic adjustment models and estimates of vertical land movement to calculate that, when the Camarones intertidal sequence was deposited, global mean sea level was 28.4 ± 11.7 m (1σ) above present. This estimate matches those derived from analogous Early Pliocene sea level proxies in the Mediterranean Sea and South Africa. Evidence from these three locations indicates that Early Pliocene sea level may have exceeded 20m above its present level. Such high global mean sea level values imply an ice-free Greenland, a significant melting of West Antarctica, and a contribution of marine-based sectors of East Antarctica to global mean sea level. Global mean sea level was 28.4 ± 11.7 m higher than at present during the Early Pliocene, at atmospheric CO2 levels of no more than 450 ppm and temperatures of 2–3 ∘C above preindustrial levels, suggests a reconstruction from Patagonia.

T he survey, interpretation and dating of paleo relative sea sevel (RSL) indicators (such as fossil coral reefs or relic beach deposits 1 ) is paramount to constraining the maximum elevation reached by global mean sea level during periods of the Earth's history warmer than the pre-industrial. The elevation of paleo RSL indicators is the only direct proxy available to estimate global mean sea level in Earth's past. Once measured, observed paleo RSL indicators must be corrected for processes causing "Departures from Eustasy" 2 (such as tectonics, mantle dynamic topography, DT, and glacial isostatic adjustment, GIA 3,4 ) to obtain paleo global mean sea level (GMSL) estimates. These are in turn important to informing models of ice sheet melting under future warmer climates 5 .
A recent global database 6 shows that about 5000 RSL indicators were preserved since the Last Glacial Maximum (30 ka). Wellpreserved and dated RSL indicators are relatively rare for older time periods: another compilation of Pleistocene RSL indicators 7 reports more than 1000 Last Interglacial (MIS 5e, 125 ka) and only around 20 MIS 11 (400 ka) RSL indicators. Only a handful of sites exist that document sea level highstands beyond one million years ago 2,[8][9][10][11] . In general, robust RSL indicators predating 400 ka are rare to find because they are poorly preserved and are most often difficult to date with precision. In addition, relating them to GMSL is difficult since they are likely affected by significant postdepositional movements. This limits our ability to gauge the sensitivity of ice caps to warmer climate conditions, such as those that characterized Earth in the Pliocene.
Some of the oldest, precisely dated and measured RSL indicators were recently reported on the island of Mallorca (Balearic Islands, Spain), in a coastal cave called "Coves d'Artá". Here, six phreatic overgrowths on speleothems mark the paleo water/air interface within the cave 9 , and are therefore closely related to paleo RSL. The highest and oldest of these formations was measured at 31.8 ± 0.25 m above mean sea level, and yielded a U-Pb age of 4.29 ± 0.39 Ma (2σ) 9 . Taking into account GIA and possible long-term deformation due to tectonics or dynamic topography, it was estimated that global mean sea level at the time of deposition of this RSL indicator was 25.1 m above present, bounded by uncertainties represented by 16th-84th percentiles of 10.6-28.3 m 9 . For the same time period, a second study 10 reported a site in the Republic of South Africa (Northern Cape Province, site Cliff Point-ZCP Section2). Here, oyster shells living in a paleo subtidal to intertidal environment constrain paleo RSL at 35.1 ± 2.2 m (1σ). The oysters were dated to 4. 28-4.87 Ma (2σ range) with strontium isotope stratigraphy (SIS). While paleo global mean sea level estimates were not calculated at this site, based on the Mallorca benchmark the authors argue that this location was affected by relatively minor vertical land movements (possibly uplift) since 5 Ma.
While indirect paleo sea level estimates spanning the last 5.3 Ma are available from oxygen isotopes [12][13][14] , the two studies cited above are arguably the only ones reporting relatively precise and well-dated direct sea-level observations for the Early Pliocene, that is regarded as a past analog for future warmer climate 15 . At this time, CO 2 was between pre-industrial and modern levels, with possibly higher peaks to 450 ppm 15,16 . During Early Pliocene interglacials, average global temperatures were 2-3 ∘ C higher than pre-industrial values 15,17 . Pliocene climate was modulated by a ca. 40 kyr periodicity in glacial/interglacial cycles with highstands and lowstands that were characterized by sea-level oscillations as high as 13 ± 5m 18 . Ice models suggest that, during the warmest Pliocene interglacials, Greenland was ice-free 19 . Similarly, they suggest that the West Antarctic Ice sheet was likely subject to periodic collapses 20 , and might have contributed as much as 7 m 21 to GMSL. Ice models and field-based evidence 22 suggest that also the East Antarctic Ice Sheet might have been smaller than today, contributing another 3 m 21 to 13-16 m 23 to GMSL.
In this study, we report a foreshore (intertidal) sequence located in the town of Camarones, along the coast of central Patagonia, Argentina (Fig. 1). Combining field data, SIS ages, GIA and DT models we conclude that this deposit formed 4. 69-5.23 Ma ago (2σ range) when sea level was 28.4 ± 11.7 (1σ) higher than today. This estimate is broadly consistent with those derived from the Republic of South Africa and Spain. Together, these three studies present a coherent picture of global mean sea level during the Early Pliocene, that likely exceeded 20 m above modern sea level.
The Patagonia geographic region includes territories belonging to the states of Argentina and Chile. Geologically, Patagonia represents the southernmost tip of the South American plate (Fig. 1a). Along the Pacific coasts of Patagonia, the Nazca and the Antarctic plates are subducting below the Andes. Towards the south, the Scotia plate moves eastward and outlines Tierra del Fuego, at South America's southern tip 24 . To the East, the Patagonian Atlantic coast is a passive margin, tectonically characterized as an extensional stress field and bordered by a wide continental shelf. The central and eastern parts of this landmass are represented by the Andean foreland, formed by a Paleozoic-Mesozoic metamorphic basement overlapped by Tertiary continental and marine sedimentary rocks, dating back to the Paleocene. These are covered by Eocene-Oligocene pyroclastic rocks and Middle Miocene fluvial sediments. Marine sedimentary rocks corresponding to Tertiary transgressions are located east of the Andean foreland 25 . In the Middle Miocene, the Chile Triple Junction migrated northward, leading to the opening of an asthenospheric window below southern Patagonia 26 . This caused a switch from subsidence to uplift, and the Patagonia region underwent a moderate but continuous uplift 27 .
Along the coastlines of Central Patagonia, several levels of paleo shorelines above modern sea level were noted by Charles Darwin in his Beagle voyage 28 , and were the subject of more than 150 years of research (See Supplementary Note 1 and Supplementary  Table 1). Studies of Pleistocene coastal sequences in Central Patagonia include outcrops of Holocene 29,30 , Pleistocene [31][32][33] , and Pliocene-to-Miocene 34,35 age. Among the latter, Del Río et al. 35 dated Early Pliocene mollusks from marine deposits few hundreds of kilometers south of the study area described in this study.
The town of Camarones lies at the northern tip of the San Jorge Gulf,~1300 km south of Buenos Aires. Within a few kilometers of Camarones, several paleo-sea level indicators have been preserved, from the Holocene 36 to the Pleistocene 31 . Already in the late 1940s, the Italian geologist Feruglio 37 identified an elevated marine terrace along a roadcut carved on the main road leading into the town of Camarones that he tentatively attributed to the Pliocene. He called this terrace, the Camarones High Terrace (originally, in Spanish, Teraza Alta de Camarones 37 ). A recent study 31 confirmed the elevation of the Camarones High Terrace at ca. 40 m above sea level, at the lower bound of the "beach barries and terrace deposits between 40 and 110 m elevation" reported by the 1:250.000 geological chart of Camarones 38 .

Results
Pliocene sea level record at Camarones and GMSL estimates.
Radiometric ages, precise GPS elevations and stratigraphic descriptions of cross-sections surveyed along the Camarones High Terrace are the subject of this paper. Along this terrace, we surveyed and dated samples from two sites, separated by less than one kilometer. One is the Roadcut, already recognized and described by Feruglio 37 . We did not find reports of the second site (that we here call Caprock, Fig. 1b) in the existing literature, although it is possible that it was included in the geological description of the High Terrace by previous authors. At both sites, we recognized a geological facies representative of sedimentation in a foreshore environment (i.e., in the intertidal zone) that marks paleo RSL with high accuracy. All data described hereafter and in Supplementary Note 2 is available in spreadsheet form from Rovere et al. 39 .
Paleo RSL. In general, Roadcut and Caprock represent sedimentation during a transgressive event on top of a raised shore platform ( Supplementary Figs. 1-2). Among the units identified within the Roadcut (Fig. 2), one (Unit Cp, see inset in Fig. 2) is composed of well-cemented fine conglomerates with rounded pebbles and shells. In particular, the uppermost part of this unit contains a dense faunal assemblage in the form of a shellbed, where we recognized 15 different species of bivalves and 11 species of gastropods (Supplementary Table 2). The bivalve shells are mostly intact and sometimes with paired valves (articulated), but not in living position. This unit was interpreted as representative of a foreshore environment, i.e., the intertidal zone. The same unit has been identified at the Caprock section, at roughly the  same elevation. The elevation of Unit Cp was measured at two points at both Roadcut and Caprock (Table 1). From these measurements, we calculate that Unit Cp has an average elevation of 36.2 ± 0.9 m (1σ) above the GEOIDEAR16 geoid 40 , which is the best approximation for present sea level in Argentina. Using modern tidal values 36 , and assuming no post-depositional movement, we calculate that the two outcrops in the area of Camarones are indicative of a paleo RSL at 36.2 ± 2.7 m (1σ) above present (see "Methods" section for details).
Age. Three oyster shells from Roadcut and Caprock were analyzed by strontium isotope stratigraphy (SIS) relative dating techniques. Using sequential leaching to target the least altered inner carbonate of each shell, we obtained multiple SIS ages on three different shells (one from Caprock and two from Roadcut; see Sandstrom et al. 11 for a detailed description of the adopted methodology). The shells yielded an age range of 4.69-5.23 Ma (n = 6, 2σ S EM ).
Glacial isostatic adjustment. The Early Pliocene intertidal units surveyed at Camarones were subject to processes that caused their past and current elevation to depart from GMSL. These include glacial isostatic adjustment (GIA) and other vertical land motions (VLMs). We calculate GIA using 36 different Earth models. For this site, we calculate a GIA correction of −14.6 ± 3.2 m (1σ) (see "Methods" section for details). This value is subtracted from the observed paleo RSL and the uncertainty propagated. This correction is a combination of effects associated with (i) the ongoing response to the last deglaciation, and (ii) Antarctic ice sheet oscillations during the early Pliocene 2 . The former contribution is −9.5 ± 3 m (1σ), which means that the Argentinian coast today experiences sea level fall due to a combination of effects associated with postglacial rebound due to the melting of the glacial Patagonian ice sheet, as well as continental levering, ocean syphoning, and rotational effects. Once fully relaxed, sea level at Camarones will therefore be lower (and a paleo sea level indicator higher) by~9.5 m than it is today. The additional contribution of −5 m is associated with the adjustment to 40 kyr oscillations in the Antarctic ice sheet. The result is that, at Camarones, GIAcorrected paleo RSL is 50.8 ± 4.2 m (1σ).
Vertical land motions. The GIA-corrected RSL elevation reported above needs to be further corrected for VLMs, that can be either due to crustal tectonics, mantle dynamic topography 41,42 or deformation associated with sediment loading/unloading 43,44 . As briefly outlined in the previous sections, Camarones is located on a passive margin, likely subject to limited tectonic influence. Dynamic topography models suggest that, since MIS 5e (125 ka), the area of Camarones was subject to uplift, with rates increasing towards the South 3 . This is in line with observations of much higher Pliocene shorelines (70-170 m above sea level 35 ) at locations 300-500 km south of Camarones (Supplementary Note 1).
A long-term slight uplift trend is also predicted by the models of Flament et al. 45 and Müller et al. 46 . Predictions in these DT models average to 4.5 ± 2.2 m/Ma (Table 2). Accounting for the age of the deposit (including 1σ uncertainties), this leads to a downward correction of our global mean sea level inference by 22.4 ± 11.0 m (1σ). As is apparent from the variation of estimates for the dynamic topography rate, this correction remains quite uncertain and the true value can possibly be even outside of this range given that it is difficult to fully explore model uncertainties (see "Discussion" section).
Global mean sea level. Using the value of VLM reported above and propagating the uncertainties related to RSL, GIA, and VLM, we calculate that, at the time of deposition of the Caprock and Roadcut outcrops, GMSL was 28.4 ± 11.7 m (1σ). We remark that there are large unknowns associated with this value. First, as described above, dynamic topography remains a process that has high uncertainties that are generally not fully quantified. Second, it is possible that, as it is the case for the US Atlantic Coastal Plain 43 , flexural response to sediment loading or tectonic deformation (that are not considered here) could also contribute to further vertical land motions in this area.

Discussion
Our results show that the intertidal units at Camarones are of Early Pliocene age (4.69-5.23 Ma, 2σ S EM ). The sedimentological and stratigraphic characteristics of the deposits analyzed in this study lead to the conclusion that they formed during a sea level highstand, when GMSL was 28.4 ± 11.7 m (1σ) higher than present. We note that there are still large uncertainties on this GMSL  estimate, which derive mostly from vertical land motion corrections, stemming from the variability of published dynamic topography predictions 45,46 . Exploring and reducing these uncertainties requires improved mapping of the mantle structure beneath Patagonia from seismic tomography, a better understanding of how wave speeds map into density variations, and improved constraints on the rheology of the subsurface. Recent advances tackle these shortcomings and promise to reduce uncertainties in the estimate of vertical land motion 47,48 . Another strategy to investigate vertical land motions at Camarones would be to use the Pleistocene shorelines at the same site to extract a long-term uplift rate for the area. We argue that such approach would lead to similarly large error bars due to uncertainties related to GIA, Pleistocene global mean sea level and the implicit assumption that uplift rates can be linearly extrapolated over these time scales 49 . Despite the uncertainties related to VLMs, there is overlap between the calculated global mean sea levels for Camarones (28.4 ± 11.7 m, 1σ) and Coves d'Artá (Spain 9 , 25.1 m, with 16th-84th percentiles of 10.6-28.3 m, Fig. 3a,b). Correcting the proxy record at Cliffs Point (South Africa 10 ) with the same GIA models used for Camarones (Table 3), results in a paleo RSL of 44.7 ± 2.7 m (1σ) above present. The DT model predictions by Müller et al. 46 , which were also used for Camarones, indicate VLMs in the range of 4.6 ± 7.8 m/Ma (1σ). This results in an average global mean sea level estimate that aligns with those obtained from the other two sites, but bounded by very large uncertainties (23.4 ± 35.8 m, 1σ), (Fig. 3b). As already underlined by Hearty et al. 10 , improving uplift estimates for this region is paramount to enable the use of RSL data in GMSL calculations.
The average global mean sea level calculated from the geological facies reported in Argentina (this study), South Africa 10 and Spain 9 is well above modern sea level. Compared to published global mean sea level estimates that are based on ice sheet models and indirect sealevel proxies (Fig. 4), it is evident that field evidence is most consistent with the highstands obtained by scaling the Lisiecki and Raymo (2004) 50 benthic oxygen isotope stack (see "Methods" section for details). Our data is also consistent with some peaks predicted by the one-dimensional ice sheet model of Stap et al. 51 . Other ice sheet model based estimates [52][53][54] significantly under predict the observed Early Pliocene sea level records presented here. The almostcontinuous Gibraltar record 12 , derived from planktic δ 18 O coupled with a hydraulic model, largely over predicts sea level observed at both Argentina and Spain suggesting that, when the Camarones outcrop was deposited, the Earth was substantially ice-free. To align with this record, the three sites in this study would have to be characterized by marked subsidence, instead of uplift as indicated by almost all dynamic topography models we considered. Early Pliocene observations from Argentina only overlap with lowstands of the Gibraltar record, which would have left regressive imprints. This is at odds with the sedimentological characteristics of the Roacut, which represents a transgressive system rather than a regressive one.
While GMSL estimates from South Africa 10 are affected by large uncertainties, their average value together with the Argentinian sealevel proxies presented in this study and those obtained from Spain 9 , suggest that Early Pliocene GMSL might have exceeded 20 m above present-day levels. Reaching the average GMSL calculated for Camarones (28.4 m) would require an ice-free Greenland (GrIS, 7.4 m sea-level equivalent 55 ), significant melting of the West Antarctic Ice Sheet (WAIS, 3.3 m sea-level equivalent 56 ) and the almost complete melting of marine sectors of the East Antarctic Ice Sheet (EAIS, 19 m sea-level equivalent 57 ). Reaching the lower end calculated for Camarones (16.7 m, 1σ below the mean) would require complete melting of the GrIS and WAIS, and melting of about 1/3 of the marine-based sectors of the EAIS. This scenario would match almost exactly a complete GrIS melting, and a contribution from Antarctica in line with the one modeled by Golledge et al. 58 . These authors calculated that the contribution of Antarctica to GMSL during an Early Pliocene (4.23 Ma) interglacial was 8.5 m, sourced primarly from WAIS and the Wilkes subglacial basin of EAIS. Reaching the upper end calculated for Camarones (40.1 m, 1σ above the mean) would require significant contributions of not only marine-based but also land-based sectors of the EAIS in addition to melting of the GrIS and WAIS. We note that geological proxies suggest that a significant melting of land-based portions of EAIS was unlikely over the past 8 million years 59 , which makes this last scenario less likely.

Conclusion
The Early Pliocene world was characterized by global annual mean temperatures of 2-3 ∘ C higher than pre-industrial, and CO 2 levels between 280 and 450 ppm 15 . In face of these relatively small differences in temperature and CO 2 , the Earth's climate was substantially different than today 16 , and ice sheets were significantly smaller. Until recently, field evidence to support the answer to the question "How high was global mean sea level in the Early Pliocene?" was elusive. In this study, we show that independent paleo sea-level indicators of similar age on three continents result in broadly similar GMSL estimates. While affected by large uncertainties, stemming mostly from vertical  The significance of the Early Pliocene and its potential role as analog for present-day and near-future warming must be taken into account as the world prepares to meet the "Paris Agreement" 60 goals and limit global warming below the 1.5 ∘ C threshold 61 .

Methods
Elevation measurements and paleo RSL estimates. We measured elevations with a differential GPS system (Trimble ProXRT receiver and Trimble Tornado antenna) equipped to receive OmniSTAR HP real-time corrections. As per technical specifications by the service provider, these corrections allow to measure, in optimal conditions, the elevation of a point with an accuracy of 0.1-0.6 m (2σ), depending on the survey conditions. We remark that, while at the Caprock outcrop there is a free view of the sky, at the Roadcut satellite reception is hindered by the vertical cliff face. This could explain, in part, the discrepancy in the two points collected at this outcrop at relatively short distance from each other. Data were originally recorded in geographic WGS84 coordinates and in height above the ITRF2008 ellipsoid. For each GPS point, we calculated heights above Mean Sea Level (orthometric height) subtracting from the measured ITRF2008 ellipsoid height the GEOIDEAR16 geoid height 40 . These geoidal elevations are the best available approximation of mean sea level in this area. GEOIDEAR16 was estimated to have an overall accuracy of 10 cm (https://www.ign.gob.ar/ NuestrasActividades/Geodesia/Geoide-Ar16). The location and elevations of Unit Cp at Roadcut and Caprock are reported in Table 1.
From these elevations, we calculate that the average elevation (μE) is 36.2 m. To calculate the elevation error (σE), we use the following formula: where N is the total number of filtered positions measured by the GPS during the survey (439, sum of "Number of filtered positions" in Table 1), σE p is the elevation error for each single point, μE p is the Height above geoid of each single point and μE is the average elevation (36.2 m) ( Table 1). On average, we calculate that the elevation of Unit Cp is 36.2 ± 0.9 m (1σ). The Unit Cp at the Roadcut and Caprock sites has been interpreted as forming in the foreshore zone, i.e., in the intertidal zone. This means that its indicative meaning 62 spans from Mean Lower Low Water (MLLW) to Mean Higher High Water (MHHW). Based on predicted tidal data for the harbor of Camarones, Bini et al. 36 report that the maximum tidal range (MHHW to MLLW) in Camarones is 5 m. We use this value (5 m) as the indicative range (IR) for a foreshore deposit in our area, and the midpoint between MHHW and MLLW (0 m) as reference water level (RWL). Then, using the formulas described in Rovere et al. 1 , we calculate paleo RSL and its associated uncertainty as follows: Using the equations above, we calculate that paleo RSL associated with Unit Cp is 36.2 ± 2.7 m. We highlight that this value does not take into account the possibility that, 5 Ma ago, tidal ranges were different than present-day ones, due to different shelf bathymetry under higher sea levels 63 .
To calculate global mean sea level (GMSL) and associated uncertainties, we used the following formulas: where μGIA is the average of the GIA models (Table 3) and μVLM is calculated as the product of mean dynamic topography rate (Table 2) multiplied by the average age of the deposit.
where σGIA is the standard deviation of GIA models shown in Table 3 and σVLM is calculated as follows: where μAge and σAge are the average and 1σ age of the deposit, and μRate and σRate are the average and 1σ rates derived from published dynamic topography models ( Table 2).
Strontium isotope stratigraphy ages. To attribute an age to Unit Cp, we used the Strontium Isotope Stratigraphy (SIS) curve published by McArthur et al. 64 (LOWESS version 5). Sr isotope ratios from carbonates are susceptible to postdepositional alteration, therefore, any significant reworking of Sr isotopes needs to be detected and discarded. Information on shell preservation was determined using 87 Sr/ 86 Sr measurements on sequentially leached shell material (assuming smaller Sr isotope variations between leaches implies better preservation 65,66 ) alongside standard screening techniques 35,67 and elemental analysis 68,69 ). A preservation index between "1" (unaltered) and "3" (highly altered) was established for each sample based on these criteria (Supplementary Note 3, Supplementary Figs. 3-7, Supplementary Tables 3-4) with samples scoring above "2.0" excluded from results. The same screening criteria have recently been used by Hearty et al. 10 and are discussed in Sandstrom et al. 11 . The latter also gives an overview of the limits and implications of SIS analyses for Plio-Pleistocene marine samples.
We selected Ostreidae species for SIS chronological constraints, primarily because these shells precipitate original calcite mineral phases, making them more robust to diagenesis than aragonitic shells. Sample screening and chemical Fig. 4 Comparison between sea-level data discussed in this study and global mean sea level derived from ice models [51][52][53][54] and indirect sea level proxies 12,50 . The blue curve shows the GMSL prediction that is used in the GIA model and based on scaling the benthic oxygen isotope record by Lisiecki and Raymo (2005) 50 following the steps described in the methods. Age ranges for observations are 2σ, while elevation ranges are 1σ for Argentina and South Africa, and 16th-84th percentiles for Spain. Horizontal black lines and graphics on the right side of the graph show total sea level equivalent for icefree Greenland (GrIS, solid line 55  processing was carried out at Lamont Doherty Earth Observatory (LDEO), and all 87 Sr/ 86 Sr measurements were made using Thermal Ion Mass Spectrometry (TIMS) on an IsotopX Phoenix at SUNY Stonybrook University (SBU) or a Finnigan Triton Plus at Lamont Doherty Earth Observatory (LDEO).
We measured three oyster shells, one from the Caprock and two from the Roadcut unit. The Caprock oyster (ACC1-A) was sampled in three different locations, with inner leaches measured on two of those splits, returning SIS ages of 4.59 Ma (3.88-4.93 Ma) and 5.21 Ma (4.96-5.44 Ma) (Fig. 5). The third sampling location was only measured for full dissolution, with an average SIS age of 4.65 Ma  Table 3). Based on these screening criteria, we exclude sample ACR1-Ctop-C, which appeared to have been altered by low 87 Sr/ 86 Sr fluids (possibly of through leaching of surrounding volcanic material from the Complejo Marifil 38 ). The remaining inner leaches that passed screening were averaged by filament to obtain an age of 4.98 +0.245/−0.295 Ma (n = 6, 2σ S EM ). In the text, this age is reported as a 2σ range, i.e., 4.69-5.23 Ma.
Glacial isostatic adjustment. To account for changes in vertical displacement and gravity field caused by GIA we use a gravitationally self-consistent sea level model, that accounts for the migration of shorelines and feedback of Earth's rotation axis 70 . We compute both the contribution to GIA from the amount of residual deformation caused by the most recent Pleistocene glacial cycles and from ice age cycles during the Pliocene.
For the first contribution we use the results from Raymo et al. 2 , who calculated the residual deformation associated with the ice model ICE-5G 71 . This ice history is paired with a suite of 36 different earth models with varying lithospheric thickness (48 km, 71 km, and 96 km), upper and lower mantle viscosities (3 × 10 20 and 5 × 10 20 Pa s for the upper mantle, and 3 × 10 21 -30 × 10 21 for the lower mantle) to calculate a mean and standard deviation in residual deformation (Fig. 6).
For the second contribution we follow the approach described in Dumitru et al. 9 by estimating ice mass variability based on the benthic stack 50 . Following Miller et al. 72 we prescribe that 75% of the benthic δ 18 O variability is due to ice volume changes (the rest being due to temperature) and a further scaling of 0.11 o / oo /10 m to  convert δ 18 O seawater into ice volume changes. These conversions are highly uncertain 13,73 , which highlights the need to obtain local sea level based ice volume estimates. Nonetheless, this scaling was used because it yielded comparable ice volume estimates to the results of Dumitru et al. 9 . To construct an ice history following this ice volume curve we only assume changes in Antarctic ice volume given evidence that continent wide expansion of northern hemisphere ice sheets did only start around 3.3 Ma 74 . However, we acknowledge that an earlier intermittent Greenland ice sheet might have existed 75 . We compute glacial isostatic adjustment using this ice history and the same suite of 36 different earth models described above. We extract local predictions of relative sea level for Argentina, Mallorca, and South Africa. To calculate global mean sea level changes we integrate the amount of water in the ocean basins as a function of time. We next calculate how this quantity has changed relative to the initial state and divide it by the oceanic area calculated at each time.
Note that this setup to calculate the GIA correction deviates slightly from the one described in Dumitru et al. 9 in three small ways, (1) we only consider one GMSL history for the Pliocene rather than a range of histories, (2) we only consider variability in southern hemisphere ice sheets and (3) we calculated GMSL as described above rather than as changes in grounded ice volume.
The GIA corrections from both processes are combined. In a last step we consider the age range for each sea level indicator and average the GIA correction during warm periods, which we define as times that had higher than average sea level over this time period 9 . The mean and standard deviation that is obtained is shown in Table 3. We also show the GIA correction calculated by Dumitru et al. 9 and note that the difference in mean GIA estimates stems mostly from our different definition of global mean sea level. For the analysis in the main text we use the GIA correction described in Dumitru et al. 9 for the datapoint from Mallorca and not the one recalculated here.
Vertical land motions. VLMs were extracted from published Dynamic Topography models 45,46 . The values extracted are reported in Table 2. Flament et al. 45 focus on the surface expression of subduction dynamics in South America. Their results are based on forward advection modeling with different tectonic surface boundary conditions. The different cases are based on different timings of slab flattening. Müller et al. 46 have a global focus and combine back advection (initialized with a seismic tomography model) and forward advection with tectonic surface boundary conditions. Their different models are based on different surface plate reconstructions and different viscosity profiles.

Code availability
The python scripts used to produce panels b and c of Fig. 3 and the main panel of Fig. 4 are available from https://doi.org/10.5281/zenodo.3689426 78 (MIT license). The computer code used to do the sea-level (GIA) calculation, written in MATLAB, is available on GitHub (https://github.com/jaustermann/SLcode).