Suvremena Tehnologija Kao Podrska Edukaciji U Turizmu

저작자표시-비영리-변경금지 2.0 대한민국 이용자는 아래의 조건을 따르는 경우에 한하여 자유롭게 l 이 저작물을 복제, 배포, 전송, 전시, 공연 및 방송할 수 있습니다. 다음과 같은 조건을 따라야 합니다: 저작자표시. 귀하는 원저작자를 표시하여야 합니다. 비영리. 귀하는 이 저작물을 영리 목적으로 이용할 수 없습니다. 변경금지. 귀하는 이 저작물을 개작, 변형 또는 가공할 수 없습니다. l l 귀하는, 이 저작물의 재이용이나 배포의 경우, 이 저작물에 적용된 이용허락조건 을 명확하게 나타내어야 합니다. 저작권자로부터 별도의 허가를 받으면 이러한 조건들은 적용되지 않습니다. 저작권법에 따른 이용자의 권리는 위의 내용에 의하여 영향을 받지 않습니다. 이것은 이용허락규약(Legal Code)을 이해하기 쉽게 요약한 것입니다. Disclaimer Ph.D Dissertation Research of Antarctic grounding line migration in high resolution using hydrostatic ice thickness derived by radar satellite data 레이더 위성 자료와 유체정역학적 얼음 두께를 이용한 남극 고해상도 지반접지선 변동 연구 August 2018 Graduate School of Earth and Environmental Sciences Seoul National University Seung Hee Kim Abstract Grounding line provides crucial information about stability of ice shelves and is also used as a boundary condition for ice sheet modelling. Monitoring grounding line migration is one of the primary examinations for determining vulnerability of ice shelves. In recent years, grounding line estimation and monitoring has been carried out using a number of techniques utilizing satellite systems and has confirmed grounding line retreats in most glaciers in West Antarctica and some areas in East Antarctica. The present study too was undertaken to monitor grounding line migration of Antarctic ice shelves. The most important aspect, however, was that monitoring must be done in high temporal and spatial resolution, and it thus required a simple and repeatable monitoring technique. Hence, it was decided that the most appropriate method was to exploit surface topography of ice because the basal variations are transmitted to the surface. The topography information can be obtained in high resolution with short time intervals due to recent development of SAR satellite systems. Finally, this study has shown that the hydrostatic ice thickness derived from surface topography can be effectively utilized to estimate grounding line position where the bottom of ice starts to diverge from the bathymetry. It required independent datasets of surface elevation of ice and sea level to calculate hydrostatic ice thickness and the bathymetry to determine where ice starts to float. Surface topography of Antarctica was successfully generated in high spatial and vertical resolution by means of TanDEM-X and TerraSAR-X i bi-static interferometry. Due to the nature of continuous variation of surface conditions, securing reliable ground reference elevation had been a major problem in adopting interferometry to generate surface topography map over snow and ice covered surface. In order to solve this problem, near-coincidental CryoSat-2 SARIn radar altimetry data were used as ground control points after careful selection using normalized cross-correlation and standard deviation. Sea level at the acquisition time of TanDEM-X and TerraSAR-X was determined by combining EIGEN-GL04C geoid and the FES2014 global tide model. Finally, the bathymetry data were employed from BEDMAP2. This technique was applied to Campbell Glacier in East Antarctica and Thwaites Glacier in West Antarctica and the results suggest that continuous observation and monitoring of grounding line is feasible with this technique using multitemporal surface topography datasets. This work not only contributes to existing knowledge of previously estimated grounding line by validating its location, but will also be able to provide a new understanding of how grounding line migrates in short time of period as well as in local regions. As an example, a detailed study on the degradation of a newly detected ice rumple at the Thwaites Ice Shelf between 2011 and 2013 was carried out. This ice rumple was located nearly 5 km offshore from the estimated grounding line and the monitoring technique successfully captured that its grounding radius progressively decreased after its first site in 2011. A simple viscoelastic deformation model was used to better understand its nature assuming that the progressive disappearance of the ice rumple was caused by unpinning of the ice from the bedrock due to ongoing ice thinning. ii Nevertheless, the main weaknesses of this technique was that the accuracy of results strongly depends on the accuracy of used input datasets and ice column density. Comparison of the results with the bathymetry datasets with different spatial and vertical accuracy confirmed that the estimation could improve significantly with better datasets. In addition, the result also confirmed that the hydrostatic ice thickness could be precise enough once an ice column density, which is less well constrained around the Antarctica, is known. In spite of its weaknesses and limitations, this study certainly showed its feasibility of grounding line monitoring using hydrostatic ice thickness. Unfortunately, an issue of ice thickness errors near grounding line was not fully addressed in this study. This is another considerable limitation of using hydrostatic potential for grounding line estimation because ice is not in full hydrostatic equilibrium state in the vicinity of grounding line, called grounding zone or hinge zone. Further work will need to be done to determine how much hydrostatic thickness can deceive in this zone. Keywords: Grounding line, ice shelf, Antarctica, synthetic aperture radar, digital elevation model, radar altimeter, hydrostatic, ice rumple, viscoelastic deformation Student Number: 2010-23131 iii Table of Contents Introduction ............................................................................. １ 1.1 Study Background ....................................................................................................... １ 1.2 Purpose of Research ................................................................................................... ４ Generation of Surface Topography of Antarctica Ice Shelf using TanDEM-X.......................................................................... ５ 2.1 Background of surface topography generation in Antarctica ................... ６ 2.2 Combined usage of TanDEM-X and CryoSat-2 for generating a high resolution surface topography ............................................................................... ９ 2.3 Surface topography of Thwaites Glacier in West Antarctica..................１９ 2.4 Discussion...................................................................................................................３１ High Resolution Grounding Line Mapping Using Hydrostatic Ice Thickness .............................................................. ３７ 3.1 Calculation of hydrostatic ice thickness of floating ice shelves ............３８ 3.2 Grounding line estimation using TanDEM-X and BEDMAP2: Campbell Glacier, East Antarctica .....................................................................４２ 3.3 Grounding line estimation using improved bathymetry and additional field measurements: Thwaites Glacier, West Antarctica ...５２ 3.4 Sensitivity analysis with surface and bathymetry elevation .................６２ Local Scale Grounding Line Changes: A Case Study of an Ice Rumple in Thwaites Ice Shelf ........................................... ７０ 4.1 Description of an ice rumple and study area ...............................................７１ 4.2 Generation of high-resolution deformation maps .....................................７５ 4.3 Interpretation of surface deformation using of a viscoelastic deformation model .................................................................................................８０ 4.4 Discussion...................................................................................................................９１ Conclusions ...................................................................... １０２ Abstract in Korea ..................................................................................... １０８ Bibliography .............................................................................................. １１２ iv List of Figures Figure 2.1. TanDEM-X intensity overlaid MODIS (12/08/2010); a composite of CryoSat-2 point data during the month of June in 2011 is plotted in colors to distinguish the acquisition dates; and locations of selected ground control points are in cyan diamonds. ...................................................................１６ Figure 2.2. High resolution digital elevation model of TanDEM-X after absolute height correction with selected CryoSat-2 SARIn data; Airborne Topographic Mapper data from Operation IceBridge is plotted over the overlapping area; flight line is numbered in the order of acquisition time; the study area is separated into three sections (green dashed box). Finally, grounding lines from ASAID (Bindschadler et al., 2011), MEaSUREs(Rignot et al., 2011) and this study are drawn. ................................................................２２ Figure 2.3 Scatter plot between CryoSat-2 heights and the corrected unwrapped heights using linear least-squares approach when a) all of the available points are used as GCPs, b) GCPs are selected with the correlation coefficients only and c) an additional threshold from RMS heights is used. The fitted curve is represented in red along with the prediction bounds in dotted line and the width is the minimum in c). .............................................２４ Figure 2.4. (a) Surface profiles of TDM_CS2 (black) and ATM (dashed red) over the grounded ice sheet; the height difference is in dashed grey line; (b) Scatter plot between the height difference and (black) backscattering coefficients in dB and (red) interferometric coherence to search for any relation with radar penetration. ............................................................................２６ Figure 2.5. Surface profiles of line 8 (top) and 3 (bottom) are plotted along with the topography map (middle) for spatial comparison of grounding lines. The black, dashed red and dashed grey line represent TDM_CS2, ATM and the height difference, respectively. Inside the middle image, the locations of intersecting points between grounding lines (red: MEaSUREs, black: this study) and the ATM flight lines are extended from the map to the plots in v red and black dashed lines, respectively. The blue dots across the map represent the ATM flight tracks. .............................................................................２８ Figure 2.6. (a) TanDEM-X intensity image from June 10, 2011; (b) Landsat ETM+ panchromatic image from December 15, 2011. The surface feature inside the dotted box moved 1 km north. The yellow dots represent the location of ATM line 1; (c) The surface profiles for line 1, 9 and 2. TDM_CS2 and ATM are shown in solid black and dashed red lines, respectively. The difference between TDM_CS2 and ATM is almost plotted in dashed grey line with a difference scale on the right axis. Dashed boxes for line 9 and 2 represent the grounded area. .......................................................................................................３０ Figure 2.7 A cross-point analysis between CryoSat-2 and ATM on a monthly basis over the grounded ice sheet region between 75.35° S and 75.75° S shows a decreasing trend of surface height at a rate of -2.38 m/year. The red and blue square indicates June and November in 2011 corresponding to TanDEM-X and ATM acquisition, respectively. The relative difference between two periods is -4.94 m. ............................................................................３４ Figure 3.1 Schematic drawing of an ice/firn two layer model for hydrostatic thickness estimation. ..................................................................................................３９ Figure 3.2 Campbell Glacier in East Antarctica. Grounding lines, estimated in 1996 by DDInSAR (Rignot et al., 2011) and in 2008-2009 by MOA (Bohlander and Scambos, 2007; Scambos et al., 2007), are overlaid in the flow velocity of the glacier and its glacier tongue. An inset indicates the location of Campbell Glacier as a red dot. ..........................................................４４ Figure 3.3 Amplitude (a), interferometric coherence (b), and unwrapped phase (c) of TanDEM-X SAR data in Campbell Glacier. The amplitude and coherence image show that the ocean is covered with sea ice. .................４７ Figure 3.4 Generated DEM of Campbell Glacier and its ice tongue. Floating regions are colored differently to show the grounding line variations during 2013 and 2016. Green region indicates advance whereas red vi regions indicates retreat. Black solid lines indicate the grounding line in 1996. ..................................................................................................................................４９ Figure 3.5 Grounding lines with respect to various ice densities. Previously estimated grounding lines in 1996 and 2009 are overlaid in red and yellow lines, respectively. ........................................................................................................５１ Figure 3.6 A line segment from point A to A’ show how ice bottom varies and crosses the seafloor at different density values...............................................５２ Figure 3.7 Extracted DEMs in 2011, 2012 and 2013. CryoSat-2 SARIn data points (red) within 15 days from the TanDEM-X acquisition dates are initially chosen and rigorously selected as GCPs using correlation and RMS between the datasets. The used GCPs are shown in blue dots. ....................................５５ Figure 3.8 Retrieved tide level from FES2014 tide model in June 2011 (a), June 2012 (b), and June 2013 (c). Tide level at the time of TanDEM-X acquisition is represented as a red diamond. ...........................................................................５６ Figure 3.9 a) A flight track of MCoRDS radar sounder of Operation IceBridge, acquired on November 4, 2011. b) Profiles of MCoRDS (black) and the estimated ice column thickness (red) along the flight line. ........................５７ Figure 3.10 The estimated grounding lines in 2011, 2012 and 2013. .............６０ Figure 3.11 Eastern part of the grounding line of Thwaites Glacier (Yellow: 2011/06/10, black: 2011/10/21, blue: 2012/06/29 and red: 2013/06/27) ..............................................................................................................................................６１ Figure 3.12 Western part of the grounding line of Thwaites Glacier (Yellow: 2011/06/10, black: 2011/10/21, blue: 2012/06/29 and red: 2013/06/27) ..............................................................................................................................................６２ Figure 3.13 Grounding line variation with respect to varying input surface elevations. The reference grounding line is estimated using a fixed ice column density of 880 kg/m3. The darker lines represent grounding line retreat as decreasing surface elevation...............................................................６４ Figure 3.14 Grounding line variation with respect to varying input bathymetry elevations. The reference grounding line is estimated using a fixed ice vii column density of 880 kg/m3. The darker lines represent grounding line advance as increasing bathymetry elevation. The flow vectors are also overlaid in green arrows. ..........................................................................................６５ Figure 3.15 Grounding line using AirGrav (white), BEDMAP2 (purple), and DDInSAR (black) are plotted over AirGrav bathymetry map. Hydrostatic anomaly using an airborne radar sounder data indicates where ice is floating (red) and grounded (blue).......................................................................６７ Figure 4.1 Grounding lines of Thwaites Glacier in 1996 (green) and 2011 (red) estimated using the DDInSAR method with European Remote Sensing (ERS) satellites (Rignot et al., 2014). The orange dotted rectangle indicates a larger ice rumple previously discussed by Tinto and Bell (2011). A newly generated digital elevation model derived from TanDEM-X data is shown within the diagonal rectangle, in which the small red feature inside the yellow dotted square indicates the smaller ice rumple considered in this study...................................................................................................................................７４ Figure 4.2 Generated high resolution surface topography maps of the Thwaites Ice Shelf on (a) 2011/06/10, (b) 2012/06/29, and (c) 2013/06/24. Dark grey denote the locations of all of the available CryoSat2 SARIn altimeter elevation data and red circles denote selected data. ７８ Figure 4.3. Detailed surface DEMs of the ice rumple from 2011–2013, its yearly surface deformation, and surface velocity (arrows): (a-c) Location of maximum elevation in each of the three years (filled circles indicate the peak in that year); (d-e) Deformation maps for the periods 2011–2012 and 2012–2013 (later year subtracted from earlier year) where cyan circles indicate maximum surface depression and black circles indicate maximum elevation in each reference period. .......................................................................８０ Figure 4.4. Schematic drawings of ice rumple (a) and the elastic deformation model with a finite spherical pressure source (b). The initial surface is not flat, thus the radial distance from the center of a sphere to the surface viii varies. ................................................................................................................................８１ Figure 4.5. Comparison between the observed and estimated deformation from the modified deformation model after the global optimization process during 2011–2012. The RMS difference between two deformation maps was 1.53 m. Profiles in the longitudinal and along-flow directions are plotted along the red dashed and dot-dashed lines, respectively.............８７ Figure 4.6. Comparison between the observed and estimated deformation after the global optimization process during 2012-2013. The RMS difference between two deformation maps was 5.55 m. Profiles in the longitudinal and along-flow directions are plotted along the red dashed and dot-dashed lines, respectively. ........................................................................................................８９ Figure 4.7. Comparison between the observed and estimated deformation from the modified deformation model after the global optimization process during 2012–2013. After the estimated deformation time was set to the measured time (0.99 years), the RMS difference decreased from 5.55 m to 2.11 m................................................................................................................................９０ Figure 4.8. Landsat 7 ETM+ images from 2003–2014 showing the gradual disappearance of the studied ice rumple in the Thwaites Ice Shelf. Crevasses and surface gradients are generally created atop an ice rumple due to surface extension and elevation increase. Such features were visible as late as 2011 but disappeared by 2013, indicating gradual ice thinning. Larger images are magnifications of selected areas indicated by red boxes. ..............................................................................................................................................９２ Figure 4.9 The viscoelastic deformation model is fitted to the observed deformation (a). The modelled deformation (c) is a summation of viscous (b) and elastic deformation......................................................................................９４ Figure 4.10. Schematic drawing of surface deformation in consecutive periods and corresponding ice thicknesses, defined as the difference between depth and radius of a spherical pressure source. ...........................................９６ Figure 4.11 Elevation and deformation profiles over the ice rumple and ix surrounding ice shelf. .................................................................................................９８ Figure 4.12. The contribution of variations of the input parameters to the output of the surface deformation model was investigated by changing the variables within 10% from the estimated values. The center location of a spherical source showed the least sensitivity, whereas the radius made the highest contribution............................................................................................... １０１ x List of Tables Table 2.1. Specification of datasets .................................................................................３５ Table 2.2. Bias and standard deviation of height differences between TDM_CS2 and ATM (TDM_CS2 minus ATM) for each line. The lines are listed in the order to their locations from the north. Regional biases and standard deviations are also calculated. ................................................................................３６ Table 4.1. Minimum, maximum, and optimized values of each unknown variable for constrained global optimization of the viscoelastic surface deformation. ..............................................................................................................................................８６ xi Introduction 1.1 Study Background Grounding line is where ice body extends into an ocean and is detached from the underlying bed to become a floating ice shelf. Its dynamics is a crucial boundary condition for ice-sheet modelling as well as an indicator of ice-shelf stability. Monitoring changes in grounding line positions has substantial importance of determining vulnerability of ice shelves, when it happens. In West Antarctica, the grounding line retreat has been the most conspicuous indicator of negative mass balance in recent years. The bathymetry in West Antarctica has a negative slope landward and the grounding line retreat in this region especially can cause high vulnerability of stability of ice shelves as its ocean base is exposed to relatively warmer sea water, Circumpolar Deep Water (Mercer, 1978). As a consequence, various researches have tried to reveal dramatic changes such as rapid flow velocity and higher surface thinning rate in the Amundsen Sea Sector of West Antarctica (Rignot, 2006; Wingham et al., 2006). These changes are closely related to the grounding line retreat. There are three well-known methods to define grounding line positions using satellite observations. First of all, grounding line position is defined and monitored using the surface elevation changes due to ocean tide. Repeat-pass laser altimetry using the Ice, Cloud, and land Elevation Satellite (ICESat) is utilized to detect and monitor relative tidal displacement of ice shelf and to １ define the grounding line which lies at the landward edge of a zone where ice shelf flexure occurs (Fricker and Padman, 2006; Fricker et al., 2009; Brunt et al., 2010). Grounding line monitoring by tidal deflection of a floating ice tongue is also possible using Double-Differential InSAR technique (Rignot et al., 2011). DInSAR technique can provide the grounding location with finer spatial resolution, but the spatial and temporal extent of available data is relatively poor. In other words, the availability of suitable datasets, tidal ranges at the time of satellite data acquisitions, and environmental conditions of ice shelves such as flow velocity and surface roughness can limit areas where this technique can be applied. In addition, the monitoring using laser altimetry was only available for 6 years during its mission lifetime and the ground track coverage is relatively sparse. Secondly, grounding line position is defined in the vicinity of the break in ice surface slope where shadows in optical satellite image occur when favorably illuminated (Scambos et al., 2007). Bindschadler et al. (2011) utilized photoclimetry technique to generate a DEM using Landsat images and extracted the location of grounding line by hands as the Antarctic Surface Accumulation and Ice Discharge (ASAID) project. This technique provided the grounding line of the entire Antarctica but the accuracy was relatively low. Lastly, grounding line position is defined by the hydrostatic potential where floating ice is in a hydrostatic equilibrium state and grounded ice is not (Fricker et al., 2002). The surface elevation of ice shelf can be converted to ice thickness using a simple relationship between the density difference between the ice and sea water. Then, the extent of the floating ice over an ice mass can be defined ２ if the surface height, ice thickness, and densities are accurately known. This is a robust technique of monitoring grounding line positions because it physically defines where the ice starts to depart from the bed. However, the greatest concern is that it strongly depends on the accuracy of the input variables such as densities and surface elevation. The most important aspect of this study is the potential of monitoring grounding line migration in high spatial and temporal resolution around Antarctica by using a simple and repeatable estimation technique. In this study, high resolution digital elevation models are generated and utilized to estimate, update, and monitor grounding lines of Antarctic ice shelves. Then, grounding line estimation is done by defining a point at which the bathymetry meets the bottom of ice shelves of which is calculated by subtracting thickness from surface. This technique is similar to the hydrostatic potential technique, but it uses the bathymetry instead of ice thickness because accurate ice thickness data is not available around the Antarctica. Ice shelf thickness is estimated from freeboard assuming a hydrostatic equilibrium state. The primary factor in calculation of bottom topography of ice shelf is an accuracy of freeboard height, a distance between ice shelf surface and sea level and ice column density. In Chapter 2, a new technique is introduced to have precise surface elevation of ice shelves using TanDEM-X bi-static interferometry. The absolute height offset, which has been a persistent problem in applying the interferometry technique for generating DEMs, is solved by providing reliable reference heights from CryoSat-2. In Chapter 3, the description of calculating hydrostatic ice thickness is provided and grounding lines of Campbell Glacier, East ３ Antarctica and Thwaites Glacier, West Antarctica are estimated with BEDMAP2 as the reference bathymetry. The grounding line of Thwaites Glacier is updated with a new regional bathymetry data for comparison. In Chapter 4, monitoring of a localized grounding line evolution is described. Then, the thesis will close with the further studies for the following step in Chapter 5. 1.2 Purpose of Research While DDInSAR has been the most widely used technique, the spatial and temporal extent for continuously monitoring long-term and short-term variation of grounding line position is still needed. It was decided to exploit surface topography of ice to estimate and extract grounding line information. Thus, the main purpose of this research is to develop an effective technique 1. To generate the most accurate surface topography of ice using TanDEM-X and TerraSAR-X bi-static SAR interferometry data, 2. To estimate locations of grounding line using generated surface topography and hydrostatic potential, 3. To monitor the spatio-temporal migration of the grounding line by using multi-temporal datasets. ４ Generation of Surface Topography of Antarctica Ice Shelf using TanDEM-X “Combined Usage of TanDEM-X and CryoSat-2 for Generating a High Resolution Digital Elevation Model of Fast Moving Ice” Definite surface topography of ice provides fundamental information for most glaciologists to study climate change. However, the topography at the marginal region of ice sheets exhibits noticeable dynamical changes from fast flow velocity and large thinning rates; and thus, it is difficult to determine instantaneous topography. In this study, the surface topography of the marginal region of Thwaites Glacier in the Amundsen Sector of West Antarctica, where ice melting and thinning are prevailing, is extracted using TanDEM-X interferometry in combination with data from the near-coincident CryoSat-2 radar altimeter. The absolute height offset, which has been a persistent problem in applying the interferometry technique for generating DEMs, is determined by linear least-squares fitting between the uncorrected TanDEM-X heights and reliable reference heights from CryoSat-2. The reliable heights are rigorously selected at locations of high normalized cross-correlation and low RMS heights between segments of data points. The generated digital elevation model with the resolved absolute height offset is assessed with airborne laser altimeter data from the Operation IceBridge which was acquired 5 months after TanDEM-X ５ and shows high correlation with biases of 3.19 m and -4.31 m at the grounding zone and over the ice sheet surface, respectively. Finally, it is expected that the combination of interferometry and altimetery with similar datasets can be applied at regions even with a lack of ground control points. 2.1 Background of surface topography generation in Antarctica The importance of observing the West Antarctic Ice Sheet has been increasing following a number of alarming scientific results of its potential vulnerability due to its exposure to relatively warmer sea water (Jenkins et al., 2010; Joughin et al., 2014; Rignot et al., 2014). In particular, observation of surface topography change can be important information for glacial remote sensing research and most glaciologists to study ice mass balance and its contribution to sea level. Also, the marginal region requires continuous surface observation to understand the significant dynamic changes such as faster flow velocity and larger thinning rates (Rignot et al., 2013). In Antarctica, the surface topography of glaciers and ice sheets is measured by satellite radar and laser altimeters, such as ERS-1/2, Envisat RA-2 and ICESat. These satellite altimeters have measured the surface height data from the early 1990s to present. The typical methods to utilize the accumulated altimeter data to extract the surface height variation are cross-track altimetry and along-track repeat altimetry (Davis and Ferguson, 2004; Legresy et al., 2006). These methods have limitations that may be potentially addressed using other methods such as radar interferometry. For example, studies that used radar ６ altimeters rejected elevation data in the marginal regions of ice sheets where rough surfaces with crevasses and rapid slope variations were present because satellite radar altimeters have a tendency to lose their accuracy in regions of complex topography (Shepherd and Wingham, 2007). Laser altimeters have greater vertical precision with smaller footprint sizes than radar altimeter but large gaps between each track must be filled by interpolation. The GLAS/ICESat 500m Laser Altimetry Digital Elevation Model of Antarctica is an example of an interpolated grid from laser point measurements and has a spatial spacing of 500m and 1km (DiMarzio et al., 2007). Besides, the DEM, especially in the marginal regions of glaciers with fast flow velocity and high thinning rate, might contain errors caused by the surface height variations during the data acquisition between 2003 and 2005. Thus, space-borne altimetry has mostly been effective in terms of measuring relative variation of the polar surface. Synthetic aperture radar interferometry is another effective tool for measuring surface topography for an extended area. Generation of surface topography using interferometry is not a new approach for the Antarctic ice bodies. However, the main limitation of radar interferometry for ice sheets and glaciers has been the temporal decorrelation effect seen from low interferometric coherence. Low interferometric coherence can decrease the accuracy of generated surface topography. Thus, special tandem missions to reduce the repeat cycle between data acquisitions had been carried out a number of times. For example, COSMO-SkyMed (Constellation of small Satellites for the Mediterranean basin Observation), conducted by the Italian Space Agency ７ (ASI), is a constellation system using four satellites to reduce the repeat cycle to a single day to minimize the temporal decorrelation. TanDEM-X (TerraSARX add-on for Digital Elevation Measurement) is a twin satellite of TerraSAR-X that almost completely eradicates the temporal difference between data acquisition (Krieger et al., 2007; Krieger et al., 2013). Another major limitation of achieving the accurate elevation grid in high resolution in the marginal region of Antarctica is the lack of coincident ground reference at the time of SAR acquisition due to its continuously moving state. In order to generate precise digital elevation model (DEM) using radar interferometry, accurate ground truth data known as ground control points (GCP) are necessary as reference. It is difficult to obtain these reference data over continuously flowing and changing ice and snow surfaces. A previous study by Baek et al. has combined ERS tandem mission data and ICESat laser altimetry to generate a DEM over King Edward VII Peninsula in West Antarctica (Baek et al., 2005). However, the time gap between ERS and ICESat was 8 years, thus the bias had to be compensated by assuming average thinning rate of the region. Although this technique provided the preliminary reference of the surface at the time, it can hardly be used when measuring the surface height variations and corresponding volume changes over time. Thus, in this study, we overcome the mentioned limitations by combined usage of single-pass synthetic aperture radar interferometry (TanDEMX/TerraSAR-X, denoted as TDM) and near-coincident radar altimetry (CryoSat-2, denoted as CS-2) for the first time. The basic idea of combining interferometry and altimetry is not completely new as in other studies such as ８ Zhou et al. (Zhou et al., 2015) in which utilized ICESat and 35-day repeat-pass interferometry to compensate the ice motion and atmospheric phase; singlepass bi-static configuration of TanDEM-X/TerraSAR-X easily solves the mentioned interferometric problem. Besides, Groh et al. (2014) presented DEM analysis over Thwaites Glacier with near-coincident lidar altimeter; however, airborne campaign at every data acquisition is impossible. Therefore, CryoSat2 in SAR Interferometric mode (SARIn) point measurements are used as ground control points (GCPs) to calculate absolute height offset. CryoSat-2 is the only available radar altimeter with an ability to provide near-coincident surface height measurements along with TanDEM-X/TerraSAR-X as discussed later in Section 2.2. The accuracy assessment is carried out with the Airborne Topographic Mapper (ATM) from the Operation IceBridge by NASA, acquired 5 months after the SAR acquisition. Finally, the generated high resolution DEM is used to estimate grounding line location of the ice shelf by utilizing hydrostatic equilibrium theorem. Section 2.2 gives detailed information about generating a DEM, datasets and study area. Section 2.3 shows a result and its assessment for a case study of Thwaites Glacier in West Antarctica. 2.2 Combined usage of TanDEM-X and CryoSat-2 for generating a high resolution surface topography The main objective of TanDEM-X mission (TDM) is to build the most accurate DEM for the globe with the estimated relative vertical accuracy of about 2m and the horizontal accuracy of 3m (Krieger et al., 2007; Krieger et al., 2013). Among several operational modes, TanDEM-X can be operated in a bi９ static mode with which TanDEM-X transmits the microwave signal and both TerraSAR-X and TanDEM-X simultaneously receive the backscattered signal. This operational mode is practically appropriate in this particular study due to the following reason. The typical drawbacks of SAR interferometric DEM over ice and snow have been the time difference between acquisitions and atmospheric artifacts. The first effect creates significant phase decorrelation that also adversely affects in the phase-unwrapping process. The second effect is observable in the Polar Regions where extreme weather conditions are common and change rapidly, although atmospheric phase delay is less critical than in low latitudes. These limitations are solved by simultaneous acquisitions of TanDEM-X and TerraSAR-X. A DEM is generated by calculating an interferogram with two SLC SAR scenes. Interferogram are multi-looked with various factors depending on wanted spatial resolution in azimuth and range directions. In this study, an interferogram is multi-looked by a factor of 2, resulting 3.89 m and 4.18 m of ground resolution in range and azimuth directions, respectively. Since two TanDEM-X Single Look Complex SAR images are provided as Co-registered Single-look Slant-range Complex (CoSSC) images, general co-registration step is not required for calculating an interferogram (Rossi et al., 2012). Then, the phase corresponding to Earth’s curvature is removed from the calculated interferogram with respect to the WGS84 ellipsoid. The Earth-curvature flattened interferogram is filtered using the Goldstein filter (Goldstein and Werner, 1998) and interferometric coherence is estimated to mask out regions １０ with low coherence. Because the residual phase still has modulo of 2 a phase unwrapping procedure is required in order to obtain an absolute phase corresponding to the ice topography. In this study, the Minimum Cost Flow method (Costantini, 1998) is applied for phase unwrapping with the coherence of greater than 0.40. After interferometric phase is unwrapped, there still exists a constant phase offset to represent the true topography, given by ztrue    [unw  off ] , (2.1) where  unw is unwrapped phase,  off is an absolute phase offset, ztrue is true topographic height represented by interferometric phase and  is the coefficients for converting phase to height (Perna et al., 2015).  is a function of geometric parameters of SAR acquisition as given by     R  sin( ) , 2    B (2.2) where  is a wavelength, R is slant range distance,  is look angle and B is perpendicular baseline (Perna et al., 2015). The height of ambiguity can be calculated by multiplying 2to  . In general, an absolute phase offset is determined using external reference topographic heights, known as ground control point (GCP), within a scene. The process takes the ground range coordinates and corresponding topographic height of GCPs are converted to slant range coordinate and corresponding unwrapped phase using SAR geometry; then a constant phase offset is determined by minimizing difference between converted phase of GCPs and １１ previously unwrapped phase of an interferogram. In addition, GCPs are then used to minimize the phase ramp caused by small baseline inaccuracies. Nonetheless, finding a constant absolute phase can be affected by a number of uncertainties, such as the accuracies of GCP positions, topographic changes between the periods of data acquisitions and systematic noise in the unwrapped interferograms. The error can induce topographic height error proportional to the geometric coefficient  such as where z gcp gcp z gcp    [unw  off ], (2.3) gcp z  z gcp  ztrue    [off  off ]     , (2.4) is the topographic height converted with an erroneous absolute gcp phase offset,  off , and  is the phase error as mentioned above. The topographic height error varies in the range direction with respect to the range, look angle and perpendicular baseline. The bigger problem with the height error induced by the erroneous phase offset is that the vertical accuracy of the generated DEM significantly declines when the GCPs are vertically biased (Perna et al., 2015). The GCPs can be vertically biased when surface height variation has occurred since the acquisition of the GCPs and a continuous measurement is not available. This can be an exceptional problem especially for the applications over the ice sheet surface because of continuous height variations as time goes by. Perna et al. (2015) tried to solve the problem by executing the initial phase offset estimation in phase domain and additionally by refining the phase error １２ using an iterative height estimation approach in topography domain. The method verified that generation of an accurate DEM was practical even without precise reference points such as corner reflectors. However, areas with high slopes have to be masked in order to retain its robustness because of serious planar misalignment. Although masking out high-sloped areas is a suitable solution, it can cause a notable complication if the study area is the marginal region of ice sheet where surface slope is relatively high. Therefore, in this study, we consider that selection of reliable GCPs must be carried out cautiously instead of masking out the entire region with high slopes. In this study, the initial phase offset is assumed to be zero such that let zunw be the topographic height derived directly from unwrapped phase expressed by zunw    unw . (2.5) zest  A  zunw  B , (2.6) Then, we assume that where zest is the estimated topographic height with A and B being two linear least squares estimators (LSE) to satisfy that N f ( zunw)  [ zest  z gcp ]2 , (2.7) 1 is the minimum where z gcp is the topographic height of N ground control points. In this way, the method is slightly bypassed by directly determining unbiased １３ estimators (A and B) to minimize height variance from reference heights using a linear least squares approach. They also demonstrate that it is possible to calculate absolute topography without accurate estimation of an initial phase offset (Perna et al., 2015). Due to the fact that the glacier topography is constantly changing, it is quite difficult to obtain reliable GCPs that are coincident with the acquisition of SAR data. This becomes even a bigger problem over the fast flowing marginal region of ice sheet because installation of any ground-based references, such as corner reflectors and GPS, may be lost due to heavy snow or horizontal displacement. In addition, the number of GCPs needs to be large and the points are to be spatially distributed throughout the image to satisfy the quality of the absolute height estimation. In order to solve the supply problem of GCPs, the usage of CryoSat-2 radar altimeter is introduced. After its launch in April, 2010, CryoSat-2 (CS-2) has been continuously providing precise surface elevation data over the Polar Region along with TanDEM-X for which was also launched in 2010. Considering 92 degrees of inclination angle, a 369 days repeat period and a smaller effective footprint of CS-2, reliable and well distributed surface elevation points can be obtained even within a single SAR image (Wang et al., 2015). Besides, unlike any other altimeters, CS-2 provides a location of the returning signal in across-track direction within 1.6 km from its nadir, mostly at the top of surface undulations of rough ice shelf (Wingham, 2005). Thus, in this study, the Level 2 product of CS-2 SARIn data, processed according to the ESA Baseline C (ftp://science-pds.cryosat.esa.int), are １４ considered as good candidates of GCPs for determining the absolute height. In order to secure GCPs well distributed within a scene, a composite of CryoSat2 tracks during the period of 1 month, centered on the SAR acquisition date, is collected. There are usually 5-6 tracks crossing a single SAR scene as shown in Figure 2.1. １５ Figure 2.1. TanDEM-X intensity overlaid MODIS (12/08/2010); a composite of CryoSat-2 point data during the month of June in 2011 is plotted in colors to distinguish the acquisition dates; and locations of selected ground control points are in cyan diamonds. １６ As mentioned above, the absolute height estimation is only robust when the GCPs are not located in the areas with high slope. Besides, GCPs may represent wrong topography due to temporal variation of the surface. Thus, a selection process is carried out with an introduction of two criteria to utilize only reliable CS-2 data as GCPs. First of all, CS-2 Level 2 data are provided with quality flags so that the points with large orbital, height and angle errors can be easily identified. Also, due to large slope and rough surfaces, problematic unwrapping could cause pointing errors leading to dislocation of the point-of-closest-approach (POCA) and significant height. Thus, the points with these systematic quality flags are disregarded. As a second selection criterion, CS-2 points are chosen based on a correlation-based filter with two thresholds. The first threshold is normalized local pixel cross-correlation between a segment of CS-2 and corresponding uncorrected topographic height, zunw , given by  ( z n CorrNorm  i 1 unw n  z i 1 (i )  zunw)( z gcp (i )  z gcp ) (i )  zunw  2 unw  z n i 1 gcp  (i )  z gcp  , (2.8) 2 where n is the number of points within a segment and overbars are mean heights of each segment. Result not only represents height similarity, but also suggests spatial correlation between the data in a sense of digital image correlation technique because CS-2 SARIn data points are spatially distributed, displaced from the １７ sub-satellite ground track. In other words, cross-correlation of a segment of points can also examine location accuracy of CS-2 data. For example, location of the POCA is determined by comparing the differential phases of the returned signal between two antennas and phase-wrapping can occur when the across track slope is great enough. This causes the POCA to appear on the other side of the nadir ground track; and thus, contains wrong topographic height. As a result, a segment with this erroneous positioning will display a low correlation coefficient and will be discarded. Although the CS-2 elevation accuracy is relatively high, GCPs determined from CS-2 at a given time cannot be used for InSAR data that were acquired at a different time. Due to the fact that CS-2 data are collected only nearcoincident, the topography may have changed significantly even during one month. For example, if the surface thinning rate is around 3~4 m per year, which can be observed in the marginal region of the West Antarctic Ice Sheet, the vertical error of the GCPs can be nearly 0.3 m. The vertical error can be even greater if the thinning rate is not constant. Also, as topographic changes may be spatially different, we anticipated that such CS-2 data would exhibit low correlation coefficients. Along with the topographic variation between data acquisitions, the correlation coefficient determines the effect of dissimilarity between radar altimetric and interferometric penetration depths in similar fashion. To distinguish the effect of different penetration depths from the topographic variation is difficult without additional external ground reference information such as corner reflectors. However, if there is a spatially varying penetration １８ depth difference, the coefficient between segments will be low. The second threshold is RMS height of a segment of CS-2 given by RMS CS 2  1 n [ z gcp (i)  z gcp ]2 ,  n i 1 (2.9) where n is the number of points of a segment. Low RMS height of a segment illustrates a relatively plane area with low slope expecting high accuracy of CS2; thus, CS-2 points with high correlation coefficients and low RMS height are extracted as reliable GCPs. 2.3 Surface topography of Thwaites Glacier in West Antarctica The region of interest of this study is the marginal region of Thwaites Glacier, located along the Amundsen Sea in West Antarctica (Figure 2.1). Since Mercer (1978) addressed the issue, the WAIS is considered vulnerable to collapse because its ocean base is exposed to relatively warmer sea water, Circumpolar Deep Water, and the base of these vast glaciers are located below sea level. Thwaites Glacier is one of these glaciers in West Antarctica and is well known for significant negative mass changes and a great deal of crevasses and steep slopes over fast flowing surface, nearly 4 km per year at the fastest (Rignot, 2001; Rignot and Jacobs, 2002; Rignot et al., 2002; Rignot et al., 2011; MacGregor et al., 2012; Joughin et al., 2014; Rignot et al., 2014). In addition to its rapid grounding line retreat, surface thinning, and decreasing extent of floating ice, the region is known for strong katabatic wind and rapid and heavy precipitation with which the surface condition can change dramatically at any time. １９ A pair of TanDEM-X and TerraSAR-X SAR Single Look Complex (SLC) data was acquired in June 10, 2011 in a descending orbit. The system uses Xband radar with center frequency of 9.6 GHz (Table 2.1). TanDEM-X and TerraSAR-X SLC data pair was acquired over the marginal region of Thwaites Glacier in HH polarization with a mean incidence angle of 44 degrees on June 10, 2011. The time difference between two SAR image acquisitions was about 7.99x10-6 second, which is considered to be simultaneous. The range and azimuth resolution of two SAR images in strip map mode were about 1.68 m and 3.30 m, respectively. It has been proven that radar altimeter can provide reliable measurement over snow and ice surface in Arctic and Antarctica (Legresy et al., 2005; Rémy and Parouty, 2009). However, some problem occurred in the regions where terrain is rough and slope is steep. The problem mainly comes from the large footprint size of typical radar altimeter systems. In order to overcome the problem, CryoSat-2 was developed to be able to collect the size of a pulse-limited footprint of 1.67 km in across-track and 313 m in along-track direction. Its high-spatial resolution has been implemented using the Doppler beam formation which allows to record the returned signal in along and across direction separately (Raney, 1998; Wingham et al., 2006). Besides, CryoSat-2 in Synthetic Aperture Interferometric mode (SARIn) operates with two receiving antennas to locate the POCA in the range direction diminishing position errors (Wingham, 2005). CryoSat-2 utilizes Ku-band with the center frequency of 13.375 GHz with 1.17m of the baseline between two antennas (Table 2.1). The height accuracy ２０ of CryoSat-2 SARIn over Antarctic ice sheet is reported to have 1.5m where surface slope is less than 0.2 degrees (Wang et al., 2015). The airborne laser altimeter (IceBridge ATM L2 Icessn Elevation, Slope, and Roughness) data from Operation IceBridge from NASA, obtained on November 4th, 2011 (Blankenship et al., 2012, updated 2013), were used as reference. It provided the location and slopes of each smoothed laser point block of which sampling frequency was 5 kHz. The airborne laser altimeter has a reported accuracy of 10 cm (Krabill et al., 2002) and it is known that the penetration depth of laser signal on snow and ice surface is negligible, and thus represents the topography of the top layer of the ice sheet surface. Figure 2.2 illustrates the generated high resolution DEM using a composite of monthly CS-2 data as reliable GCPs after a selection process. ENVI SARScape is used to generate and unwrap interferograms of TanDEM-X data. After converting it to the uncorrected height and geocoding, the absolute height estimation as well as GCP selection is done using Matlab. GCPs are selected when correlation coefficient is greater than 0.95 with RMS less than 10 m. The number of points within a segment is 5 not only to maintain enough sample number of GCPs, but also to consider that a segment makes azimuth length of around 1.5 km which is nearly the length of azimuth footprint. More than 1/3 of GCP candidates have higher than 0.95 correlation coefficients where high correlation was originally expected owing to high accuracy of CryoSat-2 data. Besides, RMS height plays a significant role in picking out bad candidates with large biases, dramatically decreasing the number of GCP candidates. Finally, a selected RMS height from the center of ２１ 1.5 km represents a slope of 0.76 degrees, indicating selected candidates are from relatively flatter regions. Figure 2.2. High resolution digital elevation model of TanDEM-X after absolute height correction with selected CryoSat-2 SARIn data; Airborne Topographic Mapper data from Operation IceBridge is plotted over the overlapping area; flight line is numbered in the order of acquisition time; the study area is separated into three sections (green dashed box). Finally, grounding lines from ASAID (Bindschadler et al., 2011), MEaSUREs(Rignot et al., 2011) and this study are drawn. Out of 368 GCP candidates, a total of 44 points are qualified under the criteria with the mean correlation coefficient of 0.97. They are mostly located over a ２２ grounded area of ice sheet; and a couple of points from floating area are seemed to be from relatively smoother area (Figure 2.1). The sampling space of the generated DEM is 5 m by 5 m, precise enough to show crevasses and large scale surface undulations over the ice shelf. Figure 2.3 shows the scatter plots between the CS-2 heights, used as GCPs, and the corrected unwrapped heights using linear least-squares approach. The fitted curve of the unbiased estimators is shown in red along with the prediction bounds in dotted red. For each figure, the number of GCPs used in the estimation differs such that in Figure 2.3a, all of the available CS-2 heights were used as GCPs and the total number was 334. The R-square value of the fitted result was 0.98; however the RMS error was 18.57 m. Figure 2.3b and 2.3c represents the cases when the correlation coefficient was the sole threshold and the coefficients and RMS heights both are combined, respectively. The width of the prediction bounds is the minimum when the correlation coefficients and RMS heights are used with the RMS error of 5.81 m. The mean difference between the CS-2 heights and the corrected unwrapped heights is nearly zero in all cases. This shows that the inclusion of the RMS height as the second threshold practically increases the accuracy of the absolute height estimation although the number of GCPs decreases. The generated DEM is directly compared with 9 ATM flight lines intersecting in longitudinal direction, west-to-east and east-to-west. Figure 2.2 shows locations of each flight line, numbered in the order of its acquisition time, ２３ Figure 2.3 Scatter plot between CryoSat-2 heights and the corrected unwrapped heights using linear least-squares approach when a) all of the available points are used as GCPs, b) GCPs are selected with the correlation coefficients only and c) an additional threshold from RMS heights is used. The fitted curve is represented in red along with the prediction bounds in dotted line and the width is the minimum in c). ２４ shown on the right side of each box. Originally provided with nadir profile plus 3 off-nadir profiles per time stamp, surface elevation acquired only from nadir direction of ATM flight lines are extracted for a single height profile with a spatial resolution 80 m; and the profile is re-processed to provide an elevation point at every 283 m, which is close to the along-track footprint size of CryoSat2. Since altimeter points are not located in the middle of each corresponding pixel of the DEM, the surface elevation of the generated DEM is bi-linearly interpolated to match exact locations of the ATM measurements. Table 2.2 shows the mean height differences between TDM_CS2 and ATM (TDM_CS2 minus ATM) of each flight line as well as of 3 sub-sites of the study area; they are 1) ice sheet, 2) grounding zone and 3) ice shelf (ice tongue) with respect to the grounding line from the MEaSUREs project. The first site contains the lines 7, 4, 6 and 5 (in order from North to South), representing the grounded ice sheet (Figure 2.4a). The profiles with an odd flight number were flown in east-to-west direction and are re-projected to represent topographic profiles in west-to-east direction for easier visualization. On average, TDM_CS2 shows lower elevation (negative bias) suggesting 4.31 m of the radar penetration if the surface had not changed during 5 months. The profiles of TDM_CS2 and ATM are also well correlated with some noticeable local height differences. ２５ Figure 2.4. (a) Surface profiles of TDM_CS2 (black) and ATM (dashed red) over the grounded ice sheet; the height difference is in dashed grey line; (b) Scatter plot between the height difference and (black) backscattering coefficients in dB and (red) interferometric coherence to search for any relation with radar penetration. ２６ Higher penetration is usually expected with lower surface density, which can be associated with fresh snow or dryer conditions. In SAR images, it can be observed as low SAR backscattering coefficients. Besides, low backscattering also leads to low coherence decreasing accuracy of interferometry. Figure 2.4b illustrates a scatter plot of the height difference versus (black) backscattering coefficients in dB and (red) interferometric coherence. More values are plotted in the negative sign as the mean difference is negative; and yet, the plot does not show any trend of height difference. The second sub-site contains line 8 and 3 at where both profiles cross the grounding line in the western and eastern side and the floating ice tongue inbetween (Figure 2.5). TDM_CS2 (red) and ATM (black) are well correlated with the mean bias of 2.19 m and 4.05 m and standard deviation of 4.31 m and 4.35 m for line 8 and 3, respectively. ２７ Figure 2.5. Surface profiles of line 8 (top) and 3 (bottom) are plotted along with the topography map (middle) for spatial comparison of grounding lines. The black, dashed red and dashed grey line represent TDM_CS2, ATM and the height difference, respectively. Inside the middle image, the locations of intersecting points between grounding lines (red: MEaSUREs, black: this study) and the ATM flight lines are extended from the map to the plots in red and black dashed lines, respectively. The blue dots across the map represent the ATM flight tracks. ２８ The positive bias of TDM_CS2 may suggest surface lowering during the 5 months of period. The thinning corresponds with other studies (Helm et al., 2014; Paolo et al., 2015) and sounds plausible; nevertheless, small bias and standard deviation infers that TanDEM-X DEM successfully extract longwavelength surface elevation near the grounding line where elevation change is the highest. At last, the third sub-site contains lines 1, 9 and 2, located over the floating ice tongue and the edge of the estimated grounding line of Thwaites Glacier (Figure 2.6). Line 1 crosses a region of disintegrating ice tongue at which was reported as the division between the western and eastern ice tongue (Kim et al., 2015). The mean height difference is the largest compared to the other flight lines, mainly due to continuous movement and calving events during the period of 5 months. Figure 2.6 illustrates such surface changes with a similar time span of 6 months with (a) an intensity of TanDEM-X and (b) a Landsat ETM+ image, respectively acquired on June 10, 2011, and December 15, 2011. As indicated in a black rectangle box inside the images, the surface feature in the box migrated about 1 km north in 6 months inferring a flow velocity of around 2 km a-1. This flow velocity is roughly close to the velocity of Thwaites Glacier around its grounding line of 3 km a-1, estimated using TanDEM-X in 2011 (Groh et al., 2014). The comparison confirms that rapid variation of surface feature is possible during a short period of time and so is the large height difference between ATM and TDM_CS2 at such an area. ２９ Figure 2.6. (a) TanDEM-X intensity image from June 10, 2011; (b) Landsat ETM+ panchromatic image from December 15, 2011. The surface feature inside the dotted box moved 1 km north. The yellow dots represent the location of ATM line 1; (c) The surface profiles for line 1, 9 and 2. TDM_CS2 and ATM are shown in solid black and dashed red lines, respectively. The difference between TDM_CS2 and ATM is almost plotted in dashed grey line with a difference scale on the right axis. Dashed boxes for line 9 and 2 represent the grounded area. ３０ Besides, the height differences may be caused by tidal variation. However, when extracting the range of tide in this area using the FES2012 tide model, the range was up to about 1 m. Therefore, the height difference between ATM and TDM_CS2, larger than 1 m, seems to be derived from the time difference between two data. Line 9 and 2 show two distinct characteristics along the tracks before and after passing the previously predicted grounding line of (Rignot et al., 2011) (Figure 2.6a). In particular, the profiles of line 2 in a dashed box are almost identical suggesting the fully or partially grounded surface. Surface profiles outside a dashed box in Figure 2.6 are located outward from the grounding line, up to 14.25 km from the west in line 9 and up to 19.33 km in line 2, experiencing rapid variation of the ice tongue during the 5 months. Line 9 represents a relatively small bias and at the same time shows a very large standard deviation. This phenomenon is due to the horizontal movement of the ice surface for 5 months and the difference of these surfaces can be displayed in the height difference plot above which changes periodically. Therefore, a mean may become small while its standard deviation becomes high. 2.4 Discussion In this study, the absolute height is determined using selected ground control points after the unwrapped phase is converted to the unwrapped height and geocoded. The estimation of zero-phase-offset is proved valid as the coefficient A (Equation 2.6) is very close to 1 otherwise the coefficient A must be other than 1 if there were any phase ramp; and this also shows that the phase ramp is ３１ negligible in this case. We decided not to use other available external DEMs in phase unwrapping, estimating the absolute phase offset or minimizing the phase ramp. We had tried to use the available DEMs, such as ASTER and GLAS/ICESat DEMs, over this region; however, the usage of these DEMs over the marginal region of Thwaites Glacier created additional phase ramps due to very poor spatial resolution and relatively extreme surface deformation since the generation of these DEMs. Besides, DEMs are not available over the floating regions. In spite of high correlation along with small bias and standard deviation between TDM_CS2 and ATM, there is a major factor needed to be taken into account. The compensated elevation of TanDEM-X represents the scattered surface of Ku-Band radar altimeter instead of X-Band interferometric phase center. Besides, the temporal difference of about 5 months between the generated DEM and ATM must be taken into consideration and corresponding seasonal variation of surface characteristics may be present. We strongly believe that the biggest contributors of the uncertainty are the 5month gap between the acquisitions of TanDEM-X and ATM as well as radar penetration. However, it is difficult to distinguish these two effects from a simple height comparison. Although the similar radar penetration was reported by Groh et al. (2014), the comparison of TanDEM-X DEM and ATM data as reported here is insufficient to yield any meaningful penetration as well as temporal variation of the surface because of significant dynamical surface changes during the 5month time difference. Groh et al. (2014) report surface lowering by several ３２ meters within one year (2011-2012) on the lower section of Thwaites glacier (Groh et al., 2014). Helm et al. (2014) also report thinning of Thwaites Glacier up to 10 m/year based on analysis of CryoSat-2 data (Helm et al., 2014). If the thinning rate (10 m/year) is assumed to be linear and spatially homogeneous, the surface decreases about 5 m during 5 months. This thinning implies that the penetration depth in this area is greater than 9 m. On the other hand, if the penetration is completely neglected, the result infers the surface increase of around 4 m within 5 months in an area where a surface lowering has been reported as previously mentioned. Although the result sounds contradictory, a monthly basis cross-point analysis between CryoSat-2 and ATM over the grounded ice sheet from 2011 to 2015 suggests that thickening in November 2011 is a possible event (Figure 2.7). While indicating surface thinning rate of -4.01 m/year, the relative mean height difference between June (blue) and November (red) in 2011 corresponding to the acquisition dates of TanDEM-X and ATM respectively is -4.94 m. What actually caused this thickening at this specific time of year needs additional analyses; however, the comparison between TDM and ATM may be explained without radar penetration. Besides, if seasonal surface variation, known as breathing is taken into consideration, then the penetration depth can vary again. Therefore, in our study, we decided not to interpret the calculated height bias as penetration depth because the actual interferometric penetration and the height variation during 5 month are indistinguishable. ３３ Figure 2.7 A cross-point analysis between CryoSat-2 and ATM on a monthly basis over the grounded ice sheet region between 75.35° S and 75.75° S shows a decreasing trend of surface height at a rate of -2.38 m/year. The red and blue square indicates June and November in 2011 corresponding to TanDEM-X and ATM acquisition, respectively. The relative difference between two periods is -4.94 m. ３４ Table 2.1. Specification of datasets Datasets TanDEM-X TerraSAR-X CryoSat-2 IceBridge ATM Acquisition Date Orbit Frequency/Wavelength 10 June 2011 Descending X-band (9.65 GHz) 09 June 2011 Descending 11 June 2011 Descending 13 June 2011 Descending 26 June 2011 Ascending 28 June 2011 Ascending 04 November 2011 Not applicable ３５ Ku-band (13.575 GHz) 532 nm Baselin e 106.06 m 1.17 m Table 2.2. Bias and standard deviation of height differences between TDM_CS2 and ATM (TDM_CS2 minus ATM) for each line. The lines are listed in the order to their locations from the north. Regional biases and standard deviations are also calculated. Line Number 1 9 2 8 3 7 4 6 5 Bias (m) -3.75 -1.05 1.00 2.19 4.05 -3.51 -5.62 -2.49 -6.01 Std (m) 18.97 21.49 7.70 4.31 4.35 7.62 6.95 6.42 7.60 Region Ice shelf Grounding zone Ice sheet Bias (m) 0.23 3.19 -4.31 Std (m) 15.50 4.42 7.27 ３６ High Resolution Grounding Line Mapping Using Hydrostatic Ice Thickness In this study, grounding line position is defined as the limit of hydrostatic potential where floating ice is in hydrostatic equilibrium state and grounded ice is not. Usually, the hydrostatic potential is calculated by comparing the measured and estimated ice thickness. Direct measurement of ice can be done by airborne ground-based radar surveys and it’s extremely difficult to obtain such datasets around the Antarctica. In contrast, thickness of the ice, denoted by hydrostatic ice thickness throughout this document, can be estimated when the relative distance between the surfaces of ice and sea level and the densities of ice and sea water are known when the hydrostatic equilibrium is satisfied. Therefore, grounding line mapping using a hydrostatic potential approach cannot be utilized when measured ice thickness is not available. In order to overcome this limitation of lack of field measurement, the bathymetry is employed instead because of its wide availability around the Antarctica. Instead of comparing the measured and estimated ice thicknesses, the bathymetry is compared with the bottom topography of ice, which can be estimated by subtracting hydrostatic ice thickness from surface elevation. Ice then can be considered floating when the estimated bottom topography is above the bathymetry, and grounded when below. This is another robust technique of monitoring grounding line position because it physically defines where the ice ３７ starts to depart from the bed. However, results can be strongly dependent on the accuracy of input datasets such as surface elevation, the bathymetry, densities of ice and sea water, etc. In this chapter, the position of grounding line is estimated by utilizing high resolution digital elevation models of glacier surface derived by combined usage of TanDEM-X and CryoSat-2 radar data as derived in Chapter 2. Then, accuracy of the results are evaluated in terms of ice density and the bathymetry. 3.1 Calculation of hydrostatic ice thickness of floating ice shelves In the absence of direct measurements, thickness of floating ice is usually calculated from surface elevation assuming a simple two layer ice column model (Figure 3.1). Flotation is determined where the mass of the displaced seawater is equal to the mass of an ice column, which contains ice and firn layer. The relationship can be written as follows: (ℎ𝑖𝑐𝑒 + ℎ𝑓 − F) 𝜌𝑠𝑒𝑎 = ℎ𝑓 𝜌𝑓𝑖𝑟𝑛 + ℎ𝑖𝑐𝑒 𝜌𝑖𝑐𝑒 (3.1) where F is ice freeboard, the elevation of the floating surface above sea level, hice, hf and ice, firn are the thicknesses and average densities of the ice and firn layers, respectively, and sea is the density of sea water. ３８ Figure 3.1 Schematic drawing of an ice/firn two layer model for hydrostatic thickness estimation. A firn layer is where the transformation of snow to ice occurs as the density of snow increases with depth. As the densification rate of a firn layer mainly depends on temperature, wind, and accumulation rate, the depth of a firn layer, defined as the depth where firn becomes ice, is spatially varying. The right hand side of the Equation (3.1) defines the total mass of an ice column with the firn layer on top of the ice layer. The total ice thickness (T) is defined as the sum of ice and firn layers (hice+ hf). It also can be defined as reduced total thickness of ice and firn column obtained by compressing the firn ３９ layer until it has the density of ice, which is about 900 kg/m3. When assumed the latter, the total ice thickness can be expressed as follows (Van den Breke et al., 2008): 𝑇= 𝜌𝑠𝑒𝑎 (𝐹 − 𝑑ℎ) 𝜌𝑠𝑒𝑎 − 𝜌𝑖𝑐𝑒 𝑑ℎ = ℎ𝑓 (1 − 𝜌𝑓𝑖𝑟𝑛 ) 𝜌𝑖𝑐𝑒 (3.2) (3.3) where dh is the firn depth correction, defined as the difference between the ice/firn column (hice+ hf) and the equivalent (T) ice thickness. The major component of ice thickness estimation comes from an incorrect average density of ice/firn column because an average ice column density is less well constrained around the Antarctica. Various studies dealing with Antarctic ice shelves suggest a column ice density ranging from 880 to 900 kg/m3 (Fricker et al., 2001). The depth of the firn layer is another huge contributor to the uncertainty which ranges from a few tens of meters to hundreds depending on the environments. Typically, the firn depth correction in the coastal region of Antarctica is known to be ranging between 12 and 20 m according to field measurement (Van den Breke et al., 2008). In this study, however, the firn depth correction is neglected, and instead, the density of ice column is varied from 700 to 918 kg/m3 because changing the ice column density can provide similar effects as using the firn depth correction. In addition, the uncertainties of thickness estimation may arise from freeboard estimation. Assuming that ice surface elevation is accurately ４０ determined by TanDEM-X interferometry in Chapter 2, sea level becomes the sole major contributor to the freeboard uncertainty. In this study, the regeneration of sea level at the time of interest is determined by combining the latest datasets of geoid and ocean tide model. For geoid, the gravity field combination model EIGEN-GL04C is utilized (Forste et al., 2008). This gravitational model is an optimal combination of GRACE and LAGEOS mission with high resolution gravitational information obtained from surface gravimetry and satellite altimetry data. It has been used in Antarctic researches with satellite altimetry as a reference geoid since its release up to date (Griggs and Bamber, 2011; Drews, 2015). In this study, EIGEN-GL04C geoid reference is obtained using Antarctic mapping tool for MATLAB (Greene et al., 2017). Ocean tide models are typically developed and distributed as gridded maps of the largest tide amplitudes or main waves. Among many ocean tide models, the FES2014 (Finite Element Solution) global tide model of which is the last version of the tide model developed in 2014-2016 and distributed from AVISO is utilized (Carrère et al., 2015). FES2014 was produced by NOVELTIS, LEGOS, and CLS Space Oceanography Division and distributed by AVISO, with support from CNES (http://www.aviso.altimetry.fr/). This model is based on hydrodynamic simulation with data assimilation of long-term satellite altimetry data and tidal gauges. One distinguishing feature of this model is that there had been special efforts to determine accurate tidal currents and to address the major non-linear tide (Lyard et al., 2006; Carrère et al., 2013). The accuracy of produced tidal signal has a strong spatial variability partly depending on ４１ distance from satellite altimetry data used for assimilation, bathymetry accuracy, and coastline complexity; the model has problems in producing the correct tidal signal in regions covered by the large ice shelves (Stammer et al., 2014). In addition, the inverse barometer effect (IBE) is one of the primary contributors to the accuracy of the reconstruction of sea level. The IBE is an elevation change associated with atmospheric pressure variability at a rate of -1 cm per hPa and has a typical height error over Antarctic ice shelves of ~0.08 m (Padman and Fricker, 2005). 3.2 Grounding line estimation using TanDEM-X and BEDMAP2: Campbell Glacier, East Antarctica Campbell Glacier (74° 25′ S 164° 22′ E) is one of the major glaciers flowing into Terra Nova Bay where Korea’s permanent Antarctic research base, Jang Bogo Station, run by the Korean Polar Research Institute (KOPRI), is closely located. Its width of terminus is approximately 110 km and the extent of the glacier forms Campbell Glacier Tongue retaining mean extent of 75.5 km2 over the years (Figure 3.2). There are two previously estimated grounding lines of Campbell Glacier in 1996 by Double Differential InSAR (DDInSAR) (Rignot et al., 2011) and in 2008-2009 by MOA (Bohlander and Scambos, 2007; Scambos et al., 2007). The disagreement between two grounding lines may have emerged from that their definitions of grounding line are different. DDInSAR locates the landward edge of a flexure zone, also known as hinge zone, where ice shelf is no longer influenced by the tidal displacement. In ４２ contrast, MOA defines a break in surface slope as a proxy for the grounding line. ４３ Figure 3.2 Campbell Glacier in East Antarctica. Grounding lines, estimated in 1996 by DDInSAR (Rignot et al., 2011) and in 2008-2009 by MOA (Bohlander and Scambos, 2007; Scambos et al., 2007), are overlaid in the flow velocity of the glacier and its glacier tongue. An inset indicates the location of Campbell Glacier as a red dot. ４４ Flow velocity of the Campbell Glacier Tongue was derived to be ~0.25 km/year by DInSAR (Han and Lee, 2015) and observation of surface elevation using ICESat GLAS showed that surface elevation has been decreasing at 0.69 m per year (Han et al., 2013; Han and Lee, 2014). This flow velocity and thinning rate is relatively insignificant in comparison with other glaciers in West Antarctica where the fastest velocity is nearly 4 km per year. In addition, the grounding line location hasn’t been moved much for nearly 20 years (Han and Lee, 2014) and the surrounding glaciers in East Ross Sea have not shown any significant changes (Konrad et al., 2018). Therefore, Campbell Glacier is chosen as an optimal site to evaluate robustness of the technique of generating TanDEM-X surface topography maps by comparing two independently generated elevation maps with enough time span in-between. Furthermore, these elevation maps can be used to estimate the grounding line of Campbell Glacier and monitor its change. 3.2.1 Digital Elevation Model of Campbell Glacier In overall, a technique of generation of a high resolution DEM using TanDEM-X was employed from Chapter 2. Two overlapping TanDEM-X bistatic data were acquired on June 21, 2013 and September 10, 2016. The perpendicular baseline for each pair is 214.59 m and 421.60 m and its corresponding ambiguity height is 30.24 m and 15.41 m, respectively. CryoSat2 data were acquired within 15 days from the acquisition of each TanDEM-X data in 2013 and 2016 to ensure enough number of candidates for ground control points. ４５ The power amplitude image shows that the scene contains both glacier and sea-ice covered ocean surface as the acquisitions were made in the austral winter (Figure 3.3a). High returned signals from the sea ice consequently caused high coherence in an interferogram (Figure 3.3b) leading an unexpected problem during the unwrapping of an interferogram that un-corrected unwrapped elevation over sea ice was higher than over ice shelf (Figure 3.3c). Thus, a linear fitting between CryoSat-2 elevation points and un-corrected unwrapped elevation caused a huge bias. In order to solve the unwrapping problem between the ice tongue and the sea ice, a decision was made to mask out the sea ice. The masking process was carried out by using the intensity only because the intensity over the sea ice is much weaker than over the glacier and ice tongue. After masking out the sea-ice covered ocean surface, a linear least-square fitting was carried out again between CryoSat-2 and un-corrected elevation data. In order to secure enough number of ground control points over the ice surface, a lower correlation threshold value and higher standard deviation of a CryoSat2 segment were used. In Chapter 2, the correlation threshold was 0.95 and the standard deviation was 10 m in Thwaites Glacier and there were enough CryoSat-2 elevation points for the absolute height calibration. In Campbell Glacier, the correlation threshold of 0.90 and the standard deviation of 20 m were applied to CryoSat-2 data. ４６ Figure 3.3 Amplitude (a), interferometric coherence (b), and unwrapped phase (c) of TanDEM-X SAR data in Campbell Glacier. The amplitude and coherence image show that the ocean is covered with sea ice. ４７ 3.2.2 Hydrostatic thickness and grounding lines in 2013 and 2016 Calculating hydrostatic ice thickness as in Chapter 3.1 and comparing the bottom topography of the ice with BEDMAP2 bathymetry, grounding line of Campbell Glacier is estimated. Tide levels were 31.103 cm on 2013/06/21 and 36.757 cm on 2016/09/10 and the density of sea water is set 1028 kg/m3. Due to the absence of field measurement, the density of the average ice column density is not known in Campbell Glacier. In general, ice density is less well constrained around Antarctica. Various studies dealing with Antarctic ice shelves suggest an average ice column density ranging from 880 to kg/m3 (Fricker et al., 2001; Wen et al., 2010). Jonathan (1994) used 917~918 kg/m3 in their study where they found a good fit in the comparison of satellite altimetry and ice thickness measurements. Thus, in this study, the firn depth correction is neglected, and instead, the average density is varied from 700 to 918 kg/m3 because changing the ice column density can provide the same effect as using the firn depth correction. As a starting point, the grounding line of Campbell Glacier is estimated using a well-used ice column density, 918 kg/m3 (Jonathan, 1994). After the ice bottom elevation is calculated using the derived ice thickness, it was compared with the bathymetry data from BEDMAP2. Generally speaking, the hydrostatic ice thickness is underestimated if the used ice column density is lower than the actual density. This case can be represented as a grounding line retreat scenario. The generated high resolution digital elevation model and estimated grounding lines of Campbell Glacier in 2013 and 2016 are shown in Figure 3.4. ４８ Figure 3.4 Generated DEM of Campbell Glacier and its ice tongue. Floating regions are colored differently to show the grounding line variations during 2013 and 2016. Green region indicates advance whereas red regions indicates retreat. Black solid lines indicate the grounding line in 1996. ４９ A solid black line represents the grounding in 1996, derived by DInSAR with ERS-1/2 datasets (Rignot et al., 2011). In order to distinguish advance and retreat of the derived grounding lines during the study period, floating ice regions are expressed in different colors. A grey area indicates an unchanged area where the ice was floating in both 2013 and 2016. Red and green area indicates retreat and advance of the grounding line between 2013 and 2016, respectively. Most of the retreat was observed in the western region of the ice tongue whereas the eastern region shows advance. However, the changes during the 3 years of observation were insignificant. The location of grounding line using hydrostatic approach varies significantly with different ice density values because ice density defines hydrostatic ice thickness. Thus, the grounding lines are again estimated using the average density ranging from 700 to 918 kg/m3 (Figure 3.5). Previously estimated grounding lines in 1996 and 2009 are also overlaid in red and yellow lines, respectively, for comparison. When the density is small (700~800 kg/m3), the ice thickness is absolutely underestimated resembling a huge retreat of the grounding line (black and darker regions in Figure 3.5). When the density is between 800 and 900 kg/m3, the grounding lines are located near the previously estimated lines inside grey regions. A line segment is drawn from the point A to A’ to see how the location of where ice bottom meets the seafloor as densities change (Figure 3.6). The grounding line in 1996 is somewhat matched only in the western region when the density is 900 kg/m3. However, the grounding lines aren’t matched at all at any density values in the eastern region. Assuming that there hadn’t been any ５０ significant thickness change since 1996, the ice column density of Campbell Glacier is approximately 900 kg/m3; however, it is not easy to conclude which ice density represents Campbell Glacier without any external measurement for validation. Figure 3.5 Grounding lines with respect to various ice densities. Previously estimated grounding lines in 1996 and 2009 are overlaid in red and yellow lines, respectively. ５１ Figure 3.6 A line segment from point A to A’ show how ice bottom varies and crosses the seafloor at different density values. 3.3 Grounding line estimation using improved bathymetry and additional field measurements: Thwaites Glacier, West Antarctica In this section, an improved high resolution bathymetry model from Operation IceBridge is utilized to estimate and update grounding line of Thwaites Glacier in West Antarctica. The same hydrostatic approach is employed using the surface elevation model, but multi-temporal datasets from 2011 to 2013 are put to use. In addition to an improved bathymetry model, the ice density and firn depth correction term are constrained by comparing with measured ice thickness from radar sounder data. To begin with, the surface topography of the marginal region of Thwaites Glacier is generated using TanDEM-X interferometry in combination with data ５２ from the near-coincident CryoSat-2 radar altimeter. A total of three TanDEMX DEMs are generated, which were acquired on June 10th 2011, June 29th 2012, and June 24th 2013 (Figure 3.7). The generated DEMs have a spatial resolution of 5 m. This is rather exceptionally high spatial resolution, which will be filtered and re-sampled to 25 m for matching with the field measurement data later. Sea level for each acquisition date is re-generated by combining geoid and tide model in order to retrieve freeboard of the ice shelf. According to FES2014 tide mode, geocentric tide levels were -22.50 cm in 2011, -21.72 cm in 2012, and -39.85 cm in 2013 (Figure 3.8). Subtracting sea level from the generated surface elevation of the ice shelf, the extracted freeboard is converted to the hydrostatic ice thickness using Equation 3.2. The used density of ice column and the firn depth correction are 913.08 kg/m3 and 45.35 m, respectively, which are constrained by matching the ice thickness from MCoRDS radar sounder data (Figure 3.9). In other words, the average ice column density and firn depth correction are retrieved so that the hydrostatic ice thickness is equivalent to the ice thickness measured by NASA’s Multichannel Coherent Radar Depth Sounder (MCoRDS) (Leuschen et al., 2010, updated 2018). The elevation, surface, bottom, and thickness for Antarctica were collected as part of Operation IceBridge funded aircraft survey campaigns. The data are provided with an along-track resolution of about 25 m and a sample spacing of about 14 m. In order to match with the generated elevation model, it was re-sampled to have a sample spacing of 25 m. Ice thickness is determined by converting the radar propagation time between the ice surface and bottom reflections which has a resolution of 17.8 m. The ５３ airborne survey was carried out on November 4th, 2011, of which there is a 5 months gap from the acquisition of TanDEM-X in 2011. ５４ Figure 3.7 Extracted DEMs in 2011, 2012 and 2013. CryoSat-2 SARIn data points (red) within 15 days from the TanDEM-X acquisition dates are initially chosen and rigorously selected as GCPs using correlation and RMS between the datasets. The used GCPs are shown in blue dots. ５５ Figure 3.8 Retrieved tide level from FES2014 tide model in June 2011 (a), June 2012 (b), and June 2013 (c). Tide level at the time of TanDEM-X acquisition is represented as a red diamond. ５６ Figure 3.9 a) A flight track of MCoRDS radar sounder of Operation IceBridge, acquired on November 4, 2011. b) Profiles of MCoRDS (black) and the estimated ice column thickness (red) along the flight line. The estimated thickness maps are then subtracted from the previously generated DEMs to construct bottom ice shelf topography maps. For comparison with the bottom topography of the ice shelf, the IceBridge Sander ５７ AIRGrav L4 Bathymetry data set of the grounding line region of Thwaites Glacier with a spatial resolution of 1 km, constructed from a series of parallel airborne survey lines, 10 km apart, is used. Errors are approximated as 70 m for the Thwaites dataset, incorporating errors in gravity measurement, radar ice thickness, ATM surface elevation, 2D model pinning point and some allowance for geological structures (Tinto and Bell, 2011). The bathymetry was interpolated to 25 m spatial resolution using spline interpolation method that is useful for interpolating unevenly distributed dataset (Sandwell and Ruiz, 1992). We expected that spline interpolation could render irregular topographic pattern of the gravity-derived bathymetry which is influenced by non-linear sediment densities. Figure 3.10-12 illustrate the grounding lines from this study, each represented in different colors. The most conspicuous changes are observable in the western and eastern edges of the catchment area in the middle of the scene. The area is where a deep trough is located and the grounding line retreat had been predicted from Tinto (Tinto and Bell, 2011). In the eastern edge of the grounding line, the grounding line had retreated ~ 1km/year at maximum whereas a smaller rate is observed in the western side. The grounding line in the middle area almost did not change at all during 2 years. Another significant observable feature of the estimated grounding lines is the evolution of an ice rumple, located 5 km offshore from the grounding line (Figure 3.12). This ice rumple was first sited in the interferogram (Rignot et al., 2014) but was never mentioned. Our result shows that the size of this ice rumple has gradually dissipated from 2011 to 2013. ５８ The grounding lines in this study show spatial migration, especially taking in place in the eastern part of the ice tongue. This gives an insight of the grounding line retreat that had been taking place as described in Rignot and Jacobs (2002). There can be numerous speculations about this, and yet the disappearance of large ice tongue and continuous rapid thinning of the glacier seemed to be the major contributor (MacGregor et al., 2012; Paolo et al., 2015). However, it must be reminded that this migration can also be observed mainly due to increase of spatial resolution. It is expected that more sequential TanDEM-X data in this region would be able to validate this result. The existence of a deep trough and its shape around the estimated grounding line may suggest the retreat might have been caused by the basal melting from the intrusion of circumpolar deep water as suggested from many publications (i.e., (Jacobs et al., 1992; Jenkins et al., 2010; Ha et al., 2014)). However, the information is not enough to conclude whether the retreat had been progressive or abrupt due to a lack of continuous observation in-between and it is expected that a continuous grounding line observation using the technique in this study would be able to fill the gap in the future to understand the stability of Thwaites Glacier. ５９ Figure 3.10 The estimated grounding lines in 2011, 2012 and 2013. ６０ Figure 3.11 Eastern part of the grounding line of Thwaites Glacier (Yellow: 2011/06/10, black: 2011/10/21, blue: 2012/06/29 and red: 2013/06/27) ６１ Figure 3.12 Western part of the grounding line of Thwaites Glacier (Yellow: 2011/06/10, black: 2011/10/21, blue: 2012/06/29 and red: 2013/06/27) 3.4 Sensitivity analysis with surface and bathymetry elevation The location of the estimated grounding line strongly depends on the accuracies of the hydrostatic ice thickness and bathymetry. The accuracy of hydrostatic ice thickness is determined by the accuracies of ice column density and surface elevation. The example in Campbell Glacier shows how the location of its grounding line can vary with respect of the used ice column densities (Figure 3.5) and that 880~900 kg/m3 is reasonable estimate of ice column density. Now, assuming a fixed ice density of 880 kg/m3 as reference, the grounding ６２ lines of Campbell Glacier are again estimated for various surface elevations (Figure 3.13). Although surface elevation of ice sheet/shelf generated using TanDEM-X and CryoSat-2 is believed to be accurate, it is helpful to know how much the location of grounding line can change by varying elevation. The grounding lines in darker colors in Figure 3.13 represent cases when the surface elevation is lowered. The maximum deviation of around 5 km is observable for 20 m elevation difference near the center of the ice flow. It is suspected that the bathymetry in this region is relatively flat or contains a trough. Likewise, spatial variation of estimated grounding line is analyzed with respect to varying bathymetry elevation (Figure 3.14). The magnitude of the variation of the bathymetry is set 10 times greater than of the surface elevation because the accuracy of the input bathymetry, BEDMAP2 in this case, is low. The result shows similar spatial variation of the grounding line as in the example with the surface elevation. ６３ Figure 3.13 Grounding line variation with respect to varying input surface elevations. The reference grounding line is estimated using a fixed ice column density of 880 kg/m3. The darker lines represent grounding line retreat as decreasing surface elevation. ６４ Figure 3.14 Grounding line variation with respect to varying input bathymetry elevations. The reference grounding line is estimated using a fixed ice column density of 880 kg/m3. The darker lines represent grounding line advance as increasing bathymetry elevation. The flow vectors are also overlaid in green arrows. In contrast to the increasing resolution of the surface topography, the spatial and vertical uncertainty of the available bathymetry still remains much higher than that of other datasets. As a result, the spatial uncertainty of the estimated grounding line caused by the vertical uncertainty of the bathymetry is much larger. Operation IceBridge campaign provides an improved bathymetry collected and generated from airborne survey measurements in Thwaites Glacier. It enabled the comparison of how much new bathymetry data can improve grounding line location. ６５ Figure 3.15 shows the estimated grounding lines using hydrostatic ice thickness, previously known grounding line using a DDInSAR technique, and the hydrostatic anomaly using an airborne radar sounder survey. Hydrostatic anomaly is calculated by subtracting a theoretical freeboard height in hydrostatic equilibrium from a measured orthometric ice surface elevation (Fricker et al., 2002). The anomaly is less than zero where ice is floating, and nonzero where ice is grounded or where there are errors in the measurements. This can be used as evidence for the grounding line position. ６６ Figure 3.15 Grounding line using AirGrav (white), BEDMAP2 (purple), and DDInSAR (black) are plotted over AirGrav bathymetry map. Hydrostatic anomaly using an airborne radar sounder data indicates where ice is floating (red) and grounded (blue). ６７ The most conspicuous differences between the grounding lines can be observed in black dashed rectangles, labeled as A, B, and C. Rectangle A shows that the grounding line estimated using hydrostatic ice thickness and improved AirGrav bathymetry successfully predicted where ice was floating whereas the grounding line estimated using hydrostatic ice thickness and BEDMAP2 and the line estimated using DDInSAR could not. Rectangle B shows a similar pattern where the grounding line using AirGrav was significantly different from the other two grounding lines. However, it is unclear which grounding line was true because of a lack of radar sounder survey crossing the region. Rectangle C shows that the grounding lines estimated using hydrostatic ice thickness successfully predicted where ice was floating whereas the grounding line using DDInSAR did not. In overall, the results clearly show that the estimation method using hydrostatic ice thickness can provide precise grounding line location. Furthermore, the estimation can be improved significantly as the quality and accuracy of input bathymetry increases. The question can arise why the grounding line using DDInSAR predicted differently in certain regions. It turns out to be that the certain areas are where a small trough is expected. Greenbaum et al. (2015) suggested that flotation over a narrow trough might not be detectible by the tidal flexure because only a small portion of full tidal deflection is reflected over such a region. Thus, the hydrostatic approach might be the only appropriate estimation method where bathymetry contains complex shapes and narrow troughs. Lastly, the assumption hydrostatic equilibrium might not be satisfied near the ６８ grounding line because of the supporting force from the grounded ice. Moving upstream from the fully floating region of ice shelf in the downstream, the grounding zone region contains an inflection point where the surface elevation is relatively lower than the surrounding area. This leads to the underestimation of the hydrostatic ice thickness. Then, the surface elevation usually increases with a large slope, depending on the basal topography, as moving upstream from the inflection point. This large gradient change is where the reduction in basal stress occurs as the ice starts float (Fricker et al., 2002). ６９ Local Scale Grounding Line Changes: A Case Study of an Ice Rumple in Thwaites Ice Shelf “Progressive degradation of an ice rumple in the Thwaites Ice Shelf, Antarctica, as observed from high-resolution digital elevation models” Ice rumples are locally grounded features of flowing ice shelves, elevated tens of meters above the surrounding surface. These features may significantly impact the dynamics of ice shelf grounding lines, which are strongly related to shelf stability. We used TanDEM-X data to construct high-resolution DEMs of the Thwaites Ice Shelf in West Antarctica and generated surface deformation maps allowing us to detect and monitor the elevation change of an ice rumple that appeared sometime between 1996 and 2011 when the grounding line of Thwaites Glacier retreated. The observed degradation of the ice rumple may relate to a loss of contact with the underlying bathymetry caused by the thinning of the ice shelf. We subsequently used a viscoelastic deformation model with a finite spherical pressure source to reproduce the surface expression of the ice rumple. Global optimization allowed us to fit the model to the observed deformation map, producing reasonable estimates of the ice thickness at the center of the pressure source. We conclude that combining the use of multiple high-resolution DEMs and the simple viscoelastic deformation model is feasible for observing and understanding the transient nature of small ice ７０ rumple, with implications for monitoring ice shelf stability. 4.1 Description of an ice rumple and study area Ice rumples are locally grounded and elevated features over which the surrounding floating ice shelves maintain their flow speed because the ice overrides the pinned bed such that horizontal divergence of the main flow does not occur (Macayeal et al., 1987; Matsuoka et al., 2015). These features are distinct from ice rises, which are locally accumulated and have their own flow dynamics including a radial ice-flow center separate from the main ice shelf (Favier and Pattyn, 2015). Ice rumples are typically elevated only tens of meters from the surrounding ice shelf surface (Matsuoka et al., 2015) and can be variably sized: the Doake Ice Rumples in the Ronne Ice Shelf are relatively extensive examples with an area of 1700 km2 (Vaughan, 1995), while smaller ice rumples can be only a few kilometers wide in other ice shelves (Fricker et al., 2009). Despite their relatively subtle physical features, ice rumples may have a significant impact on the dynamics of ice shelf grounding lines, which are closely related to ice shelf stability. Favier et al. (2012) demonstrated that an ice rumple can reduce shelf velocity and trigger the advance of the grounding line over the long term. Conversely, Tinto and Bell (2011) observed that the flow velocity of the eastern Thwaites Ice Shelf dramatically increased after the unpinning of an ice rumple. However, despite growing interest in and realization of the importance of local pinning points on ice shelf stability (Schmeltz et al., 2001; Furst et al., 2015; Berger et al., 2016; De Rydt et al., ７１ 2018), field observations are still difficult to obtain due to the presence of crevasses and rifts (Matsuoka et al., 2015). Additionally, surface responses due to short-term progressive basal perturbation have not previously been discussed in the literature. Such analysis is thus required, especially in West Antarctica where many ice shelves are experiencing rapid and notable thinning. Moreover, continuous monitoring using satellite techniques is also challenging because the detection of ice rumples usually requires a minimum of four separate SAR acquisitions for double-differential interferometry (DDInSAR). Recently, the availability of high-resolution satellite imagery has enhanced our ability to detect and observe ice rumples. One of the main objectives of TanDEM-X, a twin satellite of TerraSAR-X that orbits in a helix formation for various interferometry modes, is to construct DEMs of the globe with a vertical accuracy of 2 m and spatial resolution of 25 by 25 m. This enables the observation of ice rumples’ topographic features and the monitoring of changes when data are available in time series. In this study, I used multiple TanDEM-X datasets to study an ice rumple in Antarctica’s Thwaites Ice Shelf, previously observed via DDInSAR by Rignot et al. (2014). The pinning point appeared sometime between the initial observations in 1996 and the retreat of the Thwaites Glacier grounding line in 2011. Subsequently, I compared each DEM in an original way by anticipating surface deformation of the ice rumple because the Thwaites Ice Shelf has been experiencing rapid thinning recently (Paolo et al., 2015). Finally, I investigated the surface response due to short-term progressive basal perturbation using a simple deformation model, assuming the grounded bed underneath the ice acted ７２ as a spherical stress source. The ice rumple studied here is located on the west Antarctic coast, approximately centered at 75°23’ S and 106°42’ W with a radius of 1 km in both the longitudinal and along-flow direction in 2011. The elevation anomaly from the surrounding region was approximately 35 m, supporting its characterization as an ice rumple. The location was approximately 5 km downstream of the grounding line estimated in 2011 by Rignot et al. (2014) and had not moved since its first observation in the same study (Figure 4.1). This ice rumple was not observable in 1996 because the grounding line was located considerably in the downstream from the ice rumple (green line in Figure 1). This is a different pinning point from the one observed in the eastern ice shelf by Tinto and Bell (2011). ７３ Figure 4.1 Grounding lines of Thwaites Glacier in 1996 (green) and 2011 (red) estimated using the DDInSAR method with European Remote Sensing (ERS) satellites (Rignot et al., 2014). The orange dotted rectangle indicates a larger ice rumple previously discussed by Tinto and Bell (2011). A newly generated digital elevation model derived from TanDEM-X data is shown within the diagonal rectangle, in which the small red feature inside the yellow dotted square indicates the smaller ice rumple considered in this study. ７４ Thwaites Ice Shelf is the seaward extension of Thwaites Glacier which is one of the major glaciers in West Antarctica exhibiting significant retreat of grounding line (Rignot et al., 2014) and ice shelf margins (MacGregor et al., 2012), negative mass balance, and ice volume loss (Paolo et al., 2015) in recent decades. The intrusion of the Circumpolar Deep Water into the ocean base underneath the glaciers has been considered the predominant cause of such immense changes. The ice shelf flows at high velocity, exceeding 4 km per year at the fastest. It is also known for strong katabatic wind and rapid and heavy precipitation causing dramatic surface change at any time. In the vicinity of the grounding lines, the surface elevation ranges approximately from 50 to 300 m above the WGS84 ellipsoid. 4.2 Generation of high-resolution deformation maps I determined the local topographic variation of the floating ice shelf by differencing pairs of independently generated TanDEM-X DEMs of the Thwaites Ice Shelf. The procedure of DEM generation by calculating interferograms with two bi-static TanDEM-X data is explained in more detail in Kim and Kim (2017). TanDEM-X data were acquired on June 10, 2011, June 29, 2012, and June 27, 2013, allowing us to generate two sets of deformation maps; the time spans between each acquisition were 1.05 (385 days) and 0.99 (363 days) years. Securing reliable ground control points is key to generating DEMs using TanDEM-X data over snow-covered surfaces; the Thwaites Ice Shelf is a particularly challenging region due to its high flow velocity and abundant ７５ crevasses. I followed the methods described by Kim and Kim (2017), who showed that near-coincidental CryoSat-2 radar altimeter datasets can provide reliable reference elevations over the surface after a rigorous selection process using normalized local pixel cross-correlation between CryoSat-2 and uncorrected TanDEM-X elevation as well as the standard deviation of a segment of CryoSat-2 elevation points. I used CryoSat-2 tracks within 15 days of the TanDEM-X acquisition as possible ground control points, setting the threshold values as 0.95 and 10 m for cross-correlation and standard deviation, respectively. DEMs were generated with a 25 x 25 m grid size with a spatial resolution of ~3 m along track and ~1.2 m in radar line of sight (Rossi et al., 2012) and a vertical resolution of better than 2 m over smooth areas (Becek et al., 2016). The elevation was defined with respect to the WGS84 ellipsoid. In order to extract distinct surface elevations of the ice rumple, 2D Gaussian filtering was applied to diminish the effect of the crevasses, which can otherwise contaminate the DEM difference calculations. Elevation change was calculated by simply subtracting geocoded and layerstacked DEMs in UTM projection. Geocoding errors were considered negligible because the flight information of TanDEM-X is known to be precise and geocoding was done using the flight information. Besides, horizontal and vertical accuracy was considered satisfied because each DEM was spatially and vertically compared with CryoSat-2 SARIn points. The elevation difference contained the overall thickness change over the study period as well as the vertical offset due to mismatched tidal heights. Thus, these tidal differences had ７６ to be eliminated in order to accurately extract the surface deformation due solely to changes in the ice rumple. Tidal heights were estimated using the FES2014 tidal model, produced by NOVELTIS, LEGOS, and CLS Space Oceanography Division and distributed by AVISO, with support from CNES (http://www.aviso.altimetry.fr/). Tides heights estimated at 75°23’ S and 106°68’ W were -22.50, -21.72, and -39.38 cm at the time of each TanDEM-X acquisition in 2011, 2012, and 2013, respectively. The tide correction was done to the floating ice shelf using the grounding line in 2011. High resolution surface topography maps of the Thwaites Ice Shelf from 2011 to 2013 were generated using bi-static interferometry of TanDEM-X and TerraSAR-X SAR data with selected CryoSat-2 radar altimeter elevation data as ground control points in calibrating absolute elevation (Figure 4.2). For each topography map, there were 7 overlapping altimeter tracks with 222, 250, and 320 elevation points before the selection in 2011, 2012, and 2013, respectively (dark grey circles in Figure 4.2). The selection process then sorted out 22, 31, and 39 points with mean cross-correlation coefficients of 0.97, 9.98, and 0.98, respectively (yellow circles in Figure 4.2). Using these selected elevation data as ground control points, the correction for absolute elevation was carried out. The RMSEs between the ground control points and the corrected surface topography maps were 5.26, 2.02, and 1.96 m in 2011, 2012, and 2013, respectively. ７７ Figure 4.2 Generated high resolution surface topography maps of the Thwaites Ice Shelf on (a) 2011/06/10, (b) 2012/06/29, and (c) 2013/06/24. Dark grey denote the locations of all of the available CryoSat-2 SARIn altimeter elevation data and red circles denote selected data. ７８ The ice rumple had a circular shape with a downstream slope steeper than that upstream (Figure 4.3), resembling the simulation results of Gudmundsson et al. (1998) showing that surface undulations are in phase with basal perturbations under a high basal slip ratio. It appears that the peak of the ice rumple advanced downstream by 125 m and 340 m while decreasing in elevation by 11.67 m and 21.02 m during the 2011–2012 and 2012–2013 periods, respectively. Both deformation maps show surface depressions with a circular main shape and a weak elliptical signal in the cross-flow direction (Figure 4.3 d,e). The peak of the deformation advanced downstream by 250 m and the maximum deformation was -14.55 m and -24.23 m during the 2011–2012 and 2012–2013 periods, respectively. The fact that the surface expression was not advected downstream during the observation period shows that the ice shelf was still in contact with the basal topography. Although it did advance downstream by a few hundreds of meters each year, if the ice shelf was ungrounded the surface expression would have already disappeared from the defined map window since the overall flow velocity in this region is nearly 3 km/year. The differences in the elevation and deformation peaks are related to the asymmetric shape of the ice rumple and the nearly circular shape of the deformation pattern. There is also a dislocation between the center of deformation and the physical peak of the observed ice rumple (Figure 4.3 d,e). ７９ Figure 4.3. Detailed surface DEMs of the ice rumple from 2011–2013, its yearly surface deformation, and surface velocity (arrows): (a-c) Location of maximum elevation in each of the three years (filled circles indicate the peak in that year); (d-e) Deformation maps for the periods 2011–2012 and 2012–2013 (later year subtracted from earlier year) where cyan circles indicate maximum surface depression and black circles indicate maximum elevation in each reference period. 4.3 Interpretation of surface deformation using of a viscoelastic deformation model An ice rumple is formed when a floating ice shelf is locally grounded on the bed. As the bed in contact with the ice is acting as a local stress source, the elevation of the ice rumple is determined by the extent to which the ice shelf is grounded. Therefore, surface deformation of an ice rumple occurs when a local ８０ stress underneath the ice changes. Assuming that the ice is materially homogeneous, that its mechanical properties do not vary with direction, and that there is a linear relationship between applied stresses and strain at any point in the ice, such surface deformation can be modeled as a result of varying spherical pressure in an ideal semi-infinite elastic half-space. Although the actual shape or size of a pressure source (the extent of grounding) is unknown, the spherical source is a practical initial proposition (Figure 4.4). Figure 4.4. Schematic drawings of ice rumple (a) and the elastic deformation model with a finite spherical pressure source (b). The initial surface is not flat, thus the radial distance from the center of a sphere to the surface varies. The point dilatation model described by Mogi (1958) is the most well-known model in such cases and has been widely used in studies of volcanic magma chambers. Mctigue (1987) derived a more elaborate expression to analytically approximate deformation that includes higher-order terms accounting for the finite shape of a spherical source. The latter model improved the Mogi model by allowing the source to be shallow and introducing the size of the body. Thus, it determines not only the depth to the body but also its radius and the pressure change. Despite the half-space representing one bounding surface extending ８１ infinitely in all directions, which did not reflect the real floating ice shelf, I considered the simple elastic deformation model to be capable of providing good approximations of deformation resulting from pressure changes within the shallow ice. I assumed the linear viscoelasticity of ice by considering a combination of elastic and viscous behavior, based on the conclusion of Doake and Wolff (1985) that Antarctic ice shelves that may be influenced by the restraining force from ice rumples show a linear flow characteristic even though conventional flow laws for ice consider non-linear stress-dependent viscosity. In addition, Gudmundsson et al. (1998) concluded that rheological non-linearity can be ignored for rapidly sliding ice where the basal perturbation (with wavelengths of a few ice thicknesses) is easily transferred to the surface. I assumed that the ice shelf was rapidly sliding as the friction underneath it changed from positive to zero when moving seaward from the grounding line. In addition, I utilized the viscoelastic correspondence principle to re-derive the viscoelastic deformation of the ice rumple, expressed in the form of the corresponding elastic deformation with the shear modulus for a viscoelastic rheology. Given an exact solution to the elastic problem, any viscoelastic rheology can be implemented by taking the Laplace transform from the elastic solution and substituting the effective shear modulus containing both elastic and viscous properties. Then, the viscoelastic deformation solution can be found by taking the inverse transform. In this case, I implemented the effective shear modulus for a standard linear solid (Peltier, 1982). The elastic deformation at the free surface (z=0) in the vertical direction (uz) ８２ due to a pressure change (∆P) within a spherical source located at x0, y0, and z0 is given as follows: 𝑢𝑧 (𝑟, 𝜀, 𝜇, 𝜈) = (1 − 𝜈) − 𝛥𝑃𝛼 3 𝜇 𝑧0 (𝑟 2 + 3 𝑧02 )2 × [1 (4.1) 𝛼3 1+𝜈 15(2 − 𝜈) 𝑧02 [ − ]] 𝑧0 2(7 − 5𝜈) 4(7 − 5𝜈) 𝑟 2 + 𝑧02 where r is the radial distance from the free surface above the center of a spherical source, z0 is the depth of the center of a spherical source (a positive number), ν is Poisson’s ratio, α is the radius of the sphere, and µ is the shear modulus that can be expressed with Poisson’s ratio and Young’s modulus. I used values for the latter two (0.3 and 0.88 GPa, respectively) from Vaughan (1995), although high variability in effective Young’s modulus has been observed depending on the location and time scale of ice deformation (Schmeltz et al., 2002). The deformation here was normalized by ∆Pz0. The pressure change here is analogous to the change in uplifting stress at the bottom of the ice shelf due to the varying extent of grounding from either thinning or thickening. The Laplace transformed effective shear modulus for a standard linear solid in time-domain can be found in Peltier (1982) and is given by: μ(s) = 𝜇 𝜇(𝑠+ ⁄𝜂) (2𝜇)⁄ , 𝑠+ 𝜂 (4.2) Making the substitutions (2) for the Laplace transformed (1), the deformation ８３ problem becomes a viscoelastic deformation model in a standard linear solid rheology. Taking the inverse transform, the solution is given by: 𝑢𝑧 (𝑟, 𝜀, 𝜇, 𝜈, 𝜂, 𝑡) = 1 −𝜇𝑡 1 Δ𝑃𝑧0 (𝜈 − 1) 𝑐1 𝑓1 (𝑟) 𝑓2 (𝑟, 𝜀) × [1 − 𝑒 𝜂 ] 𝜇 2 2 𝑐1 = 𝜀3 (7 − 5𝜈) 𝑓1 (𝑟) = 1 5 (𝑟 2 + 1)2 𝑓2 (𝑟, 𝜀) = 𝜀 3 (17𝜈 − 28) + 2𝜀 3 𝑟 2 (1 + 𝜈) − (𝑟 2 + 1)(20𝜈 − 28) (4.3) (4.4) (4.5) (4.6) where  is viscosity and t is time. I used a viscosity of 50.1 TPa s (Wild et al., 2017). The free surface was represented by the topography of the ice rumple as determined by the generated DEMs, such that it was not flat in this case. I expected that having precise topography would increase the accuracy of the results, because surface deformation is known to be strongly affected by topography (the distance of the free surface from the spherical source) (Williams and Wadge, 1998). In order to determine the depth (z0), radius (α), and pressure change (∆P) of a spherical source beneath the ice rumple as well as its center (x0, y0) and time (t), I adopted a global optimization of the constrained nonlinear multi-variable problem using MATLAB software. The objective function in this case was the root-mean-square (RMS) difference between the observed deformation from DEM differencing and the estimated deformation calculated from the elastic ８４ deformation model. Constrained optimization involves finding values of a multi-variable vector that are a local minimum to an objective function subject to constraints on the variables. Each variable holds specific conditions such as inequality, equality, and boundary constraints. I adopted the interior-point algorithm with randomly selected initial values within the boundary conditions (Table 4.1). Given the observed surface deformation as the reference in an objective function, the global optimization with constraints was used to estimate the horizontal location and depth of a pressure source’s center along with its radius, pressure increment, and time. Optimization was carried out by finding the minimum of the objective function that is the RMS difference between the observed and estimated deformation. Six unknown variables and the minimum and maximum constraints of each are listed in Table 4.1. ８５ Table 4.1. Minimum, maximum, and optimized values of each unknown variable for constrained global optimization of the viscoelastic surface deformation. Variable 11-12 Maximum 11-12 x0 (m) -500 500 -22.89 -22.91 -233.67 -233.67 y0 (m) -500 500 134.72 134.71 -207.78 -207.78 z0 (m) 500 2000 1028.08 1027.06 878.48 878.45 0 500 259.39 228.80 133.03 131.37 0 -10 -1.06 -1.55 -8.53 -8.86 0 2 1.64 1.05 0.66 0.99 Radius (m) Pressure change (GPa) Time (yr) ８６ (w/o time) 12-13 12-13 Minimum (w/o time) The modified deformation model accurately reproduced the near-circular pattern observed in the surface deformation for both periods. The model predicted a pressure source at a shallow depth of 1028.08 m with a radius of 259.39 m and pressure increment of -1.06 GPa during 2011–2012 (Figure 4.5). Figure 4.5. Comparison between the observed and estimated deformation from the modified deformation model after the global optimization process during 2011–2012. The RMS difference between two deformation maps was 1.53 m. Profiles in the longitudinal and along-flow directions are plotted along the red dashed and dot-dashed lines, respectively. The RMS difference between the observed and estimated deformation was 1.53 m. There was a slight offset of the pressure source’s center location compared to the peak of observed deformation. The modified deformation model predicted that the viscoelastic deformation took place over 1.64 years (compared to the actual 1.05 years). Likewise, the model predicted a pressure source at a depth of 878.48 m with a radius of 133.03 m and pressure increment of -8.53 GPa during 2012–2013 ８７ (Figure 4.6). The RMS difference between the observed and estimated deformation was 5.55 m, significantly larger than the previous period. Again, there was a slight offset of the pressure source’s center location compared to the peak of observed deformation. The modified deformation model predicted that the viscoelastic deformation took place over 0.34 years, significantly different than the actual time (0.99 years). Given the large difference in the estimated deformation times, the optimization was re-processed with only five unknown variables (setting time as the actual elapsed time), with the results shown in Table 4.1. The modified deformation model predicted a similar result during 2011–2012. The depth of the pressure source was almost the same with a slightly reduced radius of 228.80 m and increased pressure increment of -1.55 GPa. The RMS difference was the same as in the previous estimation, 1.53 m. There was a significant improvement in the results for the 2012–2013 period where the RMS difference decreased from 5.55 m to 2.11 m (Figure 4.7). ８８ Figure 4.6. Comparison between the observed and estimated deformation after the global optimization process during 2012-2013. The RMS difference between two deformation maps was 5.55 m. Profiles in the longitudinal and along-flow directions are plotted along the red dashed and dot-dashed lines, respectively. ８９ Figure 4.7. Comparison between the observed and estimated deformation from the modified deformation model after the global optimization process during 2012–2013. After the estimated deformation time was set to the measured time (0.99 years), the RMS difference decreased from 5.55 m to 2.11 m. ９０ 4.4 Discussion The instant deformation map alone cannot uniquely resolve whether such deformation is caused by a progressive or ephemeral process. For example, surface deformation could also have been the result of tidal displacement instead of less grounding due to ice thinning. If the deformation were due to tidal displacement, the observable features of an ice rumple must have been present periodically throughout the study period. In order to validate that the deformation signal between the measurements was indeed due to progressive thinning of the floating ice, I monitored the surface features of the ice rumple using Landsat 7 ETM+ images from 2003–2014 (Figure 4.8). The images in this time series showed a gradual dissipation of the ice rumple, strongly indicating continuous thinning of the Thwaites Ice Shelf. Furthermore, that the disappearance of surface features (e.g., crevasses and surface gradient) from 2013 onwards suggests that the ice shelf has been ungrounded, removing the pressure point that had been maintaining the ice rumple. ９１ Figure 4.8. Landsat 7 ETM+ images from 2003–2014 showing the gradual disappearance of the studied ice rumple in the Thwaites Ice Shelf. Crevasses and surface gradients are generally created atop an ice rumple due to surface extension and elevation increase. Such features were visible as late as 2011 but disappeared by 2013, indicating gradual ice thinning. Larger images are magnifications of selected areas indicated by red boxes. ９２ As linear theory is invalid under high stresses, the presence of surface crevasses in the area might suggest that the ice shelf was influenced by high stress. Landsat images show that crevasses were generated at a local pinning point and high stress seemed to be affecting the ice in the early years (e.g., 2003 in Figure 4.8). However, the magnified images show that the surface crevasses present around the ice rumple since 2011 originated from the grounded region upstream of the rumple. Thus it is likely that the stress beneath the ice shelf decreased significantly due to ice thinning and was not sufficient to create crevasses at the pinning point. Therefore, I conclude that the ice shelf was experiencing low stresses underneath and that linear viscoelasticity was a reasonable approximation in this study as the TanDEM-X data were acquired in 2011, 2012, and 2013. Total deformation can be decomposed into elastic and viscous deformation because linear theory is employed. During 2011-2012 period, the modelled deformation approximately consists of 50% elastic and 50% viscous deformation (Figure 4.9). In this study, the ice viscosity and Young’s modulus were 50.1 TPa s and 0.88 GPa and the percentage of deformation can be determined with respect to viscosity and Young’s modulus for ice. A full range of deformation mixtures seems possible because ice has variable viscosity and Young’s modulus depending on stress and thickness, respectively. ９３ Figure 4.9 The viscoelastic deformation model is fitted to the observed deformation (a). The modelled deformation (c) is a summation of viscous (b) and elastic deformation. ９４ Although a low RMS difference is not necessarily an indication of the validity of the model, the model successfully fit the observed deformation with a plausible pressure. Nevertheless, this simple model alone cannot determine the actual shape of what is certainly a more complicated bathymetry, and I can only postulate the extent of grounding using the radius of a pressure source. In this context, I suggest that the initial ice thickness can be defined as the difference between depth and radius as shown in Figure 4.10. Then, thickness change can be derived as the variation of thicknesses at different time durations only when there are at least three independent observations of the surface. However, this interpretation must be approached carefully because the thickness change is valid only if a pressure source is located at the same center. ９５ Figure 4.10. Schematic drawing of surface deformation in consecutive periods and corresponding ice thicknesses, defined as the difference between depth and radius of a spherical pressure source. According to the fixed-time optimization result, the ice thickness was 798.27 m and 747.08 m in 2011 and 2012, respectively. This is quite different from the known ice thickness of the Thwaites Ice Shelf near the grounding line (~1 km), despite the low vertical resolution of the radar sounder used for such measurements (Tinto and Bell, 2011). This indicates a thickness decrease of 51.19 m during that one-year period. However, as the center of the pressure source in 2012 was located 342.49 m upstream and 210.76 m to the west from that in 2011, it is difficult to substantiate the thickness and thickness change with the obtained datasets. Our results for the surface depression and thinning of the ice rumple in the ９６ Thwaites Ice Shelf were much higher than previously reported. It is rather surprising to observe such high (> 10 m) surface depressions in an ice shelf in such a short time; Rignot et al. (2013) and Paolo et al. (2015) reported that the thinning rate of the Thwaites Ice Shelf was 6.13 m/year and 2.80 m/year, respectively. A major factor contributing to these differences may relate to the extent of the survey area, which was over 4,000 km2 (the entire ice shelf) for the other studies as opposed to a few tens of km2 focusing on a single ice rumple in our study. Elevation and deformation profiles over the ice rumple and surrounding ice shelf show two distinguishable thinning rates (Figure 4.11). The surface depression was much more significant over the ice rumple than over the ice shelf. Thus, it should not be assumed that our results represent the overall thinning of the entire ice shelf. ９７ Figure 4.11 Elevation and deformation profiles over the ice rumple and surrounding ice shelf. ９８ Despite this narrow focus on a local change, such progressive thinning may cause the un-grounding of a pinning point and the acceleration of the entire ice shelf (Tinto and Bell, 2011); it is important to realize when this has happened. Based on the Landsat imagery time series in Figure 4.8, un-grounding might have taken place sometime after 2013 due to the disappearance of the surface features. It is unclear at this point whether the high flow velocity of Thwaites Ice Shelf has been caused by the progressive unpinning of an ice rumple. However, I suspect that the ice rumple in this study did not have any significant impact on the stability of Thwaites Ice Shelf because of the timescale of the observed changes and the already high flow velocity of the ice shelf. The optimization of the multi-variable problem successfully estimated unknown parameters within the constraints of reasonable values (Table 4.1). It is intriguing that the simple viscoelastic deformation model successfully replicated the deformation of the ice rumple in the Thwaites Ice Shelf with a small RMS difference; the estimated values of the parameters were also plausible. The main cause of the remaining difference is the slight elliptical shape of the deformation, which can take a different form in other ice shelves. In such cases, deformation models for other forms of pressure sources (i.e., sills, dykes (Okada, 1985), and disks (Fialko et al., 2001)) could be used to fit the observed deformation. In order to analyze the dominant factor in the optimization result, I conducted simulations of the deformation model by varying the estimated parameters around a nominal value (Figure 4.12). The simulation outputs were induced by varying one parameter at a time while keeping all others fixed. The nominal ９９ reference values were chosen as the result of the optimization (Table 4.1), and the variations were made within 10% from the nominal values. The greatest variation in the result, i.e., the difference between the observed and estimated deformation, came from changing the radius of the spherical source (maximum variations of 1.88 and 2.51 m) from the original RMS difference; the smallest change was observed when varying the location of its center. １００ Figure 4.12. The contribution of variations of the input parameters to the output of the surface deformation model was investigated by changing the variables within 10% from the estimated values. The center location of a spherical source showed the least sensitivity, whereas the radius made the highest contribution. １０１ Conclusions In this study, high resolution digital elevation model and hydrostatic potential of floating ice are utilized to estimate, update, and monitor grounding lines of Antarctic ice shelves. The most important aspect of this study is the potential of monitoring grounding line migration in high spatial and temporal resolution by using a simple and repeatable estimation technique. In Chapter 2, a high resolution digital elevation model in the marginal region of marine glacier was successfully generated by introducing combined usage of TanDEM-X and TerraSAR-X bi-static interferometry and CryoSat-2 SARIn radar altimetry. Instead of estimating phase offset, the usage of near-coincident radar altimeter as ground control points determined absolute topographic height exploiting linear least-squares algorithm. Reliable GCPs were also carefully selected using normalized cross-correlation and standard deviation ensuring high accuracy of radar altimeter data. This technique was applied to Thwaites Glacier in West Antarctica and the generated elevation map was compared with airborne laser altimeter which was acquired 5 months after the acquisition of TanDEM-X SAR data. The result showed that a large height bias was observed on the ice tongue region at where intense activities were continuously occurring. Meantime, the result also showed different spatial trends; surface lowering was observable near the grounding line (positive bias of 3.19 m) whereas laser altimeter represented higher elevation on the grounded ice sheet suggesting radar penetration １０２ (negative bias of -4.31 m). However, the height difference between two epochs was considered not enough to conclude temporal trend of which more time series data will be required. Besides, it must be stated that the accuracy of the generated elevation model is expected to be different if the laser measurement with less time gap is available. Even without considering the height change due to the time gap, the current result contains errors arising from the radar penetration difference between TanDEM-X interferometry and CryoSat-2 altimetry as well as the difference from the laser measurement. If at least the error of the time gap can be reduced, more precise analysis will be possible. For generating a high resolution digital elevation model, the most crucial part of this application was the presence of simultaneous dataset that could be used as ground control points in calculating absolute height. Besides, the importance of the selection process of reliable ground control point must be emphasized as well. After all, this study has shown that TanDEM-X can provide highresolution elevation information even in the area without any ground truth data by using near coincident CryoSat-2 radar altimeter data. In Chapter 3, grounding line estimation was carried out by exploiting hydrostatic equilibrium of ice shelf using generated high resolution elevation map. Grounding line was defined where the bathymetry met the bottom of ice shelves of which was calculated by subtracting hydrostatic ice thickness from the surface. Hydrostatic ice thickness was calculated from freeboard height which required high resolution topography and sea level at the time of the acquisition. Sea level at the acquisition time of TanDEM-X and TerraSAR-X １０３ was determined by combining EIGEN-GL04C geoid and tide level derived by FES2014 global tide model. Finally, the bathymetry data were employed from BEDMAP2. This technique was applied to Campbell Glacier in East Antarctica and Thwaites Glacier in West Antarctica and the results suggest that continuous observation and monitoring of grounding line is feasible using multi-temporal surface topography datasets. TanDEM-X data were obtained so that there were at least two overlapping acquisitions with time interval for each ice shelf. Multiple data acquisitions provided temporal variations of the estimated grounding lines at each ice shelf. Such temporal variations were then exploited to better understand dynamics of each ice shelf. In Campbell Glacier, TanDEM-X and TerraSAR-X data in single-pass interferometry mode were acquired on June 21, 2013 and September 10, 2016. During 3 years between two data acquisitions, any significant changes in grounding line position was evident. There was a spatial difference from the grounding line estimated by DDInSAR whereas the estimated grounding line using the characteristics of the surface of the optical satellite images agreed well when the ice column density was about 880 kg/m3. Nevertheless the variations of the grounding line was also analyzed with respect to the density of ice because the average density of ice used in estimating grounding line was not accurately known. So far, generalized input datasets such as BEDMAP2, fixed ice density, and firn depth correction was applied because field measurements are rarely available in Antarctica. Although the accuracy of the grounding line estimation strongly depends on the accuracy of used datasets, newly estimated １０４ grounding lines, when compared with the previous grounding lines, showed its feasibility of detecting the spatial evolution of grounding line in time series. Thwaites Glacier is one of the fastest changing ice bodies in Antarctica; thus, there are extensive airborne and ground-based surveys available, such as Operation IceBridge. Thus, the average ice column density and firn depth correction were able to be constrained using field measurements from Operation IceBridge. There were three overlapping TanDEM-X datasets from 2011 to 2013 in Thwaites Glacier. Grounding lines were estimated using each dataset and they were compared and validated with the previously updated grounding line from DDInSAR. Comparison with the DDInSAR result showed subtle local differences where the most conspicuous difference was observed over bathymetric troughs. During 2011-2013 period, the grounding line retreat was also observed over the troughs with the maximum retreat rate of nearly 1 km/year. Since there exists new bathymetry data around Thwaites Glacier derived from airborne gravimetry survey by Operation IceBridge, grounding line was again estimated using new bathymetry data. The result confirmed that using datasets with better vertical and spatial resolution and accuracy could provide much improved and more accurate grounding line. Especially, there was a considerable change in regions where the bathymetry had a shape of bay with narrow width. Besides, an ice rumple with a radius of ~1 km was detected in this study. The pinning point was located nearly 5 km offshore from the previously estimated grounding line in 2011 and appeared sometime between 1996 and 2011 when the grounding line of Thwaites Glacier retreated. Many such pinning points １０５ exist around the margins of Antarctica and they are important features that control the flow of the ice shelves. In Chapter 4, evolution of this local feature is exploited in more details. From 2011–2013, the deformation maps showed the recent fading of this ice rumple. The deformation pattern, along with a time series of Landsat 7 ETM+ imagery, showed that the ice was still in contact with the basal topography as late as 2013 but has since been unpinned. Then, the deformation maps were used with the simple viscoelastic deformation model (widely used in volcanic studies) to interpret the surface changes in terms of pressure changes at the bottom of the ice shelf by applying an idealized spherical pressure source in Chapter 4. The deformation model was fitted with the observed deformation to yield a center location and depth of the pressure source underneath the ice. The estimated numbers were reasonable and the ice shelf thickness at the center of the spherical pressure source was also estimated using the depth and radius. This study was designed to observe and better understand the transient nature of an ice rumple’s surface elevation. Such detailed monitoring was made feasible for the first time by our use of TanDEM-X and near-coincidental CryoSat-2 data to derive DEMs with sufficient spatial and vertical resolution. Despite the extreme surface conditions of the Thwaites Ice Shelf, which is heavily crevassed and has high flow velocity, this method has shown the great potential for locating previously unknown small ice rumples or continuously monitored such features in other ice shelves. The surface depression and thinning of this ice rumple were found to be much higher than those of previously reported levels for the broader region, though this may be explained １０６ by the more local focus of this study (a few tens of km2) compared to other research on the area. Thus, the results may not represent an overall thinning process, but do suggest that significant strong melting can exist at a local pinning point. Above all, this study contributes to our understanding of the transient nature of rapidly changing ice shelves such as the Thwaites Ice Shelf. In addition, the results of the viscoelastic deformation model have demonstrated the potential to estimate the location of an uncharted pinning point as well as the ice thickness at the time of the observation. These findings are expected to be of interest to the future modeling of ice shelf instability as initial input data. As several questions remain to be answered regarding this process, more research using an increased number of datasets at different regions will be required. Finally, conceivable error sources are still present in grounding line estimation such as ice column density, bathymetry, and firn depth correction. Sadly, another considerable error source exists in the vicinity with grounding line where ice might not be fully hydrostatically at equilibrium. This issue of possible thickness errors near grounding line was not fully addressed in this study. Further work will need to be done to determine error of hydrostatic ice thickness and how this influences on grounding line estimation. Nevertheless, the results in this study confirm that DEMs derived from TanDEM-X data have sufficient spatial and vertical accuracy to allow accurate estimation and monitoring grounding line of fast-moving ice bodies; this approach should be utilized in other regions. １０７ 국문 요약문 빙하의 지반접지선은 빙붕의 안정성에 대한 중요한 정보를 제공하 고 빙상모델의 경계 조건으로도 사용된다. 따라서 접지선의 이동을 감시하는 것은 빙붕의 취약점을 검사할 수 있는 주요 척도에 속한다. 최근 다양한 위성 시스템을 이용한 기술들을 사용하여 지반접지선의 추정 및 모니터링이 이루어지고 있으며 서남극 대부분의 빙하와 몇 개의 동남극 빙하에서 지반접지선 후퇴가 진행되고 있는 것을 발견 하였다. 본 연구 또한 남극 지반접지선의 추정 및 변동을 감시하기 위해 이루어졌다. 그러나 본 연구의 가장 중요한 측면은 지반접지선의 변 동 감지는 높은 시간적 및 공간적 해상도의 모니터링이 필요하며 따 라서 단순하고 반복 가능한 기술이 필요하다는 점이다. 본 연구에서 는 빙하의 기저 변동이 얼음 표면에 전달되는 특징을 주목하여 빙하 표면의 고도 정보의 변화를 이용하는 것이 가장 적절한 접근 방법이 라고 판단하였다. 또한 최근의 SAR 위성 시스템의 개발 및 발전을 통해 고행상도 표면 고도 정보를 단시간에 얻을 수 있기 때문에 본 연구는 이러한 고해상도 표면 고도 정보를 유체정역학에 적용하여 계산한 얼음의 두께가 지반접지선의 위치를 추정하는데 유효하게 이 용될 수 있음을 보여주고자 하였다. １０８ 유체정역학적 얼음의 두께를 계산하기 위해서는 얼음과 바다 표면 고도의 독립적인 자료가 필요하며 얼음이 지반으로부터 떨어지면서 얼음과 해수의 밀도 차이로 인하여 떠오르기 시작하는 위치를 결정 하기 위해서는 해저 지형이 필요하다. 본 연구에서는 TanDEM-X 와 TerraSAR-X bi-static 간섭기법을 이용하여 눈과 얼음으로 덮여 있는 남극 빙하의 디지털 표고 모델을 높은 공간 해상도 및 수직 해 상도로 제작하였다. 얼음의 흐름과 강설로 인하여 지속적으로 변화 하는 남극의 성격 때문에 신뢰할 수 있는 지상기준고도 확보의 어려 움은 이러한 간섭기법을 활용한 표고 모델 제작에 있어서 가장 큰 문제였다. 이 문제를 해결하기 위해 SAR 위성자료가 확보된 시기와 가장 근접한 시기에 확보된 CryoSat-2 SARIn 레이더 고도계 데이 터를 상호 상관 값과 표준 편차를 이용하여 신뢰할 수 있는 자료만 을 신중하게 선택한 후 지상기준고도로 사용하였다. 바다 표면 고도 는 TanDEM-X와 TerraSAR-X 취득시의 고도를 추출하기 위해서 EIGEN-GL04C 지오이드와 FES2014 해양 조수 모델을 융합하여 추정하였다. 마지막으로, 해저 지형은 BEDMAP2 데이터를 사용하 였다. 위의 자료들을 시계열로 확보한 후 동남극의 Campbell 빙하 와 서남극의 Thwaites 빙하의 지반접지선을 추정하였으며 그 결과 접지선의 지속적인 관찰과 감시가 고해상도로 가능함을 확인할 수 있었다. 본 연구는 이전에 추정 된 지반접지선의 위치를 검증함으로써 기 １０９ 존에 알고 있는 정보를 확인할 뿐만 아니라 단주기의 자료 확보를 통한 지반접지선의 단기적 변동 관측과 고해상도 자료를 이용함에 따른 국지적 변동 관측을 가능하게 하였다. 이를 이용하여 얼음의 변화를 새로운 관점으로 이해할 수 있도록 도움을 줄 수 있었다. 그 예로 2011 년에서 2013 년 사이에 Thwaites Ice Shelf 에서 발견 된 ice rumple 에 대한 자세한 연구가 이루어졌다. 이 ice rumple 은 추정된 지반접지선에서 약 5km 떨어진 지역에 위치하고 있었으 며 2011 년에 처음 관측된 이후로 접지 반경이 점차 감소하는 것을 관찰할 수 있었다. 얼음의 두께가 서서히 감소하기 때문에 접지되어 있는 면적의 감소로 인하여 이러한 현상이 발생한다고 가정을 하고 점탄성 변형 모델 (viscoelastic deformation model) 을 사용하여 얼 음 두께의 진행상황을 해석하고자 하였다. 이러한 결과에도 불구하고 본 연구에서 제안된 기술의 주요 약점 은 결과의 정확도가 얼음의 평균 밀도와 사용되는 입력 데이터의 정 확도에 강하게 의존한다는 점으로 요약할 수 있다. 따라서 유체정역 학적 얼음두께의 추정 시 사용된 얼음의 평균 밀도는 정확하게 알려 져 있기 않기 때문에 밀도를 변화시키면서 지반접지선 추정이 어떻 게 변하는 지를 분석하였으며 현장 조사를 통하여 확보된 얼음의 밀 도를 이용하였을 때 기존에 알려진 지반접지선과 비교적 일치하는 것을 확인하였다. 또한 서로 다른 공간과 수직적 정확도를 갖는 해 저 지형 데이터 세트의 결과 비교하였을 때 더 나은 데이터를 사용 １１０ 한다면 결과의 정확도가 크게 개선되는 것을 확인할 수 있었다. Thwaites 빙하처럼 광범위한 현장 자료가 존재하는 지역에서는 보 다 정확한 지반접지선 추정 및 관측이 가능할 것으로 기대된다. 안타깝게도 본 연구는 얼음이 지반접지선 근처에서 유체정역학적 평형 상태에 있지 않기 때문에 발생하는 얼음 두께 오차의 문제는 제대로 연구가 진행되지 않았다. 따라서 얼음 평균 밀도에서 발생 하는 오차와 함께 가장 큰 오차 원인으로 고려된다. 유체정역학적 얼음 두께에서 발생하는 오차가 지반접지선 추정에 어떠한 영향을 어느 정도 미칠 수 있는지는 추가적인 연구가 필요할 것이며 이러한 오차를 줄여나간다면 지반접지선 추정의 정확도를 더욱 증가시킬 수 있을 것으로 고려된다. １１１ Bibliography Baek, S., O. I. Kwoun, A. Braun, Z. Lu and C. K. Shum (2005). "Digital elevation model of King Edward VII Peninsula, West Antarctica, from SAR interferometry and ICESat laser altimetry." IEEE Geoscience and Remote Sensing Letters 2(4): 413-417. Becek, K., W. Koppe and S. H. Kutoglu (2016). "Evaluation of Vertical Accuracy of the WorldDEM (TM) Using the Runway Method." Remote Sensing 8(11). Berger, S., L. Favier, R. Drews, J. J. Derwael and F. Pattyn (2016). "The control of an uncharted pinning point on the flow of an Antarctic ice shelf." Journal of Glaciology 62(231): 37-45. Bindschadler, R., H. Choi, A. Wichlacz, R. Bingham, J. Bohlander, K. Brunt, H. Corr, R. Drews, H. Fricker, M. Hall, R. Hindmarsh, J. Kohler, L. Padman, W. Rack, G. Rotschky, S. Urbini, P. Vornberger and N. Young (2011). "Getting around Antarctica: new high-resolution mappings of the grounded and freelyfloating boundaries of the Antarctic ice sheet created for the International Polar Year." Cryosphere 5(3): 569-588. Blankenship, D. D., S. D. Kempf, D. A. Young, J. L. Roberts, T. van Ommen, R. Forsberg, M. J. Siegert, S. J. Palmer and J. A. Dowdeswell (2012, updated 2013). "IceBridge Riegl Laser Altimeter L2 Geolocated Surface Elevation Triplets, Version 1." (Boulder, Colorado USA. NASA National Snow and Ice Data Center Distributed Active Archive Center.). Bohlander, J. and T. Scambos (2007). "Antarctic coastlines and grounding line derived from MODIS Mosaic of Antarctica (MOA), Boulder, Colorado １１２ USA: National Snow and Ice Data Center." Digital media. Brunt, K. M., H. A. Fricker, L. Padman, T. A. Scambos and S. O'Neel (2010). "Mapping the grounding zone of the Ross Ice Shelf, Antarctica, using ICESat laser altimetry." Annals of Glaciology 51(55): 71-79. Carrère, L., F. Lyard, M. Cancet, A. Guillot, N. Picot and S. Dupuy (2015). "FES 2014: a new global tidal model". Proceedings of the Ocean Surface Topography Science Team Meeting, Reston, FL, USA. Carrère, L., F. Lyard, M. Cancet, A. Guillot and L. Roblou (2013). "FES 2012: a new global tidal model taking advantage of nearly 20 years of altimetry". 20 Years of Progress in Radar Altimatry. Costantini, M. (1998). "A novel phase unwrapping method based on network programming." IEEE Transactions on Geoscience and Remote Sensing 36(3): 813-821. Davis, C. H. and A. C. Ferguson (2004). "Elevation change of the Antarctic ice sheet, 1995-2000, from ERS-2 satellite radar altimetry." IEEE Transactions on Geoscience and Remote Sensing 42(11): 2437-2445. De Rydt, J., G. H. Gudmundsson, T. Nagler, J. Wuite and E. C. King (2018). "Recent rift formation and impact on the structural integrity of the Brunt Ice Shelf, East Antarctica." Cryosphere 12(2): 505-520. DiMarzio, J., A. Brenner, R. Schutz, C. Shuman and H. Zwally (2007). "GLAS/ICESat 500 m laser altimetry digital elevation model of Antarctica." National Snow and Ice Data Center (NSIDC). １１３ Doake, C. S. M. and E. W. Wolff (1985). "Flow Law for Ice in Polar Ice Sheets." Nature 314(6008): 255-257. Drews, R. (2015). "Evolution of ice-shelf channels in Antarctic ice shelves." Cryosphere 9(3): 1169-1181. Favier, L., O. Gagliardini, G. Durand and T. Zwinger (2012). "A threedimensional full Stokes model of the grounding line dynamics: effect of a pinning point beneath the ice shelf." Cryosphere 6(1): 101-112. Favier, L. and F. Pattyn (2015). "Antarctic ice rise formation, evolution, and stability." Geophysical Research Letters 42(11): 4456-4463. Fialko, Y., Y. Khazan and M. Simons (2001). "Deformation due to a pressurized horizontal circular crack in an elastic half-space, with applications to volcano geodesy." Geophysical Journal International 146(1): 181-190. Forste, C., R. Schmidt, R. Stubenvoll, F. Flechtner, U. Meyer, R. Konig, H. Neumayer, R. Biancale, J. M. Lemoine, S. Bruinsma, S. Loyer, F. Barthelmes and S. Esselborn (2008). "The GeoForschungsZentrum Potsdam/Groupe de Recherche de Geodesie Spatiale satellite-only and combined gravity field models: EIGEN-GL04S1 and EIGEN-GL04C." Journal of Geodesy 82(6): 331346. Fricker, H. A., I. Allison, M. Craven, G. Hyland, A. Ruddell, N. Young, R. Coleman, M. King, K. Krebs and S. Popov (2002). "Redefinition of the Amery Ice Shelf, East Antarctica, grounding zone." Journal of Geophysical ResearchSolid Earth 107(B5). Fricker, H. A., R. Coleman, L. Padman, T. A. Scambos, J. Bohlander and K. １１４ M. Brunt (2009). "Mapping the grounding zone of the Amery Ice Shelf, East Antarctica using InSAR, MODIS and ICESat." Antarctic Science 21(5): 515532. Fricker, H. A. and L. Padman (2006). "Ice shelf grounding zone structure from ICESat laser altimetry." Geophysical Research Letters 33(15). Fricker, H. A., S. Popov, I. Allison and N. Young (2001). "Distribution of marine ice beneath the Amery Ice Shelf." Geophysical Research Letters 28(11): 2241-2244. Furst, J. J., G. Durand, F. Gillet-Chaulet, N. Merino, L. Tavard, J. Mouginot, N. Gourmelen and O. Gagliardini (2015). "Assimilation of Antarctic velocity observations provides evidence for uncharted pinning points." Cryosphere 9(4): 1427-1443. Goldstein, R. M. and C. L. Werner (1998). "Radar interferogram filtering for geophysical applications." Geophysical Research Letters 25(21): 4035-4038. Greenbaum, J. S., D. D. Blankenship, D. A. Young, T. G. Richter, J. L. Roberts, A. R. A. Aitken, B. Legresy, D. M. Schroeder, R. C. Warner, T. D. van Ommen and M. J. Siegert (2015). "Ocean access to a cavity beneath Totten Glacier in East Antarctica." Nature Geoscience 8(4): 294-298. Greene, C. A., D. E. Gwyther and D. D. Blankenship (2017). "Antarctic Mapping Tools for MATLAB." Computers & Geosciences 104: 151-157. Griggs, J. A. and J. L. Bamber (2011). "Antarctic ice-shelf thickness from satellite radar altimetry." Journal of Glaciology 57(203): 485-498. １１５ Groh, A., H. Ewert, R. Rosenau, E. Fagiolini, C. Gruber, D. Floricioiu, W. A. Jaber, S. Linow, F. Flechtner, M. Eineder, W. Dierking and R. Dietrich (2014). "Mass, Volume and Velocity of the Antarctic Ice Sheet: Present-Day Changes and Error Effects." Surveys in Geophysics 35(6): 1481-1505. Gudmundsson, G. H., C. F. Raymond and R. Bindschadler (1998). "The origin and longevity of flow stripes on Antarctic ice streams." Annals of Glaciology, Vol 27, 1998 27: 145-152. Ha, H. K., A. K. Wahlin, T. W. Kim, S. H. Lee, J. H. Lee, H. J. Lee, C. S. Hong, L. Arneborg, G. Bjork and O. Kalen (2014). "Circulation and Modification of Warm Deep Water on the Central Amundsen Shelf." Journal of Physical Oceanography 44(5): 1493-1501. Han, H., Y. Ji and H. Lee (2013). "Estimation of annual variation of ice extent and flow velocity of Campbell Glacier in East Antarctica using COSMOSkyMed SAR images." Korean Journal of Remote Sensing 29(1): 45-55. Han, H. and H. Lee (2014). "Tide deflection of Campbell Glacier Tongue, Antarctica, analyzed by double-differential SAR interferometry and finite element method." Remote Sensing of Environment 141: 201-213. Han, H. and H. Lee (2015). "Tide-corrected flow velocity and mass balance of Campbell Glacier Tongue, East Antarctica, derived from interferometric SAR." Remote Sensing of Environment 160: 180-192. Helm, V., A. Humbert and H. Miller (2014). "Elevation and elevation change of Greenland and Antarctica derived from CryoSat-2." Cryosphere 8(4): 15391559. １１６ Jacobs, S. S., H. H. Helmer, C. S. M. Doake, A. Jenkins and R. M. Frolich (1992). "Melting of Ice Shelves and the Mass Balance of Antarctica." Journal of Glaciology 38(130): 375-387. Jenkins, A., P. Dutrieux, S. S. Jacobs, S. D. McPhail, J. R. Perrett, A. T. Webb and D. White (2010). "Observations beneath Pine Island Glacier in West Antarctica and implications for its retreat." Nature Geoscience 3(7): 468-472. Jonathan, B. (1994). "A comparison of satellite-altimetry and ice-thickness measurements of the Ross Ice Shelf, Antarctica." Annals of Glaciology 20: 357364. Joughin, I., B. E. Smith and B. Medley (2014). "Marine Ice Sheet Collapse Potentially Under Way for the Thwaites Glacier Basin, West Antarctica." Science 344(6185): 735-738. Kim, J. W., D. J. Kim, S. H. Kim, H. K. Ha and S. H. Lee (2015). "Disintegration and acceleration of Thwaites Ice Shelf on the Amundsen Sea revealed from remote sensing measurements." Giscience & Remote Sensing 52(4): 498-509. Kim, S. H. and D. J. Kim (2017). "Combined Usage of TanDEM-X and CryoSat-2 for Generating a High Resolution Digital Elevation Model of Fast Moving Ice Stream and Its Application in Grounding Line Estimation." Remote Sensing 9(2). Konrad, H., A. Shepherd, L. Gilbert, A. E. Hogg, M. McMillan, A. Muir and T. Slater (2018). "Net retreat of Antarctic glacier grounding lines." Nature Geoscience 11(4): 258-+. １１７ Krabill, W. B., W. Abdalati, E. B. Frederick, S. S. Manizade, C. F. Martin, J. G. Sonntag, R. N. Swift, R. H. Thomas and J. G. Yungel (2002). "Aircraft laser altimetry measurement of elevation changes of the greenland ice sheet: technique and accuracy assessment." Journal of Geodynamics 34(3-4): 357-376. Krieger, G., A. Moreira, H. Fiedler, I. Hajnsek, M. Werner, M. Younis and M. Zink (2007). "TanDEM-X: A satellite formation for high-resolution SAR interferometry." IEEE Transactions on Geoscience and Remote Sensing 45(11): 3317-3341. Krieger, G., M. Zink, M. Bachmann, B. Brautigam, D. Schulze, M. Martone, P. Rizzoli, U. Steinbrecher, J. W. Antony, F. De Zan, I. Hajnsek, K. Papathanassiou, F. Kugler, M. R. Cassola, M. Younis, S. Baumgartner, P. Lopez-Dekker, P. Prats and A. Moreira (2013). "TanDEM-X: A radar interferometer with two formation-flying satellites." Acta Astronautica 89: 8398. Legresy, B., F. Papa, F. Remy, G. Vinay, M. van den Bosch and O. Z. Zanife (2005). "ENVISAT radar altimeter measurements over continental surfaces and ice caps using the ICE-2 retracking algorithm." Remote Sensing of Environment 95(2): 150-163. Legresy, B., F. Remy and F. Blarel (2006). "Along track repeat altimetry for ice sheets and continental surface studies." In Proc. of the Symposium on 15 Years of Progress in Radar Altimetry. Leuschen, C., P. Gogineni, F. Rodriguez-Morales, J. Paden and C. Allen (2010, updated 2018). "IceBridge MCoRDS L2 Ice Thickness, Version 1." Boulder, Colorado USA. NASA National Snow and Ice Data Center Distributed Active Archive Center. １１８ Lyard, F., F. Lefevre, T. Letellier and O. Francis (2006). "Modelling the global ocean tides: modern insights from FES2004." Ocean Dynamics 56(5-6): 394-415. Macayeal, D. R., R. A. Bindschadler, S. Shabtaie, S. Stephenson and C. R. Bentley (1987). "Force, Mass, and Energy Budgets of the Crary Ice Rise Complex, Antarctica." Journal of Glaciology 33(114): 218-230. MacGregor, J. A., G. A. Catania, M. S. Markowski and A. G. Andrews (2012). "Widespread rifting and retreat of ice-shelf margins in the eastern Amundsen Sea Embayment between 1972 and 2011." Journal of Glaciology 58(209): 458-466. Matsuoka, K., R. C. A. Hindmarsh, G. Moholdt, M. J. Bentley, H. D. Pritchard, J. Brown, H. Conway, R. Drews, G. Durand, D. Goldberg, T. Hattermann, J. Kingslake, J. T. M. Lenaerts, C. Martin, R. Mulvaney, K. W. Nicholls, F. Pattyn, N. Ross, T. Scambos and P. L. Whitehouse (2015). "Antarctic ice rises and rumples: Their properties and significance for ice-sheet dynamics and evolution." Earth-Science Reviews 150: 724-745. Mctigue, D. F. (1987). "Elastic Stress and Deformation near a Finite Spherical Magma Body - Resolution of the Point-Source Paradox." Journal of Geophysical Research-Solid Earth and Planets 92(B12): 12931-12940. Mercer, J. H. (1978). "West Antarctic Ice Sheet and Co2 Greenhouse Effect - Threat of Disaster." Nature 271(5643): 321-325. Mercer, J. H. (1978). "West Antarctic ice sheet and CO 2 greenhouse effectA threat of disaster." Nature 271(5643): 321-325. １１９ Mogi, K. (1958). "Relations between the eruptions of various volcanoes and the deformations of the ground surfaces around them." Bull. Earthq. Res. Inst. 36: 99-134. Okada, Y. (1985). "Surface Deformation Due to Shear and Tensile Faults in a Half-Space." Bulletin of the Seismological Society of America 75(4): 11351154. Padman, L. and H. A. Fricker (2005). "Tides on the Ross Ice Shelf observed with ICESat." Geophysical Research Letters 32(14). Paolo, F. S., H. A. Fricker and L. Padman (2015). "Volume loss from Antarctic ice shelves is accelerating." Science 348(6232): 327-331. Peltier, R. (1982). Dynamics of the ice age Earth. Advances in geophysics, Elsevier. 24: 1-146. Perna, S., C. Esposito, P. Berardino, A. Pauciullo, C. Wimmer and R. Lanari (2015). "Phase Offset Calculation for Airborne InSAR DEM Generation Without Corner Reflectors." IEEE Transactions on Geoscience and Remote Sensing 53(5): 2713-2726. Rémy, F. and S. Parouty (2009). "Antarctic ice sheet and radar altimetry: A review." Remote Sensing 1(4): 1212-1239. Raney, R. K. (1998). "The delay/Doppler radar altimeter." IEEE Transactions on Geoscience and Remote Sensing 36(5): 1578-1588. Rignot, E. (2001). "Evidence for rapid retreat and mass loss of Thwaites １２０ Glacier, West Antarctica." Journal of Glaciology 47(157): 213-222. Rignot, E. (2006). "Changes in ice dynamics and mass balance of the Antarctic ice sheet." Philosophical Transactions of the Royal Society aMathematical Physical and Engineering Sciences 364(1844): 1637-1655. Rignot, E., S. Jacobs, J. Mouginot and B. Scheuchl (2013). "Ice-Shelf Melting Around Antarctica." Science 341(6143): 266-270. Rignot, E. and S. S. Jacobs (2002). "Rapid bottom melting widespread near Antarctic ice sheet grounding lines." Science 296(5575): 2020-2023. Rignot, E., J. Mouginot, M. Morlighem, H. Seroussi and B. Scheuchl (2014). "Widespread, rapid grounding line retreat of Pine Island, Thwaites, Smith, and Kohler glaciers, West Antarctica, from 1992 to 2011." Geophysical Research Letters 41(10): 3502-3509. Rignot, E., J. Mouginot and B. Scheuchl (2011). "Antarctic grounding line mapping from differential satellite radar interferometry." Geophysical Research Letters 38. Rignot, E., J. Mouginot and B. Scheuchl (2011). "Ice Flow of the Antarctic Ice Sheet." Science 333(6048): 1427-1430. Rignot, E., D. G. Vaughan, M. Schmeltz, T. Dupont and D. MacAyeal (2002). "Acceleration of Pine Island and Thwaites Glaciers, West Antarctica." Annals of Glaciology, Vol 34, 2002 34: 189-194. Rossi, C., F. R. Gonzalez, T. Fritz, N. Yague-Martinez and M. Eineder (2012). １２１ "TanDEM-X calibrated Raw DEM generation." ISPRS Journal of Photogrammetry and Remote Sensing 73: 12-20. Sandwell, D. T. and M. B. Ruiz (1992). "Along-Track Gravity-Anomalies from Geostat and Seasat Altimetry - Gebco Overlays." Marine Geophysical Researches 14(3): 165-205. Scambos, T., T. Haran, M. Fahnestock, T. Painter and J. Bohlander (2007). "MODIS-based Mosaic of Antarctica (MOA) data sets: Continent-wide surface morphology and snow grain size." Remote Sensing of Environment 111(2-3): 242-257. Schmeltz, M., E. Rignot and D. MacAyeal (2002). "Tidal flexure along icesheet margins: comparison of InSAR with an elastic-plate model." Annals of Glaciology, Vol 34, 2002 34: 202-208. Schmeltz, M., E. Rignot and D. R. MacAyeal (2001). "Ephemeral grounding as a signal of ice-shelf change." Journal of Glaciology 47(156): 71-77. Shepherd, A. and D. Wingham (2007). "Recent sea-level contributions of the Antarctic and Greenland ice sheets." science 315(5818): 1529-1532. Stammer, D., R. D. Ray, O. B. Andersen, B. K. Arbic, W. Bosch, L. Carrere, Y. Cheng, D. S. Chinn, B. D. Dushaw, G. D. Egbert, S. Y. Erofeeva, H. S. Fok, J. A. M. Green, S. Griffiths, M. A. King, V. Lapin, F. G. Lemoine, S. B. Luthcke, F. Lyard, J. Morison, M. Muller, L. Padman, J. G. Richman, J. F. Shriver, C. K. Shum, E. Taguchi and Y. Yi (2014). "Accuracy assessment of global barotropic ocean tide models." Reviews of Geophysics 52(3): 243-282. Tinto, K. J. and R. E. Bell (2011). "Progressive unpinning of Thwaites １２２ Glacier from newly identified offshore ridge: Constraints from aerogravity." Geophysical Research Letters 38. Van den Breke, M., W. J. Van de Berg and E. Van Meijgaard (2008). "Firn depth correction along the Antarctic grounding line." Antarctic Science 20(5): 513-517. Vaughan, D. G. (1995). "Tidal Flexure at Ice Shell Margins." Journal of Geophysical Research-Solid Earth 100(B4): 6213-6224. Wang, F., J. L. Bamber and X. Cheng (2015). "Accuracy and Performance of CryoSat-2 SARIn Mode Data Over Antarctica." IEEE Geoscience and Remote Sensing Letters 12(7): 1516-1520. Wen, J. H., Y. F. Wang, W. L. Wang, K. C. Jezek, H. X. Liu and I. Allison (2010). "Basal melting and freezing under the Amery Ice Shelf, East Antarctica." Journal of Glaciology 56(195): 81-90. Wild, C. T., O. J. Marsh and W. Rack (2017). "Viscosity and elasticity: a model intercomparison of ice-shelf bending in an Antarctic grounding zone." Journal of Glaciology 63(240): 573-580. Williams, C. A. and G. Wadge (1998). "The effects of topography on magma chamber deformation models: Application to Mt Etna and radar interferometry." Geophysical Research Letters 25(10): 1549-1552. Wingham, D. (2005). "CryoSat: A mission to the ice fields of Earth." Esa Bulletin-European Space Agency(122): 10-17. １２３ Wingham, D., C. Francis, S. Baker, C. Bouzinac, D. Brockley, R. Cullen, P. de Chateau-Thierry, S. Laxon, U. Mallow and C. Mavrocordatos (2006). "CryoSat: A mission to determine the fluctuations in Earth’s land and marine ice fields." Advances in Space Research 37(4): 841-871. Zhou, Y., C. X. Zhou, F. H. Deng, D. C. E, H. Y. Liu and Y. M. Wen (2015). "Improving InSAR elevation models in Antarctica using laser altimetry, accounting for ice motion, orbital errors and atmospheric delays." Remote Sensing of Environment 162: 112-118. １２４