GNSS Observations of GIA‐Induced Crustal Deformation in Lützow‐Holm Bay, East Antarctica

The Japanese Antarctic Research Expeditions have conducted long‐term Global Navigation Satellite System (GNSS) observations at a series of bedrock outcrops along the coast of Lützow‐Holm Bay, East Antarctica, to elucidate the solid Earth response to ice mass changes. These GNSS observations capture both the viscoelastic response of the solid Earth to Antarctic Ice Sheet changes since the Last Glacial Maximum, termed glacial isostatic adjustment (GIA), and elastic deformation associated with recent surface mass changes due to snowfall variations. Here, we extract the GIA signals from the Lützow‐Holm Bay GNSS data by applying an elastic deformation correction based on the mass fluctuations observed in the Gravity Recovery and Climate Experiment and the satellite altimetry. Our results indicate that the region is currently experiencing approximately 1 mm/yr of subsidence, which is large enough to influence GIA model estimations. The resultant GIA‐induced deformation response in Lützow‐Holm Bay is 1–6 mm/yr of uplift.

Antarctic GIA is poorly constrained due to the lack of in situ geological and/or geodetic observations. This means that the key parameters in modeling the GIA response, such as lithospheric thickness and mantle viscosity, remain poorly constrained. Numerous ice-history and mantle-viscosity models have been proposed to predict Antarctic GIA, even though they are largely dependent on these poorly constrained model parameters (e.g., Whitehouse et al., 2019). Therefore, it is necessary to evaluate which GIA model best represents the true crustal deformation to accurately capture and understand current AIS mass changes.
Geodetic observations, such as Global Navigation and Satellite System (GNSS) and gravity measurements, measure present-day deformation and provide constraints that should be satisfied by the GIA model. Many studies have employed GIA-related geodetic observations to compare and evaluate GIA model predictions in Antarctica (e.g., Thomas et al., 2011), Greenland (e.g., Khan et al., 2016), Scandinavia (e.g., Kierulf et al., 2014), and North America (e.g., Sella et al., 2007) all of which have undergone massive ice-sheet-scale retreat/collapse in the past.
Large GNSS networks have been installed in West Antarctica, where the ice sheet is rapidly melting, whereas there are few GNSS sites across East Antarctica. One of the reasons for the poor coverage in East Antarctica is that most Antarctic research stations are concentrated in West Antarctica, especially the Antarctic Peninsula. However, the GIA-induced crustal deformation response in East Antarctica is an essential research target due to the much greater mass of the East Antarctic Ice Sheet compared with that of the West Antarctic Ice Sheet.
Lützow-Holm Bay, which extends ∼150 km from north to south, is located between 37°E-40°E and 69°S-70°S. The Japanese Antarctic Research Expedition (JARE) has been conducting geodetic observations around Lützow-Holm Bay (Shibuya et al., 2003) since 1993. Ohzono et al. (2006) reported on GNSS observations from a series of outcrop sites in the region until 2004, and Kazama et al. (2013) reported on absolute gravimetry at Syowa Station and Langhovde. Furthermore, geomorphological studies have revealed the detailed melting history of the ice sheet at some regional outcrop sites since the LGM (Kawamata et al., 2020).
Here we interpret the GIA signal in Lützow-Holm Bay, East Antarctica, using GNSS and GRACE observations. We estimate the vertical deformation rate induced by recent and past ice-sheet mass changes in this area from the GNSS position time series. We then calculate the elastic deformation caused by the recent GRACE-derived ice-sheet mass changes. The GRACE observations indicate a surface mass gain in Dronning Maud Land, East Antarctica, which includes the study area; this gain appears to be attributed to increased snowfall. The crustal deformation signal derived from the GNSS observations includes the elastic deformation response due to this current redistribution of mass. The results demonstrate that it is necessary to apply an elastic deformation correction to the GNSS solutions to accurately estimate the GIA-induced deformation response.
Finally, we compare the current deformation rate, which is derived from a GNSS analysis with an elastic correction, with numerical predictions from four GIA forward models the ANU (Lambeck et al., 2014), ICE-5G (Peltier, 2004), ICE-6G (Peltier et al., 2015), and IJ05 (Ivins et al., 2013 models) to reveal the GIA-induced deformation in Lützow-Holm Bay.
We analyzed the GNSS data obtained from four outcrop GNSS sites (Langhovde, Skarvsnes, Padda, and Rundvågshetta) along the coast of Lützow-Holm Bay and International GNSS Service (IGS) SYOG, the continuous GNSS station at Syowa Station.
SYOG was established in 1995 and has since been maintained by the Geospatial Information Authority of Japan; a continuous time series with no long interruptions has been acquired to date. The antenna and reciever replacement of SYOG was conducted for multi-GNSS observation at the end of 2013. Conversely, the outcrop GNSS sites are discontinuous campaign GNSS sites, which were mainly installed by the JARE geoscience group. JARE has been conducting GNSS campaigns across the outcrop sites since 1998. Ohzono et al. (2006) first reported on the geodetic results from these sites using the 1998-2003 campaign GNSS observations. Campaign GNSS observations have been acquired at the outcrop sites for over two decades. However, campaign observations generally have the following problems during the geodetic analysis and interpretation: (a) It is necessary to obtain precise antenna height information at each GNSS campaign site, as inaccurate antenna heights will increase the vertical uncertainties. (b) It is difficult to distinguish annual and/ or seasonal variations from the long-term crustal deformation trend when the campaign observations are only acquired once or twice a year. An automated observation system has been gradually installed at these outcrop sites to provide continuous GNSS measurements, thereby resolving the abovementioned problems. The four outcrop sites at Langhovde, Skarvsnes, Padda, and Rundvågshetta are sites where geodetic data have been successfully acquired after automation. Each automated outcrop site is equipped with an on-off timer that is set to power the GNSS receiver (GEM-1), which is connected to a JAVA GrAnt-G3T antenna, for 24 h each week. A lead-acid battery is also installed at each site to allow the acquisition of continuous observations, even during the polar night, with a solar power generation system meeting the power needs during the summer season. The maintenance of this system and data collection are conducted once or twice a year, mainly during the JARE summer observation activities.
The position of each GNSS site with respect to IGS14 was estimated using the "c5++" multi-technique space geodetic analysis software (Otsubo et al., 2019), which employs the precise point positioning (PPP) static method, and using the antenna model IGS14.atx. We used the ionosphere-free combinations of GPS L1/L2 measurements and did not apply the higher-order ionosphere correction. The station position was estimated every 24 h, and the receiver clock offset was estimated every 30 s. The IGS Final products were used for the GNSS satellite orbit and clock information. The analysis was conducted after removing the solid Earth tides of the major tidal constituent (from semi-diurnal to 18.6 years) and oceanic tidal loading using the EOT11a model (Savcenko et al., 2012). The non-tidal ocean and atmospheric loading deformations were HATTORI ET AL. Note. The GNSS results show the linear trend and estimation errors in each direction, which were estimated via the least squares method based on the PPP solutions at each GNSS site. The "Corrected Δu" column provides the vertical component of the GNSS result after applying the elastic deformation correction. ANU and ICE-5G indicate the GIA model prediction ranges for the ANU (Lambeck et al., 2014), ICE-5G (Peltier, 2004), and ICE-6G (Peltier et al., 2015) ice models with different Earth models, respectively. Abbreviations: GIA, glacial isostatic adjustment; GNSS, Global Navigation Satellite System; PPP, precise point positioning.

Table 1
Summary of the GNSS Results corrected using the International Mass Loading Service data set (Petrov, 2015). We used the Global Mapping Function (Boehm et al., 2006) for the troposhperic delay correction. Figure 1b shows the velocity at each site estimated by fitting its GNSS position time series. The velocity of SYOG was estimated using Hector software (Bos et al., 2013) with flicker noise and white noise model. The time series of SYOG shown in Figure 1b was corrected for the offset due to the equipment update in 2013. The velocities of other outcrop sites were estimated using a linear regression model, with the annual and semi-annual components after removing outliers that were more than 10 cm away from the median. where t i is the observation epoch, x(t i ) is the position at t i , x 0 is the intercept, v is the estimated velocity, doy is the day of year, a, b, c, and d is amplitudes of each periodic terms.
The SYOG velocity was estimated based on the daily solutions. The velocity at each outcrop site was estimated using the weekly solutions and some campaign observations before automation. Table 1 summarizes the measured velocities and associated uncertainties at each site.
Although we estimated the velocities for the east-west (E-W) ( Figure S1), north-south (N-S) ( Figure S2), and vertical (U-D) components during our GNSS analysis, we focused only on the vertical velocity results in our interpretation, and did not discuss the horizontal components. The reason for this omission of the horizontal velocity results is twofold: (a) it is difficult to separate the variations in the horizontal signals due to plate rotation and GIA from the GNSS survey results acquired in this narrow region , particularly when the GIA-induced horizontal deformation in this region is expected to be small (Turner et al., 2020); and (b) the GIA component is dominant in the vertical direction. It should be noted that the vertical velocities were estimated at each site, even though the velocity magnitude varied among the sites. Furthermore, the estimated vertical velocities increased with increasing latitude and proximity to the AIS at four of the sites (SYOG, Skarvsnes, Padda, and Rundvågshetta). The estimated deformation (both vertical and horizontal) at Langhovde is different from that at the other outcrop sites. It is therefore possible that accurate crustal deformation could not be extracted for the Langhovde GNSS observations due to fitting errors, such as the entire data set being overinfluenced by the early campaign observations.

Elastic Deformation
We calculated the elastic deformation due to the current AIS mass variation at each GNSS site. It is necessary to apply the appropriate elastic deformation correction to accurately extract the GIA-induced crustal deformation because the vertical displacement velocity obtained via the GNSS analysis is the sum of the GIA-induced viscoelastic component and the elastic component, which reflects recent changes in surface mass.
GRACE has recently observed a widespread gain in surface mass over the Drowning Maud Land area, including Lützow-Holm Bay (Velicogna et al., 2020). This mass gain indicates a recent increase in regional snowfall.
We used the GRACE mascon solution (JPL GRACE Mascon Release 06: Watkins et al., 2015;Wiese et al., 2018), which has been corrected with a GIA effect based on the ICE6G-D model (Peltier et al., 2018), and the satellite altimetry data (Schröder et al., 2019) as the recent mass load change in the elastic deformation calculation. We calculated the mass change from the surface elevation change observed satellite altimetry with firn density model (Ligtenberg et al., 2011) and the density of ice (917 kg/m 3 ). Since the firn density corresponds to the lower limit and the ice density corresponds to the upper limit, the actual mass change is to be between these two assumptions. A mesh was created to divide the Antarctic into uniform 2.8 × 2.8 km grid cells, and a disk load that corresponded to the surface mass change at each grid cell was set on each mesh. The elastic deformation due to each disk load was calculated in the reference frame of the center of mass of the solid earth (CE) via Farrell's Green's function method (Farrell, 1972), with the load Love number calculated from the Preliminary Reference Earth Model (Dziewonski & Anderson, 1981) using the density and elastic constants of the Earth's interior (Martens et al., 2019). We computed the Green's function using load Love numbers up to the 10,000th order. The GRACE mascon solution shows that the region around Lützow-Holm Bay experienced a constant surface mass gain of ∼10 cm/yr (water equivalent) from 2002 to 2017. The calculated deformation indicates that this region is experiencing ∼1 mm/yr of subsidence due to elastic deformation. The differences in the total deformation among the GNSS sites approach 5 mm over the 15-years analysis period across the Lützow-Holm Bay region.
As with GRACE, the satellite altimetry data shows a continuous surface elevation increasing. Although the calculated elastic deformation varies greatly depending on the assumed density, the similar time series HATTORI ET AL.
10.1029/2021GL093479 6 of 10  (Wiese et al., 2018) was used as the surface mass variation in both (a and b). The native resolution is 3° × 3°, and the data product is interpolated to a 0.5° × 0.5° grid.
(c) Spatial distribution of surface elevation change in Lützow-Holm Bay estimated from satellite altimetry (Schröder et al., 2019). (d) Time series of surface elevation change (SEC) in the grid surrounded by the dotted black line in (a/c) shown by 2d-histogram (bin-width: 2 month × 5 cm) with (upper) and the corresponding elastic deformation at SYOG (the northernmost aite) and Rundvågshetta (the southernmost site) under two assumptions of snow density: the firn density model (Ligtenberg et al., 2011) and the density of ice 917 kg/m 3 (lower).
were shown. The elastic deformation rate with the firn density model is about 0.5 mm/yr, which is about half of the GRACE-derived elastic deformation, and using ice density assumption shows the similar rate as GRACE.
The GNSS-derived vertical displacement velocity after applying the elastic deformation correction is listed for each GNSS site in Table 1. Here we consider these corrected values to be the GIA-induced vertical velocity due to past changes in ice-sheet mass. The results highlight the GIA-induced uplift trends across this region.

Discussion
The elastic deformation due to recent ice-sheet mass changes shows that the Lützow-Holm Bay region is currently subsiding at ∼0.5-1 mm/yr. We applied this elastic correction to the GNSS-derived vertical deformation rate, and obtained >1 mm/yr of uplift as the estimated GIA-induced signal at each GNSS site. Our results reveal that elastic deformation has a significant effect on the GIA signal, because the calculated elastic deformation and GNSS-derived vertical deformation possess similar amplitudes. Therefore, it is important to apply the elastic deformation correction to ensure the accurate evaluation of the GIA signal.
Our GNSS results after applying the elastic deformation correction indicate that the Lützow-Holm Bay region is experiencing 1-6 mm/yr of uplift (Table 1 and Figure 3). The estimated uplift rates show an increase in the vertical displacement velocity to the south, except for Langhovde. It is assumed that this spatial pattern reflects the past AIS configuration.
HATTORI ET AL.
10.1029/2021GL093479 7 of 10 The dark-green dots indicate the GNSS estimations with the elastic deformation correction due to the recent GRACEderived surface mass variations (Section 3), which are assumed to be the GIA signals. The blue, orange, red, and lightblue lines show the GIA predictions based on the ANU (Lambeck et al., 2014), ICE-5G (Peltier, 2004), ICE-6G (Peltier et al., 2015), and IJ05 (Ivins et al., 2013) which is combined with the Northern Hemisphere part of ICE-6G, with different Earth models, respectively. Estimated and modeled vertical deformation velocities at each GNSS site (upward is positive).
The elastic deformation due to recent variations in surface mass, as estimated from GRACE, has a coarse resolution because the spatial resolution of the GRACE mascon solution is only about 300 km (Wiese et al., 2018). Consequently, it does not accurately reflect the difference in the local mass change distribution (<300 km) near each GNSS site. This coarse resolution makes it difficult to accurately interpret the result from each GNSS site. On the other hand, although the elastic deformation calculated from the altimeter has a high spatial resolution, the deformation rate depends on the assumed snow density when converting the height change into mass change. However, it does not have a significant impact on discussions concerning the entire Lützow-Holm Bay region, which shows an average subsidence rate of ∼1 mm/yr.
The Langhovde GNSS uplift rate is much higher than other sites. Although the cause of this high rate is not well understood, one plausible reason is that the observation period and data coverage are not long/ large enough to detect the long-term trend. The same reason will explain why the horizontal component at Langhovde differs from those at the other sites. The SYOG result is the most reliable because it is obtained from the longest (and most continuous) observation duration among the GNSS sites. It is therefore necessary to acquire longer GNSS time series at the outcrop sites to more accurately understand local crustal deformation.
Finally, we compared the vertical deformation rate obtained in this study with the GIA predictions based on the ANU (Lambeck et al., 2014), ICE-5G (Peltier, 2004), ICE-6G (Peltier et al., 2015), and IJ05 (Ivins et al., 2013) ice models, which are well-established GIA models. The predicted current vertical displacement at each GNSS site was estimated using a typical three-layer viscosity model that was characterized by the elastic thickness, and upper (above 670 km depth) and lower-mantle viscosities. We used an elastic thickness of 50-100 km, and upper-and lower-mantle viscosities of 2-10 × 10 20 Pa s and 10 21 -10 23 Pa s, respectively. At secular timescales, it is necessary to correct for shifts between the reference coordinate systems (CE-CM correction) when comparing GNSS-derived deformation (CM) with GIA model predictions (CE) (Riddell et al., 2020). We applied CE-CM corrections of −0.34 mm/yr (Sun & Riva, 2020) to GNSS results and compare them on CE. Figure 3 compares the GNSS vertical displacement velocities obtained in this study with the two GIA model predictions. The GIA model predicted displacement velocities that were approximately in the range of −1 to +3 mm/yr for all of the sites. Large discrepancies are observed between the GNSS-derived GIA signals and the predictions of ANU and ICE5G models. ICE6G and IJ05 model predicted the similar deformation rate to the GNSS results. Although we observed an increase in the displacement velocity with increasing latitude in the GNSS results, this spatial characteristic was not clearly observed in the GIA model predictions.
Here, we considered the GNSS-derived vertical velocities to represent GIA-induced deformation after applying elastic deformation corrections. However, it is noted that the GIA signal extracted in this study does not take into account other deformation factors, such as post-seismic deformation. Furthermore, our results would contain elastic deformation due to very local load variations near each site and the GNSS analytical errors caused by our limited observations.
Lützow-Holm Bay is considered an aseismic region, with no reported earthquakes that could have potentially caused large post-seismic deformation since the 1980s (Kanao, 2014). King and Santamaría-Gómez (2016) reported that the post-seismic deformation caused by a large earthquake (M W 8.1) near Balleny Island on March 25, 1998, is ongoing across a large part of Antarctica. However, the post-seismic deformation caused by this earthquake is relatively small in Lützow-Holm Bay, due to its distance from the epicenter, and may act as an offset in the vertical components throughout Lützow-Holm Bay. Therefore, care must be taken when discussing the absolute value of the vertical deformation, whereas the spatial distribution of the vertical velocity in Lützow-Holm Bay is not significantly affected, thereby allowing us to discuss past AIS changes and the Earth's internal structure based on these spatial differences across the region.
The abovementioned problems highlight the need to acquire additional geodetic observations at more sites to obtain accurate GIA signals from in situ observations. It is imperative that elastic deformation corrections are considered in greater detail to make an accurate comparison of the GNSS results and GIA model predictions. It is also worth noting that the GIA predictions for the ice-sheet margins are strongly dependent on ice-sheet models. Therefore, additional GIA model predictions that employ other ice models (e.g., W12, Whitehouse et al., 2012) are required to further investigate the observed discrepancies in the GNSS-derived and modeled GIA signals.

Conclusions
We analyzed the GNSS observations from SYOG at Syowa Station and four outcrop GNSS sites (Langhovde, Skarvsnes, Padda, and Rundvågshetta) in the coastal area of Lützow-Holm Bay via the PPP static method to investigate the response of the region to post-LGM retreat and recent mass changes. We applied an elastic deformation correction due to recent Antarctic surface mass fluctuations to the GNSS vertical displacements to obtain the GIA-induced uplift velocities for the region.
Our results show that the elastic deformation response in the Lützow-Holm Bay region is about 1 mm/yr of subsidence due to the recent increase in snow cover; the GIA-induced deformation response is 1-6 mm/yr of uplift after applying an elastic deformation correction to the GNSS observations.
Although it is desirable to apply more detailed corrections, such as post-seismic deformation, to discuss the absolute values of the vertical deformation rates, the ICE-6G model yielded the most consistent results among the three GIA models simulated in this study. However, none of the models clearly showed the N-S spatial characteristics that were detected by the GNSS results. This discrepancy between the model and observational results highlights the need to discuss such spatial features using higher-resolution models of past AIS mass placement and the Earth's internal structure.

Data Availability Statement
The SYOG rinex data were obtained through the online archives for the Crustal Dynamics Data Information System (CDDIS), NASA Goddard Space Flight Center (https://cddis.nasa.gov/archive/gnss/data/ daily/). The other outcrop sites rinex files used in this study are available from the NIPR site (http://geo. nipr.ac.jp/?page_id=1035). By clicking on site on the map, a link of rinex files obtained at the outcrop site will be shown.