On the geophysical processes impacting palaeo-sea-level observations

Past sea-level change represents the large-scale state of global climate, reflecting the waxing and waning of global ice sheets and the corresponding effect on ocean volume. Recent developments in sampling and analytical methods enable us to more precisely reconstruct past sea-level changes using geological indicators dated by radiometric methods. However, ice-volume changes alone cannot wholly account for these observations of local, relative sea-level change because of various geophysical factors including glacio-hydro-isostatic adjustments (GIA). The mechanisms behind GIA cannot be ignored when reconstructing global ice volume, yet they remain poorly understood within the general sea-level community. In this paper, various geophysical factors affecting sea-level observations are discussed and the details and impacts of these processes on estimates of past ice volumes are introduced.


Introduction
Growth and decay of ice sheets influence global climate via changing ocean chemistry and atmospheric circulation (e.g., Broecker and Denton 1989;Yokoyama and Esat 2011;Matsumoto and Yokoyama 2013). An accurate understanding of the timing and magnitude of past sealevel change are, therefore, a critical component of paleoclimate studies. Mass exchanges between terrestrial ice sheet complexes and glaciers and the ocean basins, since the end of the last glacial maximum (LGM; ca. 20 ka; Clark et al. 2009;Yokoyama et al. 2018) changed the magnitude and distribution of the load applied to Earth's surface ( Fig. 1: LGM Ice: ICE-6G). These changes in load in turn induce slow deformation of the Earth with corresponding gravitational and rotational effects. The process by which this occurs is known as glacio-hydro-isostatic adjustment (most commonly contracted to glacio-isostatic adjustment or GIA), and complicates the quantification of the present-day changes in global mean sea level and projections of future global mean sea level (Clark et al. 2016).
GIA is the most ubiquitous geophysical process that impacts relative sea level (RSL) and accurate calculation of GIA effects is critical to reconstructing global mean sea level (GMSL) from observations of RSL change. Areas directly proximal to ice sheets that are highly affected by isostatic effects are referred to as "near-field". In the "far-field", the dominant component of sea level is the volumetric effect of the water added to (or removed from) the oceans, and this contribution is spatially invariant and, therefore, is known variously as eustatic, barystatic, ice-equivalent or global mean sea level (throughout this paper, we will use the latter term and the acronym GMSL). The "intermediate-field" is the region in which isostatic effects are not the dominant component of sea level, but are large enough to prevent the formation of a high stand at any point during the Holocene.
At present, there are significant uncertainties associated with GIA modelling in the near-and intermediate fields that require large, high-quality observational data sets with broad temporal and geographical coverage to resolve them. For example, much of our knowledge of changes in GMSL and ice sheet volume during the deglaciation, since the LGM has relied heavily on the contribution of material obtained by marine coring and drilling, but there still remain a number of important Yokoyama and Purcell Geosci. Lett. (2021) 8:13 yet-unresolved questions with significant ramifications for near future sea-level change. Of particular interest are the scale and geographical distribution of on-going GIA effects and past rates of mass loss from the various sectors of the Antarctic Ice Sheet (AIS) in response to paleoclimatic forcing.
On longer time-scales, such as during the Mid-Pleistocene, Pliocene and Miocene, sea-level changes at particular locations are heavily influenced by plate tectonic and mantle dynamic effects. As crustal plates drift over regions of vigorous mantle activity, their elevations may be changed by tens to hundreds of meters. This phenomenon may be much larger in amplitude than the potential change in GMSL. Validation of numerical models of mantle dynamic effects against observational data is only possible, where these processes are dominant over all other contributions to relative sea level (Mitrovica et al. 2020;Richards et al. 2020). In general, however, there are  (Peltier et al. 2015)] limited opportunities to distinguish the effects of mantle dynamic processes from other geophysical processes with sufficient accuracy to validate numerical models of effects for the majority of Earth's surface.
Yet, another complicating factor is the additional isostatic effect due to changes in surface loading caused by erosion and sedimentation. The magnitude of these load changes can be considerable, in some instances amounting to hundreds of meters and applied over broad wavelengths (Browning et al. 2006). These sedimentary loading effects must be considered when interpreting older paleo sea-level indicators or isostatic effects at points close to very active sediment transport systems (e.g., Ferrier et al. 2017;Dalca et al. 2013).
Below, we introduce the primary geophysical factors affecting inferences of GMSL from RSL observations.

Indirect measures
Microfossils, such as foraminifera, record past global ice volume in the chemical composition of their calcium carbonate shells (known as "tests"). The global ocean and atmosphere circulation systems are driven from the low latitudes by a surplus of solar energy and a deficit at high latitudes, resulting in large-scale moisture transport from the equator to the polar regions. Of the two most abundant stable isotopes of oxygen, 16 O and 18 O, evaporated seawater contains a slight (0.5-0.8%) overabundance of the lighter oxygen isotope, because at constant temperature, lighter isotopes move faster and escape the ocean-atmosphere barrier at a higher rate (Fig. 2). Moisture-laden clouds repeatedly vaporize to ascend to higher altitudes or shed some moisture through condensation and rain, as they move poleward. The heavier oxygen isotopes tend to condense preferentially as precipitation occurs. The net result is that water vapor in clouds becomes progressively more depleted in heavy oxygen isotopes, becoming approximately 3-4% enriched in light isotopes over Greenland and 4-5% enriched over Antarctica. When snow is formed, heavy isotopes preferentially condense, driving the net light isotope enrichment to slightly lower values (by approximately half a percent). The exact magnitude of this fractionation depends on the air temperature during snow formation and provides a unique record of atmospheric temperatures, or climate, over time (Dansgaard 1964;Uemura et al. 2005). The residual heavy-oxygen enrichment in the oceans is less pronounced, because the total amount of water in the ocean is many orders of magnitude larger, yet it is still measurable with high-precision analytical equipment. During the LGM, approximately 130 m of sea-level equivalent ice was stored on land, enriching the oceans in 18 O by about 0.1-0.15% (Fairbanks 1989;Adkins et al. 2002;Yokoyama et al. 2000).
Open ocean foraminiferal oxygen isotope measurements using deep-sea sediment cores can thus provide a continuous record of past climate and sea-level changes (e.g., Lisiecki and Raymo 2005). Since temperature changes and ice sheet volumes both impact oxygen isotope values, the isotope signals resulting from each Fig. 2 Cartoon showing the hydrological cycle of the Earth causing δ 18 O changes in precipitation at different latitudes. Because of these processes, polar ice sheets hold relatively lighter isotopes of oxygen in ice as much as 3-5% compared to the open ocean of these effects must be disentangled. Carefully considering the uncertainty effects discussed below the oxygen isotope record provides detailed records of GMSL change through time. Changes in the amplitude of ice sheet volumes during glacial-interglacial transitions increased substantially at around 1 Ma (1 million years before present) in deep-sea sediments (Fig. 3), this event is known as the "Mid-Pleistocene Transition (MPT)". In addition, around this time (~ 1.25 to ~ 0.7 Ma), the duration of glacial-interglacial cycles switched from 40,000 to 100,000 years (Fig. 3). The cause of this transition is still under debate, but several mechanisms have been put forward (e.g., Clark et al. 2009;Abe-Ouchi et al. 2013).
In the case of the Red Sea and Mediterranean Sea, substantial bodies of water are almost entirely enclosed by land and their sole connection with the broader ocean basin takes the form of narrow shallow channels (Babel-Mandeb and the Straits of Gibraltar, respectively, both of which are < 300 m in water depth). Physical and chemical conditions within such basins are sensitive to the scale and rate of water exchange with the broader ocean which in turn is sensitive to water depth within the connecting channels. Geochemical signatures from within the margin basins may then be used to reconstruct sea-level change within the channel. Using this technique, foraminiferal oxygen isotopes have been used Fig. 3 Sea-level record from the last 5.5 Million years. The long-term record is produced by an indirect method, namely δ 18 O obtained from deep sea sediments (Lisiecki and Raymo 2005). The sea level since the last interglacial is reconstructed based on direct methods, such as corals (b). The last 35 ka of detailed sea level reconstructions using data obtained from Northeastern and Northwestern Australia is depicted in (c). Rapid rises and falls of sea level can be identified Webster et al. 2018;Ishiwa et al. 2019) Yokoyama andPurcell Geosci. Lett. (2021) 8:13 to identify millennial to orbital scale sea-level variations in these regions including Japan Sea (e.g., Siddall et al. 2003;Yokoyama et al. 2007;Rohling et al. 2017). Indirect methods are useful for reconstructing sea-level change in the more distant past. However, on these scales model uncertainties associated with GIA and mantle dynamic effects are generally large. While methods relying on benthic oxygen isotope records from the open ocean are little impacted by these effects, basin sills are significantly influenced, so that basin isolation oxygen isotope analysis becomes increasingly uncertain over longer time-scales. On longer time scales, the assumption that the lithosphere is elastic becomes less robust. There is also greater scope for different time scales of viscoelastic relaxation in this case. While the number of model parameters required to accurately represent GIA increases as we move into the remote past, the observational data set becomes sparser and of lower quality. The geomorphological indicators associated with ice sheet geometry and sea-level change through glacial cycles are very likely to be erased by subsequent glaciation and inundation episodes (e.g., Yamane et al. 2015). The ice-and waterloading histories associated with previous glacial cycles, therefore, become increasingly uncertain, as more glacial cycles are considered.
Similarly, observational constraints on mantle dynamic processes become looser and the geographical and temporal patterns of mantle processes become less welldetermined over longer time-scales (e.g., > 100,000 years).
The assumption that dynamic topographic processes are linear or even monotonic becomes less defensible on these time scales. There are some individual regions in which observational controls have been preserved and a reconstruction of dynamic topography is possible in these instances, but this is not generally the case (e.g., Mitrovica et al. 2020). Thus, GMSL reconstructions are subject to large uncertainties in both amplitude (~ 10-30 m) and timing (~ several thousand years) because of the limitations inherent in the radiometric dating techniques applied (radiocarbon is limited to the last ca. 50,000 years) and uncertainties associated with the oxygen isotope signals recorded by foraminifera due to seawater temperature fluctuations (Grant et al. 2014;Raymo et al. 2018;Miller et al. 2020).

Direct measures
Geologists have used a broad array of fauna and flora as direct sea-level indicators or sea-level index points (SLIPs), particularly corals, benthic foraminifers, ostracods, diatoms and mangroves from coastal regions, and sedimentary/stratigraphic features (marine terraces, ooids, and so forth) (Yokoyama et al. 2006;Tam and Yokoyama 2021). Microfossil-based reconstruction of sea level utilizes various taxa's affinity for particular salinities. Thus, horizons containing marine-freshwater transitions can be identified through microfossil assemblage analysis. For the Holocene, diatoms (whose frustules are silicabased) are particularly useful (see Fig. 4a), because they are often well preserved and inhabit a relatively narrow salinity range, particularly at the boundary between fresh and brackish water . As such, Holocene sea-level information has been extracted from diatoms at decimeter spatial scales (e.g., Yokoyama et al. 1996Yokoyama et al. , 2012. Except for uplifted regions, however, pre-Holocene diatom frustule preservation is often relatively poor, particularly in the tropics, because seawater tends to be undersaturated with respect to silica. Therefore, analysis of foraminifera and ostracod (Fig. 4) assemblages provide better information prior to the Holocene, because they are composed of CaCO 3 (calcite) (e.g., Scott et al. 1996;Murray 1991;DeDeckker 1988;Cann et al. 2006;Cann et al. 2006, Murray-Wallace et al. 2021. Micropaleontological methods are often complementary and may produce precise results when used in unison to study sealevel history. In tropical regions, reef-building corals are one of the most reliable direct measures of past sea level (e.g., Camoin et al. 2012;Yokoyama and Esat 2015;Woodroffe and Webster 2014). Their habitats are close to the sea level because of the symbiotic relationship with the algae zooxanthellae requires their habitat to be within the photic zone. Coral skeletons are composed of CaCO 3 (aragonite), so they can be radiometrically age-dated by both radiocarbon and uranium series methods (Yokoyama et al. 2001a, b;Yamane et al. 2019). Both tectonically uplifted reefs as well as drowned reefs have revealed millennial-(e.g., Yokoyama et al. 2001a, b) to decadal-scale sea-level variations (e.g., Sieh et al. 1999;Natawidjaja et al. 2006;Woodroffe et al. 2012;Hallmann et al. 2018). Mortality of corals at the top of the colony, at sea level, is often higher than on the periphery, which may create an atoll-like structure called a "microatoll" (Fig. 5). Microatolls have been reported to capture decimeter-scale sealevel fluctuations (Sieh et al. 1999;Natawidjaja et al. 2006;Hallmann et al. 2018;Woodroffe et al. 2012;Dechnik et al. 2019).
It should be noted that the formation of coral reefs, established marine ecosystems, erosion notches, rock platforms and similar sea-level indicators require that sea level at that point be relatively stable for an extended period and that conditions be favourable for their formation and subsequent preservation. As a result, sea-level indicators may not have an opportunity to form if sea-level changes rapidly (e.g., Murray-Wallace and Woodroffe 2014). Where sea-level indicators do form, their survival through to the present depends on appropriate environmental conditions prevailing. As a result, the sea-level record in many localities is incomplete. For regions in which the surviving SLIPs are widely separated in time, it is not possible to reliably reconstruct sea-level changes for the periods between the SLIPs and care must be exercised not to arbitrarily exclude the possibility of more complex sealevel histories in these areas than is suggested by a simple interpolation of the observational data.
In the far field, GIA-loading effects are small and the uncertainties associated with GIA effects is reduced to less than 10%, so that the correction of paleo-sealevel indicators to extract GMSL is readily achieved. Analysis of far-field sea-level records has permitted detailed reconstructions of global GMSL (e.g., Yokoyama and Esat 2011; Lambeck et al. 2002); however, data constraining LGM sea levels remain sparse. What data there is suggests GMSL was quite dynamic in this period with GMSL falling and rising by 20 m in the course of 4 ka (Yokoyama et al. 2001a(Yokoyama et al. , b, 2018Ishiwa et al. 2019). More observations are needed to better constrain the GMSL record beyond 18 kaBP.

Glacio-hydro-isostatic adjustment: relative and global mean sea level
The sea level observed from any particular geological archive, as compared to present-day sea level, includes the sum of all the imprinting geophysical factors [GIA, sediment compaction/erosion/deposition, tectonic effects, dynamic sea level (deviation from the geiod); (Cazenave and Emy 2011), wave climate, etc.] and is known as relative sea level (RSL). In this section, we describe the traditional GIA factors affecting observed RSL and derive the equations for GMSL. The change in RSL (Δζ rsl ) is defined by the equation: where ζ rsl (t 0 , φ) is present-day sea level at locality φ and ζ rsl (t, φ) is the elevation of the fossil shoreline. Farrell and Clark (1976) developed the theory for GIAinduced sea-level change on a self-gravitating earth in response to global mass redistribution due to growth and decay of former ice sheets. If we assume that the melt water from ice sheets is uniformly distributed over the oceans at time t, then a first approximation of sea-level change is given by where ρ ice and ρ water are the densities of ice and water, respectively (917 and 1027 kg/m 3 : Lambeck et al. 2014). ΔV ice is the global, land-based ice volume difference between time t and the present. A ocean is the surface area of the ocean at time t. Δζ gmsl is the change in GMSL, one of the most significant parameters in modelling paleosea-level change. Δζ gmsl provides only a zeroth-order approximation of sea level.
When ice sheets melt, the corresponding surface load is reduced and the lithosphere beneath the former ice sheet rebounds as previously displaced mantle material returns beneath the formerly depressed lithosphere. At the same time, meltwater loads the oceanic lithosphere and induces further mantle flow from oceanic lithosphere to beneath the continental lithosphere. Because of the redistribution of surface load and resulting displacement of underlying mantle material, Earth's gravitational field and rotational moment of inertia are influenced. As a consequence (and as discussed in Whitehouse 2018), the total relative sea-level change will be a combination of: (a) changing ocean volumes, (b) crustal displacement due to the changing load, (c) change in the shape of gravitational potential surfaces due to internal and external mass redistribution, and (d) changes in Earth's axis of rotation (Fig. 6).
Component (b) describes the deformation induced by changes in surface load. On short time scales, the deformation consists of a small amplitude near-instantaneous displacement that may be treated as essentially elastic. This is followed by a longer term relaxation process of much larger amplitude. On these longer time scales, the deformation is well-approximated as viscous a b Fig. 6 Figure showing schematics of various geophysical processes related to sea-level changes. a Crustal loading underneath of ice sheets in the "near-field" or center of ice sheet increases and mantle material flows beyond the ice sheet to form a peripheral bulge in the "intermediate-field" during ice sheet advances. Sites far from the ice sheet regions, where direct loading effects due to the change in ice sheet geometry is negligible are called "far-field". For these locations, changes in the geometry and magnitude of water load influence the change in RSL. b The opposite is true during the waning phase of global ice sheets. Crustal uplift and peripheral bulge collapse occur in the near-and intermediate-fields due to mantle flow. In the far-field, a sea-level high-stand results from ocean siphoning due to the combination of bulge collapse and mantle flow from beneath the newly-loaded ocean basins. Earth's rotational axis shifts because of mass exchanges between ice sheets and oceans and the effect of this process on Earth's moment of inertia which in turn alters the centrifugal force applied to the water within the ocean basins. Reduction or increase in ice sheet mass produces a corresponding fall or rise in near-field sea level due to the change in gravitational attraction of the ice sheet on the water in the ocean basin flow. This suggests that Earth's response to changes in surface loading is effectively viscoelastic on time scales of thousands to hundreds of thousands of years. Sophisticated theoretical treatments of radially symmetric Maxwell or Burger solids were developed by Peltier (1974), Yuen and Peltier (1982), and Wu and Peltier (1982). Recent advances in computational resources have facilitated investigations of inhomogeneous mantles (using 3-dimensional code, e.g., van der Wal et al. 2015) and non-linear rheologies (e.g., Wu and Wang 2008). While more sophisticated rheological models such as these are likely to more closely approximate the behavior of the mantle, it can be formally demonstrated that there is insufficient resolution in the sea-level record to constrain more than the coarse broad-scale character of the mantle's response characteristics (Paulson and Richards 2009). As a result, GIA inversions based on these more sophisticated mantle rheologies are unlikely to yield unique solutions.
The combined effects of (b), (c) and (d) are the isostatic contribution to RSL change. Because of Earth's viscous nature (Fig. 7), crustal displacements persist long after the melting of the ice sheets. For example, even though most of the last-glacial ice sheets decayed by about 7000 years ago during the last deglaciation, relative sealevel change continues to the present, and this effect in turn impacts modern sea-level change (Fig. 8).
In view of these various contributions, the sea-level change at time t at a site φ can be written schematically as (e.g., Farrell andClark 1976, Nakada andLambeck 1987;Mitrovica and Peltier 1991;Milne and Mitrovica 1998;Lambeck et al. 2003) where Δζ rsl is the change in relative sea level, Δζ iso is the total isostatic effect associated with the glacial rebound process and Δζ geo is any additional geophysical contribution.
(3)  Tushingham and Peltier (1991), c Mitrovica and Peltier (1991), d and e models for Europe from , f model for the Australian coast from , g Nakada and Lambeck (1989), h Ricard et al. (1989) These two terms include the gravitational, rotational and surface elevation effects caused by deformation under the ice and water loads, respectively. Because the magnitude of the water-load depends on the change in sea level which in turn depends on Δζ water (t, φ) the sea-level equation derived by Farrell and Clark (1976) is implicit and requires multiple iterations to achieve a convergent solution. The total isostatic component must also incorporate the effect of variations in the gravitational attraction between ice-sheets and ocean water and selfattraction of the water contained in the ocean basins. It must additionally include the corresponding influence of surface mass redistribution on Earth's rotation. These effects are incorporated into what we call the rigid Earth component, Δζ rig (e.g., Nakada and Lambeck 1987;Johnston 1993). This term is a function of both ice sheet and ocean geometry through time and is particularly significant in the near-field. With these definitions the isostatic component may be re-written: In practice, these three components are inter-dependent and the separation implied by this adopted notation is not formally rigorous, we have adopted this form for the sea-level components as a matter of convenience when discussing their relative significance.
With this expansion of the isostatic term Eq. 3 becomes A more complete discussion of these terms and their formal relationships has been presented by previous papers (e.g., Nakada and Lambeck 1987;Johnston 1993;Mitrovica and Milne 2003;Lambeck et al. 2003;Kendall et al. 2005).

Gravitational and rotational effects
In addition to deriving a self-consistent integral equation describing sea-level change, Farrell and Clark (1976) also demonstrated that when GMSL change is rapid (in comparison with deformation processes), the dominant near-field component is the change in direct gravitation exerted by the ice sheets. As the ice sheets lose mass their gravitational influence on ocean water decreases and in the near field sea-level falls (Figs. 6 and 9). There is a similar effect due to self-gravitation of water in the ocean basins which acts to increase ocean water depth toward the center of the oceans.
While the change in the volume of water contained in the ocean basins remains the dominant component of far-field sea-level change on longer time scales, there are second-order volumetric effects. In particular, ocean area may be impacted by shoreline migration and ocean volume may change as a result of ocean floor deformation and gravitational and rotational effects on ocean water depth.
The configuration and magnitude of the water load are sensitive to shoreline position, which in some cases may migrate on the order of thousands of kilometers (e.g., Lambeck et al. 2003). Grounding line position of the ice sheet further impacts the geometry of not only the water load, but also the ice load (Milne 1998).
As an ice sheet accumulates, the pressure applied to Earth's surface induces mass flow in the under-lying mantle from beneath the ice sheet to outside the ice sheet. The resulting crustal uplift combines with associated (4) gravitational effects to lower sea level in a broad annulus surrounding the ice sheet. In GIA literature, this region is called the peripheral bulge. In many cases, an ice sheet's peripheral bulge is submerged beneath the ocean. After the ice sheet retreats, crustal subsidence of a submerged peripheral bulge increases ocean volume (Johnston 1993) through the process known as ocean siphoning (Mitrovica and Milne 2002;Milne et al. 2009; see also Figs. 6 and 8). This effect is the dominant mechanism driving the widely observed fall in sea level since the mid Holocene (e.g., Hallman et al. 2018;Yokoyama et al. 2019a, b;Fukuyo et al. 2020).
At any point on Earth's surface, the gravitational acceleration experienced by a mass is partially offset by centripetal force due to Earth's rotation. Consequently, changes in the orientation Earth's rotational axis or its speed of rotation will alter the centripetal component of acceleration experienced by the water column and thus influence sea level (Mitrovica et al. , 2011. On short time scales, the dominant rotational effect is due to the redistribution of surface mass between ice sheets and oceans. On longer time scales, Earth's rotation is further influenced by load-induced mass redistribution within The value plotted in each panel is the dimensionless quantity (�ζ rsl (φ) − �ζ gmsl (φ))/�ζ gmsl (φ) , where Δζ gmsl (φ) is global mean sea level added to the oceans over a short (effectively instantaneous) time interval. Contours are in increments of 0.25 for positive and 1.0 for negative values. The far field signal is uniquely characteristic of the source region the mantle. As a result of the intricate feedbacks between the geometry of the water load and sea-level change, solving the sea-level equation requires multiple iterations to accurately determine water load geometry (e.g., Kendall 2005).
Relative sea-level indicators from intermediate locations constrain both the collapse of the peripheral bulge and the direct gravitation signal of the ice sheets (Fig. 8). These sites are critical in partitioning the total GMSL signal between the various ice sheets (e.g., Nakada et al. 2018); however, the rise in sea level since the LGM has inundated most sea-level indicators from this critical region. In the context of modern sea-level change, accurate modelling of on-going GIA effects plays an important role in the attribution of recent sea-level changes between individual ice sheets and glaciers.
Over short time scales, the deformation signal is relatively small. This allows the regional distribution of glacier and ice sheet mass loss to be inferred from patterns of global sea-level change through a process called "sealevel fingerprinting" (as illustrated in Fig. 9). This possibility was quickly pursued by Clark and Lingle (1977), and as computational resources improved, more detailed efforts incorporating the rotational component became feasible (e.g., Mitrovica et al. 2011). Exploiting this methodology, Hay et al. (2013) developed a Kalmann filter processing technique to extract time series of ice sheet contributions of sea-level change from tide gauge data sets. More recently, the fingerprinting methodology has been extended to also include the gravitational signal of the surface mass redistributions observed by satellite gravimetry (e.g., Hsu and Velicogna 2017).
"Fingerprinting" requires a unique signal that is dependent only on the approximate location of the meltwater sources involved. For that to occur, the ice sheets concerned must discharge their mass into the oceans over a relatively short time scale. If ice sheet mass loss takes place over a time scale on which significant viscoelastic deformation can occur then an inherent ambiguity will exist between rheological model and ice sheet parameters. On these longer time periods, standard inversion methods may be used to constrain the various model parameters, but this cannot be considered to be "fingerprinting".

Limitations of GIA modelling
Because RSL measures water depth relative to the underlying crust, the GIA effects described in the previous two sections must be properly constrained to relate sea-level indicators to contemporary GMSL. The precise character of GIA-induced sea-level change depends crucially on ice sheet history and earth response parameters. However, observational constraints on both these parameters are subject to significant uncertainties.
Geomorphological indicators of ice sheet extent are generally restricted to the last glacial cycle (in most cases, they are in fact restricted to the period since the most recent deglaciation), and are generally only preserved in regions and epochs, where the ice sheet is stable and only rarely provide direct constraints on ice thickness (e.g., Clark et al. 2009;Yamane et al. 2015). The geomorphological data available to constrain the ice sheet geometry cannot simply be assumed to be complete (e.g., Fabel et al. 2006;Heyman et al. 2011;Briner et al. 2014).
It is necessary, therefore, to supplement the glaciological indicators with near-field sea-level observations which are strongly dependent on the magnitude and timing of the change in ice load and the relaxation characteristics of the underlying mantle. However, sea-level indicators only occur in coastal regions and under particular environmental conditions (e.g., Mackintosh et al. 2011;White et al. 2011a;Yamane et al. 2011;Takano et al. 2012). Where they do form they may be subsequently eradicated or obscured by subsequent ice sheet re-advance, inundation or erosional effects (Heyman et al. 2011;White et al. 2011b). In practical terms, near-field sealevel indicators most commonly post-date the local glacial maximum (Behrens et al. 2019;Sproson et al. 2021). The accuracy of a GIA-based ice sheet reconstruction, therefore, decreases as one moves further back in time beyond the period for which near-field RSL observations are available. As a general rule, if near-field indicators first become numerous and wide-spread at time X ka BP, then the resulting ice sheet reconstruction cannot be treated as robust for any epoch earlier than X + 2 ka BP. As an example, the North American RSL record begins to become more substantial at approximately 16 ka BP and thus, the character of the North American ice sheet complex for the period before 18 ka BP is extremely uncertain (Lambeck et al. , 2017. This limitation is unlikely to be remedied as the presence of the ice sheet inherently suppresses SLIP formation. The rheological properties of the Earth at depth are also relatively poorly understood (e.g., Jain et al. 2018). Viscosity is dependent on composition, water content and temperature (Karato et al. 1986). In addition, there is evidence that the constitutive equation of the mantle may be non-linear in which case material viscosity is also stress-and time-dependent (Wu and Wang 2008). However, sea-level observations provide only a limited ability to resolve detailed structure within Earth's mantle, especially given the uncertainties associated with the loading history (Paulson and Richards 2009).
In early theoretical treatments of GIA (e.g., Peltier 1974), the Earth was modelled as being a radially stratified, spherically symmetric visco-elastic solid. Regional studies using this theory obtained different rheological parameters for different localities, suggesting that lateral variation might be significant (Lambeck et al. 1998(Lambeck et al. , 2003(Lambeck et al. , 2010. There have been significant, and rapidly accelerating, efforts to develop fully threedimensional rheological models for use in GIA modelling (e.g., Milne et al. 2018;Geruo et al. 2013) and to quantify the sensitivity of sea level to rheological and ice model parameters (e.g., Al-Attar and Tromp 2014; Crawford et al. 2018;Paulson and Richards 2009;Mitrovica et al. 2018). While three-dimensional rheological models can be readily employed in forward modelling of GIA effects, and are being used with increased frequency, inversion of observations to constrain three-dimensional rheology is problematic. The simple number of additional parameters introduced in a fully three-dimensional rheology is a severe numerical and theoretical burden, while as discussed above, the resolution of the observational data set is not sufficiently high to allow reconstruction of fine structure within the mantle (Nakada et al. 2018).
Treating rheology as an unknown complicates GIAbased ice sheet inversions, most immediately by enlarging the parameter space to be searched, but also by introducing possible trade-offs between ice sheet and rheological parameters. A large load applied to a strong mantle or a small load on a weak mantle could both be consistent with an observed uplift signal. Similarly, a small load applied for a long time or a large load applied briefly will also be difficult to discriminate between. There is an additional complexity in the case, where the SLIP record is restricted to only the later part of the relaxation curve (as is often the case). There is, therefore, ambiguity between ice and earth model parameters (Whitehouse 2018) and an inherent uncertainty associated with GIA model results that has not often been made explicit (e.g., Lambeck et al. 2017).
Assuming that rheology is known (e.g., Peltier et al. 2015) requires that the ice model adjust to take up any observational misfit resulting from the inaccuracy of the assumed Earth model, introducing a bias into ice model outputs. At any point in time, ice sheet geometry has a great many degrees of freedom in spatial extent and ice thickness, which are only loosely constrained. With these variables it is, in general, possible to find an ice model that broadly matches the available observational data sets for a given earth model. However, the accuracy of the resulting ice model is limited by that of the underlying rheological model, which in many cases is untested.
Conversely, many studies present an ice model derived from glaciological or climatological constraints and search for rheological parameters that, when combined with the ice sheet history, match the observations (e.g., Whitehouse et al. 2012). The accuracy of the underlying glaciological models remains open to question, with different studies producing very different ice sheet reconstructions (e.g., Philipon et al. 2006;Simon et al. 2010;Ivins et al. 2013;Boer et al. 2017). If cold-based conditions apply, then it may happen that no geomorphological indicators are produced (Sugden 1978). In that case, there is a strict limit on how accurate a glaciological inversion based on geomorphological indicators can be. Assuming a given ice model is accurate in turn biases the rheology parameter search to compensate for any inaccuracy in the ice model. As long as the proposed ice model is approximately correct, a series of rheological models can be found which will broadly match the observations. The more complex the rheological model, the closer the agreement with the observations can be made.
Effective separation of ice load and rheological parameters requires a dense spatial and temporal distribution of accurate sea-level indicators to constrain the wavelength of sea-level trends and their change through time.
The poorer the quality of the sea-level indicators or the sparser their distribution, the less robust the inferred ice sheet inversion will be, and the broader the range of permissible rheological models with which it may be combined. Sea-level records from the ice sheet margin during its retreat phase directly constrain the timing of the retreat and the character of the post-glacial relaxation (Yokoyama et al. 2016a, b;Valletta et al. 2018;White et al. 2019;Behrens et al. 2019). These indicators are particularly valuable in the reconstruction of the ice sheet.
Despite the implicit uncertainties in GIA-based ice sheet reconstructions, there has for many years been a tendency in at least some of the GIA modelling literature to present results for a single ice model (assumed to be correct) and a single rheology model (also assumed to be correct) with no discussion of the uncertainties associated with either. This lack of transparency has led many reading this literature to overestimate just how reliable the resulting model outputs are. While there have been some recent efforts to analyze and quantify the uncertainties associated with GIA models (e.g., Stocchi et al. 2018;Dendy et al. 2017;Rohling et al. 2017;Kuchar et al. 2020), the theory related to GIA model uncertainties is still maturing. At present, there is no consistently applied methodology for exploring the uncertainties associated with GIA models. It is, however, necessary to consider that solutions for ice sheet geometries and rheology obtained from GIA inversions are unlikely to be unique. Any combination of a single ice model with a single rheological model should not be accepted as authoritative (e.g., IPCC 2007). A range of variations in ice sheet geometry and mantle response should be considered when presenting GIA model results.
The frequency of collaborations between modelers and observational scientists is steadily increasing and producing increasingly detailed understanding within the sealevel community of the relative strengths and limitations of the component disciplines (e.g., Brenderyen et al. 2020;Bracegirdle et al. 2019). Collaborations of this sort play an important in the members of all communities acquiring a deeper understanding of the complexities involved in modelling, observing, projecting and interpreting sealevel change.

Antarctic ice volume since the LGM
On the time scale of the last deglaciation, sedimentary and dynamic topography effects are small. As discussed above sea-level indicators from tectonically stable farfield regions need only be corrected for GIA to extract changes in GMSL, and for these sites, the errors in the GIA calculation can be demonstrated to be small. In the near-field, however, the uncertainty in ice load and rheology, and the trade-off between them are significant (Nakada et al. 2016;Lambeck et al. 2014). In these circumstances, a large, accurate data set is required to achieve a robust inversion that correctly separates ice load geometry and rheology. Fortunately, such data sets exist for the northern hemisphere ice sheets, but Antarctica's data coverage is very poor. In particular, the inundation of the Antarctic continental shelf has made it extremely difficult to determine the location of the ice sheet margin during the deglaciation process and the post-glacial relaxation that has occurred at these points (Prothro et al. 2020;Bart et al. 2018).
Submerged near-field RSL indicators provide valuable constraints on the extent of the ice sheets over the continental shelves and are likely to be particularly significant for the Antarctic Ice Sheet (AIS), since its history is especially poorly known. A number of attempts have been made to reconstruct histories for the AIS using glaciological constraints (e.g., Whitehouse et al. 2012;Philippon et al. 2006), but these models produce widely divergent estimates of AIS volume ranging from as much as 20 m to as little as 10 m of GMSL. However, sea-level records from far-field sites (e.g., Lambeck et al. 2014;Yokoyama et al. 2018) and reconstructions of the northern hemisphere ice sheets suggest that the Antarctic ice sheet must have contained at least 18 m of GMSL at its maximum extent (Peltier et al. 2015) and very possibly 27 m or more (e.g., Lambeck et al. 2017). Though large, this latter value is consistent with other geophysical data (Nakada et al. 2018). A larger collection of near-field RSL and ice coverage data from the Antarctic continental shelf is required to better discriminate between these various hypotheses (e.g., Yokoyama et al. 2016a, b;Behrens et al. 2019;Simms et al. 2019;Prothro et al. 2020).
The recent suggestion of Gowan et al. (2021) that there is no "missing ice" problem to be resolved seems to be based on a misinterpretation of the LGM sea level constraints employed in that study. In particular, the Bonaparte Gulf data DeDeckker and Yokoyama 2009;Ishiwa et al. 2016Ishiwa et al. , 2019 is characterised by persistent estuarine conditions that require sea level at this location to have fallen at or very close to the level of the basin sill. Similarly, subaerial exposure is evident for some the Great Barrier Reef observations of Yokoyama et al. (2018), inconsistent with their interpretation as corresponding to mean sea level as is adopted by Gowan et al. (2021).
The uncertainty in estimates of the AIS geometry through its most recent deglaciation is of immediate concern in the context of modern sea-level change. The stability of the various sectors of the AIS to climate forcing is a significant unknown in projections of future sea-level change. In ddition, the magnitude of the Antarctic GIA signal directly impacts estimates of modern AIS mass loss (e.g., Golledge et al. 2019). Sea-level records from the Antarctic continental shelves have the potential to greatly assist modelers attempting to distinguish between various retreat scenarios, which cannot presently be differentiated, because the observational data are too sparse. In this context, marine drilling programs are likely to be particularly valuable (Yokoyama et al. 2019b).

Last interglacial sea-level uncertainty
The Last Interglacial (LIg; ca. 120 ka) occurred during Marine Isotope Stage (MIS) 5, and sea-level records from this epoch have been the focus of active investigation (e.g., Dutton et al. 2015), because the LIg has been traditionally considered an analog for the present interglacial, has better geochronological control than other potential analogs such as MIS 11, and has the potential to illustrate the impact of a two-degree increase in temperature on GMSL. However, the interpretation of LIg sea-level indicators and GMSL is greatly complicated by a number of factors. Nakada and Lambeck (1987) and Lambeck et al. (2012) showed that the elevation of LIg sea-level records relative to modern sea level depends on the point during the LIg the SLIP was formed. The duration of past interglacials and the deglaciation phase preceding each varied with differing combinations of orbital parameters. The mantle relaxation that occurs through each interglacial thus varies accordingly and SLIPs will occur at different elevations relative to modern sea level based on their time of formation. In addition, Lambeck et al. (2010), Dendy et al. (2017) and Rohling et al. (2017) showed that the configuration of the major ice sheets from MIS 6 to MIS 4 also impacts the modern elevation of LIg sea-level indicators. Kopp et al. (2009), Dutton and Lambeck (2012) and Rovere et al. (2016) presented global data bases of LIg sea-level indicators, the majority of which were found to lie between 0 and 9 m above modern sea level. The observed variation in LIg SLIP elevation has widely been interpreted as suggesting that GMSL was not constant throughout the course of the LIg but was consistently well above present sea level (e.g., Rohling et al. 2019).
The suggestion of rapid sea-level fluctuations through the course of the LIg was disputed by Barlow et al. (2018) who suggested that the uncertainties associated with LIg sea-level indicators do not permit a robust sealevel curve reconstruction for this period and that there is, therefore, no robust evidence of sea-level variation through the LIg. However, this interpretation relies heavily on the resolution of the LIg observational data being sufficient to independently reflect the changes apparently evidenced by the oxygen isotope data. The rate at which geomorphological indicators of ice sheet geometry and RSL can form in response to rapidly evolving conditions directly limits the resolution that can be achieved. Where such indicators are able to form, their subsequent preservation through the last glacial cycle to the present further limits the spatial and temporal coverage of our observational data set. Another element of Barlow et al. 's argument is that glaciological modelling suggests that large-scale ice sheet accumulation is too slow to permit rapid fluctuations in GMSL. There is, however, mounting observational evidence that GMSL can rise and fall over relatively short periods (e.g., Yokoyama et al. 2018;Ishiwa et al. 2019;Fujita et al. 2019;Cutler et al. 2003;Yokoyama et al. 2001a, b;Deschamps et al. 2012).
Where LIg SLIPs fall below sea level, they become subject to a variety of erosional and sedimentation effects. Even where such features survive to the modern era SLIPs below modern sea level are difficult to locate and access. They are, additionally, often subject to geochemical conditions that make accurate dating impossible (e.g., Diagenesis, see Gallup et al. 1994;Stirling et al. 1995;Yokoyama and Esat 2015). All of which acts to inherently skew the observational record. This implicit selection bias should be carefully considered when interpreting the SLIP record from the LIg.
For sea-level observations from the LIg and earlier, the effects of sediment transport can accumulate over these longer time-scales to produce very significant changes in surface load (e.g., Rabineau et al. 2006). Substantial land motion may also occur over these time scales as a result of sediment compaction or the accumulated effect of tectonic processes (e.g., Pico 2020; Mitrovica et al. 2020). These processes are generally slow enough that on the time-scale of the most recent deglaciation, they may be neglected, but they may become significant on longer time scales and contaminate the sea-level record for the LIg.
The broad range of SLIP elevations and different inferred GMSL values reported by the studies listed above can potentially be reconciled by considering the effects of mantle dynamics on surface elevation (Fig. 10). Lithospheric elevation is strongly influenced by thermomechanical work done by the underlying mantle (e.g., Hager et al. 1985). Changes in the thermal conditions within the mantle and the motion of the continental plates relative to the mantle produce a change in the surface expression of these effects and a corresponding influence on long-term relative sea-level change. As crustal locations move toward regions with hot, upwelling underlying mantle they experience uplift and sea-level falls. This effect is slow enough to be negligible on shorter time-scales (< 0.05 mm/y), but can amount to 4-6 m on time scales of a 100,000 years with sea-level perturbations of different magnitudes and signs in different locations.
Theoretical and numerical models of mantle thermomechanical processes have been developed (e.g., Tan et al. 2006;Spasojevic et al. 2010), but these models do not produce consistent results (Guerri et al. 2016) and are still in the early stages of validation against observational evidence (Mitrovica et al. 2020;Rovere et al. 2015) which is sparse. Model uncertainties for these processes remain high and until model accuracy is significantly improved, the interpretation of LIg shoreline indicators from widely separated regions is problematic . At present, only RSL indicators in close proximity to one another (separated by length scales too small for  Fig. 10 Schematic of some of the processes associated with dynamic topography adapted from Austermann and Forte (2019). Since sea level is representative of geoid position relative to surface elevation, it is influenced by long-term mantle flow and relative motion between the mantle and the overlying lithosphere mantle dynamic effects to vary significantly) can be used to constrain the progress of sea level during the LIg with confidence. Widely geographically separated sea-level indicators from the LIg may be subject to differential mantle thermodynamic effects that perturb the SLIPs in opposite directions and are extremely difficult to observationally quantify.

Conclusions
Increased precision in sea-level observations, especially from the last glaciation to the present, has resulted from recent advancements in sampling and the analytical techniques applied to the result samples. Improved data have provided immediate constraints on the rate and magnitude of the interactions between climate-forcing and icesheet response. It has also permitted greater insight into paleo-ecological conditions (such as ocean salinity and temperature profiles, nutrient supply, turbidity etc.) during the geological past (e.g., Webster et al. 2018;Sanborn et al. 2020;Yagioka et al. 2019;Braga et al. 2019). The observational record continues to expand, and as more data become available, it seems to indicate that relative sea level is highly variable and very dynamic, at times contradicting the results of ice sheet modelling efforts (e.g., Yokoyama et al. 2018;Ishiwa et al. 2019). However, this variability suggests that more comprehensive observational coverage is necessary to fully understand sealevel change before and during the LGM at a resolution sufficient to more tightly constrain GMSL.
The effect of geological processes becomes increasingly significant on longer timescales and proximal to terrestrial ice sheets. The magnitude of the change in GMSL between the present and the LIg is comparable to the potential amplitude of mantle thermo-mechanical effects and this process must be considered when interpreting sea-level data sets from this period.
Continuing efforts to collect data from the near-, intermediate-and far fields of the major ice sheets will lead to the development of more detailed and accurate geophysical models. These improved models will permit us to better understand the rate and magnitude of ice sheet mass-loss resulting from climate forcing under a variety of basal conditions. This in turn will improve our understanding of the potential trajectory of future sea level and more effectively plan for the various contingencies that may emerge.
Models of GIA have improved significantly since their initial introduction and have permitted a more detailed understanding of sea-level dynamics. While the modelling community has begun to formalize the uncertainties attached to GIA modelling products, this process remains in its early stages and it remains frustratingly common for GIA models to be presented with no indication of how reliable their associated outputs are. This makes it somewhat difficult for researchers from outside the GIA community to determine which model might be best suited to their particular needs. It should be clearly recognized by the modelling community and their colleagues that GIA model results should be interpreted and applied with some caution, as there are fundamental uncertainties associated with GIA inversion techniques.
Using GIA-based inversion techniques to reconstruct ice sheet geometry requires a large amount of observational data. The data set employed must ideally be of high quality and have broad spatial and temporal coverage if separation of ice sheet and rheological model parameters is to be achieved. This requirement is not always met; sea-level indicators can only be generated at land-sea boundaries and require certain conditions to form. They may subsequently be eroded by glacial action or inundated as sea-level rises. Historically, land-based sea-level indicators have been easier to identify and study, whereas submerged sea-level indicators remain frustratingly sparse due to the practical difficulties involved in their study. This introduces an implicit bias into most sea-level data sets toward land-based studies that may be offset by a systematic program of ocean floor surveys and scientific ocean drilling. Communication between the observational data and modelling communities is necessary to forward these efforts (e.g., PALSEA of PAGES).