The Global Fingerprint of Modern Ice-Mass Loss on 3-D Crustal Motion

Crustal motion generated by rapid ice-mass loss from Earth's glaciers and ice sheets has previously been considered in Global Navigational Satellite System (GNSS) analyses and numerical models across regions of ice retreat. However, the fingerprint of ice-mass loss is not limited to glaciated but is characterized by a global pattern of 3-D crustal deformation. We compute “far-field” vertical and horizontal deformation rates that occurred in response to early 21st century mass flux from the Greenland and Antarctic Ice Sheets, global glaciers and ice caps, and associated ocean loading. We demonstrate that mass changes in the Greenland Ice Sheet and high latitude glacier systems each generated average crustal motion of 0.1–0.4 mm/yr across much of the Northern Hemisphere, with significant year-to-year variability

However, the impacts of ice-mass loss on the Earth system have a global reach. So-called "far-field" effects (which we loosely define as those that occur at a distance beyond 1000 km from the ice-mass loss) have long been recognized in GIA studies of sea-level change (e.g., Lambeck et al., 2014;Peltier & Andrews, 1976;Wu & Peltier, 1983) and crustal deformation (e.g., James, 1991;James & Lambert, 1993;Latychev et al., 2005;Mitrovica et al., 1994b). Moreover, in the context of modern global change, the rapid melting of ice sheets and glaciers drives significant, global-scale geographic variability in sea-level change (and its two bounding surfaces -the sea surface and solid earth surface) that have come to be known as sea-level fingerprints (e.g., Clark & Lingle, 1977;Mitrovica et al., 2001;Plag & Juettner, 2001).
The study of Tamisiea et al. (2003) was the first to predict the global scale fingerprints of modern ice-mass loss from the Greenland and Antarctica Ice Sheets on vertical land motion (henceforth VLM). More recently, Riva et al. (2017) compiled estimates of ice-mass changes in Greenland, Antarctica and mountain glaciers through the last century, and predicted VLM of several tenths of a millimeter per year across large continental areas. Patterns of VLM associated with modern ice-mass loss have since been observed in GNSS datasets, and have been estimated and corrected for in these data as well as in tide gauge measurements (Schumacher et al., 2018;Frederikse et al., 2019). Both predictions and observations of VLM from modern ice-mass loss indicate large temporal variability, over time periods of several years, as well as significant spatial variability (Frederikse et al., 2019;Riva et al., 2017). Finally, Ludwigsen et al. (2020) demonstrated that the VLM signal of Arctic ice-mass loss reached amplitudes of 0.5-1.0 mm/yr over the 2003-2015 time period across high latitude polar regions.
Despite the recent emphasis placed on accurately modeling and analyzing VLM driven by modern icemass loss, no previous studies have estimated the accompanying horizontal crustal motions. This omission is surprising given that published work on GIA-induced crustal deformation (e.g., James, 1991;James & Lambert, 1993;Latychev et al., 2005;Mitrovica et al., 1994b) indicate that horizontal motions can exceed their vertical counterparts in far-field regions, and, additionally, uncertainties in GNSS estimates of the former are typically a factor of 2-3 smaller than the latter. Generating models that match observations of horizontal motion due to ongoing GIA has been an outstanding problem in several regions. For example, while calculations in Scandinavia based on 1-D viscoelastic Earth models have accurately reconciled 3-D crustal deformation rates estimated from the BIFROST GNSS network (e.g., Johansson et al., 2002;Milne et al., 2001;Milne et al., 2004), analogous GIA model-data comparisons for North American sites away from tectonic influence have yielded good fits to vertical rates but large misfits in the horizontal component of deformation (Argus et al., 1999;Sella et al., 2007). There is evidence that these misfits may arise from limitations in simplified Earth models used in the GIA simulations (Latychev et al., 2005;Peltier & Drummond, 2008), but a full understanding of the cause of these residuals has yet to emerge.
Herein, we present global predictions of 3-D crustal motions in response to mass changes in polar ice sheets and glaciers over the early 21st century, though we focus our discussion on the far-field signal. Our calculations include the associated sea-level changes in specifying the full surface mass load.

Predictions of 3-D Crustal Motion
Our predictions of three-dimensional deformation arising from surface ice-mass and ocean loading adopt a spherically symmetric (i.e., depth varying) Earth model. The calculations are based on elastic Love number theory (Farrell, 1972;Mitrovica et al., 1994a). We adopt the elastic and density structure of the Earth given by the seismic model PREM (Dziewonski & Anderson, 1981) and incorporate rotational effects using now-standard approaches (Milne & Mitrovica, 1998;Mitrovica et al., 2005). Gravitationally self-consistent sea-level changes are included using the pseudo-spectral algorithm described by Kendall et al. (2005) with a truncation at spherical harmonic degree and order 512. Crustal motions are referenced relative to an origin located at the center of mass of the deformed planet (Mitrovica et al., 1994a). To gain insight into the physics of far-field effects, we individually consider horizontal and vertical surface motions generated by recent mass loss from three of the main contributors to modern-day sea-level rise: GrIS, AIS and a global database of mountain glaciers and ice caps. Total crustal deformation rates can be computed by simply summing these contributions, and are discussed at the end of Section 2. Figure 1a shows the predicted average annual surface motion in response to recent GrIS mass loss, computed using a published estimate of average mass balance from satellite altimetry datasets. Specifically, we adopt mass changes in GrIS (2003GrIS ( -2018 from an analysis of NASA's ICESat and ICESat-2 data given in Smith et al. (2020). This allows prediction of the average crustal motions due to ice and ocean mass changes over the 16 years time period. As explored in detail in recent literature, within south and west Greenland and bordering the ocean, average crustal uplift rates on the order of several mm/yr are predicted (e.g., Simpson et al., 2011;Spada et al., 2012). In contrast, in northeast Greenland, an area of ice accumulation, sites are predicted to subside at rates up to 1 mm/yr. At slightly further distance, for example, across Iceland, Baffin Island and other sites in Arctic Canada, the average predicted uplift rates lie in the range of 0.4-1.0 mm/yr. There is a further decrease in magnitude, to 0.0-0.4 mm/yr, across the rest of Canada and 0.0-0.2 mm/yr in Fennoscandia and northern Europe. Our global predictions are characterized by minor subsidence (0.1 mm/yr) across the equator and uplift at high southern latitudes.

Greenland
Significant tangential motions are predicted across North America, the North Atlantic Ocean and northern Eurasia, with large north-south gradients in magnitude. Tangential motions largely follow a northward trend in the Atlantic Ocean, turning to the east and west above 40°N. Moving from northeast Canada toward the western US, the magnitude of the average horizontal rates decreases gradually from 0.3 mm/yr to less than 0.05 mm/yr. Similarly, average horizontal rates decrease from maximum values near 0.2 mm/yr in Norway to less than 0.05 mm/yr in southern Europe.
The panels of Figure 1b focus on the variability of these trends through time over North America and Europe. These predictions are calculated using maps of Greenland ice-mass change derived from GRACE gravity field coefficients (CSR, Release 6.0, University Of Texas Center For Space Research (UTCSR), 2018). Following the methodology of Harig and Simons (Bevis et al., 2019;Harig & Simons, 2012), we transform the GRACE geopotential coefficients into surface mass density, subtracting a modeled signal associated with GIA, and project the resulting fields onto a Slepian basis to reduce the signal-to-noise ratio. This allows for the calculation of gridded yearly ice histories for the period 2003-2015.  (2006) ice-mass loss. In years of high mass loss, VLM rates in northern (mainland) Canada and Scandinavia reach peaks of 0.9 and 0.7 mm/yr respectively, decreasing rapidly to the south. These VLM rates are more than double those predicted in years of low ice-mass loss from Greenland. Tangential motion also varies dramatically, reaching peak rates of 0.45 mm/yr in northernmost regions and magnitudes of 0.2-0.3 mm/yr across large sections of North America and Europe in 2012, as compared to rates of 0.1 mm/yr across both continents in years of low mass loss. We note that for both high and low mass loss years, predicted horizontal motions across most of Canada and parts of Europe exceed the predicted vertical motions that have been a focus of recent geophysical literature (Frederikse et al., 2019;Ludwigsen et al., 2020;Riva et al., 2017).

Glaciers and Ice Caps
Figures 2 and 3 explore crustal motion generated by global glacier and ice cap ice-mass loss from 2003-2013. The predictions are calculated using a database of yearly mass balance for 17 glacier systems generated from observation-validated climate modeling (Marzeion et al., 2012). We pair the mass balances with static glacier outlines from the Randolph Glacier Inventory (RGI Consortium and others, 2017) to establish a complete history of glacier (and ice cap) evolution, and solve for sea-level variations to specify the full surface Yearly spatial ice-mass loss is inferred from an analysis of GRACE gravity data (see text). Vertical motion is given by the red-blue color scale and tangential motions are represented as arrow vectors. Tangential motions in the vicinity of Greenland are excluded in order to highlight the far-field response. mass load. The global crustal motion field predicted from this exercise is largely dominated by ice-mass changes in the Arctic and northern North America.
The pattern of deformation directly surrounding the Arctic glaciers, and their effect on northern Europe, Asia and North America is shown in Figure 2. Outside of regions in close proximity to ice-mass changes, peak crustal uplift rates reach 0.1-0.3 mm/yr across northern North America and Eurasia, depending on the yearly average glacial mass balance of Arctic glaciers. Average tangential motions in these regions are largely directed northward toward the Arctic glaciers and have magnitudes of 0.1 mm/yr, although there is some influence of ice-mass loss from Alaskan glaciers (Figure 3). Predictions of horizontal motion show significant year-to-year variations in magnitude and direction. In years of maximum ice-mass loss, tangential motions across North America can reach 0.15 mm/yr (excluding the region directly surrounding the Alaskan glaciers). All predictions show large north-south gradients in tangential motion across the North American and Eurasian continents.  localized, outward-directed spoke-like pattern typical of post-glacial uplift, whereas in both Alaska and Himalayas, the pattern of horizontal motion reflects the more complex ice-mass loss geometries in these regions. Figure 3 demonstrates that outside of the localized regions of mass loss, horizontal motion is largely dominated by northward-directed motion towards the Arctic glacier systems.
COULSON ET AL.

Antarctica
Finally, we predict crustal deformation rates generated by a mass flux history across the Antarctic Ice Sheet constructed by Martín-Español et al. (2016). These are shown in Figure 4. The Martín-Español et al. (2016) ice history was derived by combining observations of satellite gravimetry and altimetry, as well as surveying using GNSS, over the period 2003-2013. We find that these ice-mass changes led to far-field crustal motions that were smaller than those driven by GrIS or high latitude glacier systems. While relatively significant mass loss from ice sheets and glaciers occurred across West Antarctica and the Antarctic Peninsula, these losses were partially compensated by mass gain in East Antarctica between 2003(Martín-Español et al., 2016. The variation in geometry and net magnitude of mass change across the AIS was also highly variable through time; for example, while some areas of both the East and West Antarctic ice sheets accumulated mass during years of low net ice-mass loss (e.g., 2005), they experienced significant mass loss during later years (e.g., 2010). The net effect is that the magnitude of predicted crustal motions averaged across the 2003-2013 time period is muted. The largest far-field crustal motions generated by AIS mass changes occur at sites across the southern Pacific Ocean (Figures 4 and 5).

Summary of Total 3-D Crustal Motion
Sections 2.1-2.3 and the associated figures provide the global 3-D deformation field induced by mass loss from GrIS, AIS and a global database of glaciers and ice caps. Figure 5 shows each of these individual contributions and their total from 2003-2013, for a selection of five sites. The temporal trends in both predicted vertical and tangential motions exhibit considerable year-to-year variability that reflects the mass balance inferred for the various regions of ice. For example, during years of high ice-mass loss from GrIS, predicted horizontal motions at sites across the northern hemisphere are up to a factor of 3 larger than those during low mass loss years.
As demonstrated in Figures 1-4, both horizontal and vertical deformation rates are highest at sites close to the area of mass loss. Deformation rates at Norilsk in Figure 5 illustrate this effect; these rates show a large contribution from glacial mass loss since the site is proximal to the Russian Arctic Island and Svalbard glacier systems, which are experiencing rapid retreat. Crustal rates from 2003-2013 at London, Boston and Norilsk are heavily impacted by the signal from mass loss from GrIS. While GrIS-induced deformation rates decrease towards the equator, these three sites demonstrate the non-negligible sensitivity to GrIS ice-mass loss across the Northern Hemisphere at relatively large distance from the ice sheet. As shown in Figure 4, ice-mass loss from AIS has a relatively muted deformation signal, even at proximal sites (e.g., Falklands).
COULSON ET AL.

Discussion and Conclusions
Crustal motions are a key observable in a multitude of geodynamic applications, including estimating icemass changes, constraining Earth models, and isolating GIA signals in the crustal deformation field, as well as in a range of tectonic studies. The predictions presented herein illustrate the range of signals induced by early 21st century ice-mass loss that have previously been overlooked in GNSS and tide gauge studies and therefore have the potential to improve a range of geophysical analyses.
Several studies have investigated the signature of recent ice-mass loss on far-field vertical land motion (e.g., Ludwigsen et al., 2020;Riva et al., 2017). Despite some differences in methodology, our predictions show good agreement with the spatial variation and magnitude of vertical deformation rates presented by Riva et al. (2017). Moreover, our predictions support the earlier study's conclusion that it is important to account for the VLM signal in GNSS data at sites well away from present-day ice cover, particularly in the Northern Hemisphere where motions generated by mass loss from GrIS and Arctic glaciers are large. In the present study, we have extended this earlier work to show the first global predictions of the full 3-D crustal motion field induced by recent ice-mass changes in the polar ice sheets and mountain glaciers and ice caps. Using reconstructions of Greenland ice history over the period 2003-2018, we predict average vertical motions of up to 0.4 mm/yr across the Northern Hemisphere. Tangential rates of 0.05-0.3 mm/yr are predicted across most of Canada and the US, and rates in the range 0.05-0.2 mm/yr are predicted in Fennoscandia and Europe. The ice-mass loss from Arctic glaciers also produces pervasive horizontal motions with magnitude up to 0.15 mm/yr (average 2003-2013) across the high latitudes. Both of these sites of ice-mass change drive large gradients in horizontal and vertical motion across the northern hemisphere, with significant year-toyear variability in their magnitude and direction. These far-field horizontal rates have not been accounted for in any study of GNSS data, despite the higher accuracy of horizontal versus vertical rate determinations.
Crustal motion induced by recent ice-mass loss will have an impact on the terrestrial reference frames derived from global geodetic measurements. A model for these signals is generally not included in the determination of the International Terrestrial Reference Frame (ITRF) (e.g., Altamimi et al., 2011Altamimi et al., , 2016. Future analyses of GNSS networks must therefore account for reference-frame errors by including adjustments to the a priori reference frame (e.g., Hill et al., 2010;Kierulf et al., 2021).
The magnitude of both horizontal and vertical rates driven by modern ice-mass loss suggests that a reappraisal of some previous results may be warranted. As an example, analyses of GNSS data from sites across the BIFROST network have inferred viscoelastic structure beneath Fennoscandia from the estimated 3-D COULSON ET AL. 10.1029/2021GL095477 8 of 11 crustal motion rates (e.g., Kierulf et al., 2014;Milne et al., 2004) and have used the vertical component of these rates to correct local tide gauge (sea-level) records for the contaminating signal of GIA (e.g., Davis et al., 1999;Johansson et al., 2002). The combined fingerprint of elastic crustal deformation due to recent changes in the mass of the GrIS, AIS and global glaciers and ice caps (shown in Figure S1) shows spatially systematic variations in average deformation rates of 0.3 mm/yr for the vertical component and 0.1-0.2 mm/yr for the horizontal component, across Fennoscandia. A recent analysis of BIFROST GNSS observations by Kierulf et al. (2021) yields uncertainties for long-term rates at individual stations of comparable or smaller magnitude than these predicted variations (typically  0.1 mm/yr for the horizontal component and 0.1-0.5 mm/yr for the vertical component). Given the dense distribution of sites, accuracy of the GNSS measurements, and unique spatial and temporal variations in deformation associated with the present-day mass changes, inclusion of both components of this deformation in a network analysis therefore has the potential to significantly alter results of geodetic analyses. This would impact, for example, any inferences of viscoelastic Earth-model parameters. The unique time history of modern ice-mass flux relative to other contributing signals will aid in the detection of these processes in geodetic data.

Acknowledgments
Sophie Coulson and Jerry X. Mitrovica acknowledge funding from the Star-Friedman Challenge and Lemann Brazil Research Fund, through Harvard University. Sophie Coulson was also supported by Frank Knox Memorial Fellowships. Mark J. Hoggard acknowledges support from the Australian government's Exploring for the Future program. We thank Chris Rollins and one anonymous reviewer for constructive reviews that improved the final manuscript.