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. In this study, we used TanDEM-X data to construct high-resolution DEMs of the Thwaites ice shelf in West Antarctica from 2011 to 2013. We also generated surface deformation maps which allowed us to detect and monitor the elevation changes of an ice rumple that appeared sometime between the observations of a grounding line of the Thwaites glacier using Double-Differential Interferometric SAR (DDInSAR) in 1996 and 2011. The observed degradation of the ice rumple during 2011–2013 may be related 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. Our conclusion is 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 rumples, with implications for monitoring ice shelf stability.


Introduction
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 [1,2].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 [3].Ice rumples are typically elevated only tens of meters above the surrounding ice shelf surface [2] and can vary in size.For example, the Doake ice rumples in the Ronne Ice Shelf are relatively large with an area of 1700 km 2 [4], while smaller ice rumples may only be a few kilometers wide in other ice shelves [5].
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. [6] demonstrated that an ice rumple can reduce shelf velocity and trigger the advance of the grounding line over the long-term increasing stability.Conversely, Tinto and Bell [7] 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 [8][9][10][11], field observations are still difficult to obtain due to the presence of crevasses and rifts [2].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 Interferometric SAR (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 digital elevation models (DEMs) of the globe with absolute vertical accuracy of better than 10 m and a spatial resolution of 12 × 12 m [12].This enables the ice rumples' topographic features to be observed, as well as the monitoring of changes when data are available in time series.
In this study, we used multiple 2011-2013 TanDEM-X datasets to study an ice rumple in Antarctica's Thwaites ice shelf, previously observed via DDInSAR by Rignot et al. [13].The pinning point appeared sometime between the initial observation in 1996 and the last observation in 2011 by DDInSAR.Subsequently, we compared each DEM, derived from TanDEM-X datasets, in an original way by anticipating surface deformation of the ice rumple, because the Thwaites ice shelf has recently been experiencing rapid thinning [14].Finally, we investigated the surface response due to short-term progressive basal stress perturbation using a simple deformation model, assuming that the extent of grounding underneath the ice acted as a spherical stress source.

Study Location
The ice rumple studied here was located on the Thwaites ice shelf in West Antarctica, centered approximately 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 location was approximately 5 km downstream of the grounding line estimated in 2011 by Rignot et al. [13] and had not moved since its first observation in the same study (Figure 1).This ice rumple was not observable in 1996 because the grounding line was located considerably 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 [7].
Thwaites ice shelf is the seaward extension of Thwaites Glacier, which is one of the major glaciers in West Antarctica exhibiting a significant retreat in grounding line [13] and ice shelf margins [18], negative mass balance, and ice volume loss [14] in recent decades.The intrusion of the circumpolar deep water into the ocean base underneath the glaciers has been considered to be the predominant cause of such immense changes.The ice shelf has been flowing at high velocity, exceeding 4 km per year at the fastest.It is also known for its strong katabatic wind and rapid and heavy precipitation, which can cause dramatic surface changes at any time.In the grounding lines' near vicinity, the surface elevation can range approximately from 50 to 300 m above the WGS84 ellipsoid.
approximately at 75°23′ S and 106°42′ W with a radius of 1 km in both the longitudinal and alongflow direction in 2011.The location was approximately 5 km downstream of the grounding line estimated in 2011 by Rignot et al. [13] and had not moved since its first observation in the same study (Figure 1).This ice rumple was not observable in 1996 because the grounding line was located considerably 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 [7].

High-Resolution Deformation Map
We 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 by Kim and Kim [19].TanDEM-X data were acquired on 10 June 2012, 29 June 2012, and 27 June 2013, allowing us to generate two sets of deformation maps.The timespans between each pair of DEMs 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.We followed the methods described by Kim and Kim [19], where near-coincidental CryoSat-2 radar altimeter datasets are able to 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.We 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 × 25 m grid size from raw data with a spatial resolution of ~3 m along track and ~1.2 m in radar line of sight [12] and a vertical resolution of better than 2 m over smooth areas [20].The elevation was defined and geocoded with respect to the WGS84 ellipsoid.Then, a 2D Gaussian smoothing kernel with a standard deviation of 6 was applied to diminish the effect of the crevasses, which can otherwise contaminate the DEM difference calculations, and to extract distinct surface elevations of the ice rumple.
Elevation change was calculated by simply subtracting geocoded DEMs in UTM projection.Geocoding errors were considered negligible because the flight information of the TanDEM-X is known to be precise, and geocoding was done using this flight information.Horizontal and vertical accuracy was considered to be satisfied because each DEM was spatially and vertically compared with the CryoSat-2 SARIn points.The overall difference in elevation contained the surface elevation change proportional to the thickness change over the study period as well as the vertical offset due to the mismatched tidal heights.Thus, these tidal differences had to be eliminated in order to accurately extract the surface deformation, which was solely due to changes in the ice rumple.Tidal heights were estimated using the FES2014 tidal model, produced by the NOVELTIS, Laboratoire d'Etudes en Géophysique et Océanographie Spatiales (LEGOS), and Collecte Localisation Satellites (CLS) Space Oceanography Division and distributed by Archiving, Validation and Interpretation of Satellite Oceanographic data (AVISO), with support from the Centre national d'études spatiales (CNES) (http://www.aviso.altimetry.fr/).Tidal 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.

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 may occur when there are changes to the extent and/or center location of grounding underneath the ice.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 2).
In such cases, the point dilatation model described by Mogi [21] is the most well-known and has been widely used in studies of volcanic magma chambers.McTigue [22] derived a more elaborate expression to analytically approximate deformation that included 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, we considered the simple elastic deformation model to be capable of providing good approximations of deformation resulting from pressure changes within the shallow ice.
We assumed the linear viscoelasticity of the ice by considering a combination of elastic and viscous behavior, based on the conclusion of Doake and Wolff [23] that Antarctic ice shelves which may be influenced by the restraining force from ice rumples display a linear flow characteristic, even though conventional flow laws for ice consider non-linear stress-dependent viscosity.In addition, Gudmundsson et al. [24] 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.We assumed that the ice shelf was sliding rapidly, as the friction underneath changed from positive to zero when moving seaward from the grounding line.
In addition, we 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 it with 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, we implemented the effective shear modulus for a standard linear solid [25].
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 it with 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, we implemented the effective shear modulus for a standard linear solid [25].The elastic deformation at the free surface (z = 0) in the vertical direction (u z ) due to a pressure change (∆P) within a spherical source located at x 0 , y 0 , and z 0 is given as follows [22]: where r is the radial distance from the free surface above the center of a spherical source, z 0 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.We used values for the latter two (0.3 and 0.88 GPa, respectively) from Vaughan [4], although high variability in an effective Young's modulus has been observed depending on the location and time scale of the ice deformation [26].The deformation here was normalized by ∆Pz 0 .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 a time-domain has been found by Peltier [25] and is given by: where η is viscosity and s is the Laplace transform variable.
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: , and where t is the duration time for viscous deformation.We used a viscosity of 50.1 TPa [27].
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.We expected that having a 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) [28].

Optimization
Basically, this is the inverse problem where the left-hand side of the Equation ( 3) is known and the model parameters are unknown.In order to determine the depth (z 0 ), radius (α), and pressure change (∆P) of a spherical pressure source beneath the ice rumple as well as its center (x 0 , y 0 ), we 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.We adopted the interior-point algorithm with randomly-selected initial values within the boundary conditions (Table 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 was the RMS difference between the observed and estimated deformation.Five unknown variables, and the minimum and maximum constraints of each, are listed in Table 1.A total of 50 optimization processes were repeated for measuring the stability of results.

Results
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 3).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 (shown by dark grey circles in Figure 3).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 (shown by yellow circles in Figure 3).Using these selected elevation data as ground control points, the correction for absolute elevation was carried out.The root-mean-square errors (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 3).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 (shown by dark grey circles in Figure 3).The selection process then sorted out 22, 31, and 39 points with mean crosscorrelation coefficients of 0.97, 9.98, and 0.98, respectively (shown by yellow circles in Figure 3).Using these selected elevation data as ground control points, the correction for absolute elevation was carried out.The root-mean-square errors (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.The ice rumple had a circular shape with a downstream slope steeper than that upstream (Figure 4a-c), resembling the simulation results of Gudmundsson et al. [24] showing that surface undulations are in phase with basal perturbations under a high basal slip ratio.The elevation anomaly from the surrounding region was approximately 35 m, supporting its characterization as an ice rumple.It The ice rumple had a circular shape with a downstream slope steeper than that upstream (Figure 4a-c), resembling the simulation results of Gudmundsson et al. [24] showing that surface undulations are in phase with basal perturbations under a high basal slip ratio.The elevation anomaly from the surrounding region was approximately 35 m, supporting its characterization as an ice rumple.The optimization searches the best approximation of parameters so that the modified deformation model can reproduce the near-circular pattern observed in the surface deformation for both periods.The optimization result predicted a pressure source at a shallow depth of 1072.42 ± 0.93 m, with a radius of 360.78 ± 15.16 m and pressure increment of −0.86 ± 0.16 GPa during 2011-2012 (Figure 5).The minimum RMS difference between the observed and estimated deformation was 1.13 m.Likewise, the result predicted a pressure source at a depth of 832.03 ± 1.04 m, with a radius of 148.27 ± 13.51 m and pressure increment of −8.59 ± 0.38 GPa during 2012-2013 (Figure 6).The RMS The optimization searches the best approximation of parameters so that the modified deformation model can reproduce the near-circular pattern observed in the surface deformation for both periods.The optimization result predicted a pressure source at a shallow depth of 1072.42 ± 0.93 m, with a radius of 360.78 ± 15.16 m and pressure increment of −0.86 ± 0.16 GPa during 2011-2012 (Figure 5).The minimum RMS difference between the observed and estimated deformation was 1.13 m.Likewise, the result predicted a pressure source at a depth of 832.03 ± 1.04 m, with a radius of 148.27 ± 13.51 m and pressure increment of −8.59 ± 0.38 GPa during 2012-2013 (Figure 6).The RMS difference between the observed and estimated deformation during this period was 2.11 m.During both periods, there was a slight offset of the pressure source's center location compared to the peak of observed deformation.

Discussion
The instant deformation map alone cannot uniquely resolve whether such a 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

Discussion
The instant deformation map alone cannot uniquely resolve whether such a 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

Discussion
The instant deformation map alone cannot uniquely resolve whether such a 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 was 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, we monitored the surface features of the ice rumple using Landsat 7 ETM+ images from 2003-2014 (Figure 7).The images in this time series showed a gradual dissipation of the ice rumple, strongly indicating continuous thinning of the Thwaites ice shelf.Furthermore, 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.Nonetheless, the ice shelf might have been in contact with the pinning point even after the disappearance of surface features, as intermittent ice contact to the pinning point could be possible due to ice shelf thickness fluctuations [29].
Remote Sens. 2018, 9, x FOR PEER REVIEW 8 of 14 was 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, we monitored the surface features of the ice rumple using Landsat 7 ETM+ images from 2003-2014 (Figure 7).The images in this time series showed a gradual dissipation of the ice rumple, strongly indicating continuous thinning of the Thwaites ice shelf.Furthermore, 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.Nonetheless, the ice shelf might have been in contact with the pinning point even after the disappearance of surface features, as intermittent ice contact to the pinning point could be possible due to ice shelf thickness fluctuations [29].Furthermore, the fact that the surface expression was not advected downstream during the observation period suggests that the ice shelf was still in contact with the basal topography.Although the peak of the surface topography did advance downstream by hundreds of meters each year, the surface expression would have already disappeared from the defined map window if the ice shelf was ungrounded.This is due to the overall flow velocity in this region which is nearly 3 km/year.The advance of the peaks was probably due to unknown irregular shapes of the basal topography in contact with the ice, producing spatially variable pressure as the ice thinned.The differences in the elevation and deformation peaks are related to the asymmetric shape of the ice rumple, as well as the nearly circular shape of the deformation pattern (Figure 4d,e).Surface undulations were also observable on the surface in the downstream from the ice rumple in Figure 7.These undulations advect downstream at the flow velocity of the ice shelf and cause positive (or negative) deformations when they are not aligned between the observations, as shown in Figure 4d,e.An elliptical deformation pattern in Figure 4e was likely due to this surface undulation.
As the linear theory is invalid under high stress, 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 have affected the ice in the early years (e.g., 2003 in Figure 7).However, the magnified images show that the surface crevasses present Furthermore, the fact that the surface expression was not advected downstream during the observation period suggests that the ice shelf was still in contact with the basal topography.Although the peak of the surface topography did advance downstream by hundreds of meters each year, the surface expression would have already disappeared from the defined map window if the ice shelf was ungrounded.This is due to the overall flow velocity in this region which is nearly 3 km/year.The advance of the peaks was probably due to unknown irregular shapes of the basal topography in contact with the ice, producing spatially variable pressure as the ice thinned.The differences in the elevation and deformation peaks are related to the asymmetric shape of the ice rumple, as well as the nearly circular shape of the deformation pattern (Figure 4d,e).Surface undulations were also observable on the surface in the downstream from the ice rumple in Figure 7.These undulations advect downstream at the flow velocity of the ice shelf and cause positive (or negative) deformations when they are not aligned between the observations, as shown in Figure 4d,e.An elliptical deformation pattern in Figure 4e was likely due to this surface undulation.As the linear theory is invalid under high stress, 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 have affected the ice in the early years (e.g., 2003 in Figure 7).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, we 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.Nevertheless, the interpretation may be more accurate when a non-linear viscoelasticity is employed, although it is very likely that the problem can be especially difficult to solve.It is expected that further work beyond the scope of linear viscoelasticity is required to advance our assumption.
Although a low RMS difference is not necessarily an indication of the validity of the model, we were able to fit the observed deformation with reasonable numbers.Nevertheless, this simple model alone cannot determine the actual shape of what is certainly a more complicated bathymetry, and we can only postulate the extent of grounding using the radius of a pressure source.In this context, we suggest that the initial ice thickness can be defined as the difference between depth and radius, as shown in Figure 8. Thickness change can then be derived as the variation of thicknesses at different time durations, though 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.
Remote Sens. 2018, 9, x FOR PEER REVIEW 9 of 14 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, we 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.Nevertheless, the interpretation may be more accurate when a non-linear viscoelasticity is employed, although it is very likely that the problem can be especially difficult to solve.It is expected that further work beyond the scope of linear viscoelasticity is required to advance our assumption.Although a low RMS difference is not necessarily an indication of the validity of the model, we were able to fit the observed deformation with reasonable numbers.Nevertheless, this simple model alone cannot determine the actual shape of what is certainly a more complicated bathymetry, and we can only postulate the extent of grounding using the radius of a pressure source.In this context, we suggest that the initial ice thickness can be defined as the difference between depth and radius, as shown in Figure 8. Thickness change can then be derived as the variation of thicknesses at different time durations, though 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.According to the optimization result, the ice thickness was 711.64 ± 14.25 m and 683.76 ± 12.48 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 [7].This indicates a thickness decrease of 36.17 ± 17.27 m during that one-year period.However, as the center of the pressure source in 2012 was located 341.96 m upstream and 210.10 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 that 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. [30] and Paolo et al. [14] reported that the thinning rate of the Thwaites ice shelf was 6.13 and 2.80 m/year, respectively.A major factor contributing to these differences may be related to the extent of the survey area, which was over 4000 km 2 (the entire ice shelf) for the other studies as opposed to a few km 2 focusing on a single ice rumple in our study.Besides, strong circulation near the grounding According to the optimization result, the ice thickness was 711.64 ± 14.25 m and 683.76 ± 12.48 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 [7].This indicates a thickness decrease of 36.17 ± 17.27 m during that one-year period.However, as the center of the pressure source in 2012 was located 341.96 m upstream and 210.10 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 that 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. [30] and Paolo et al. [14] reported that the thinning rate of the Thwaites ice shelf was 6.13 and 2.80 m/year, respectively.A major factor contributing to these differences may be related to the extent of the survey area, which was over 4000 km 2 (the entire ice shelf) for the other studies as opposed to a few km 2 focusing on a single ice rumple in our study.Besides, strong circulation near the grounding area might have caused this large localized basal melting.Thus, it should not be assumed that our results represent the overall thinning of the entire 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 [7], meaning it is important to realize when this has happened.Based on the Landsat imagery time series in Figure 7, un-grounding might have taken place sometime after 2013 due to the disappearance of the surface features.It is unclear at this point whether progressive unpinning of the ice rumple caused an increase in ice velocities of the Thwaites ice shelf.However, we suspect that the ice rumple investigated in this study did not have a significant impact on the stability of the 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 1).However, there are still some discrepancies, such as the modelled and measured ice thicknesses.Uncertainties within Young's modulus and viscosity are possible reasons for such deviations because we only adopted the rheology of the ice shelf based on previous research.It is suspected that the used ice rheology might be different from the actual state of the Thwaites ice shelf or is spatially varying, especially at the exact location of local grounding.Nevertheless, it is still 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 [31], and disks [32]) could be used to fit the observed deformation.
In order to analyze the dominant factor in the optimization result, we conducted simulations of the deformation model by varying the estimated parameters around a nominal value (Figure 9).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 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 deformations, 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.
area might have caused this large localized basal melting.Thus, it should not be assumed that our results represent the overall thinning of the entire ice shelf.
Despite this narrow focus on a local change, such progressive thinning may cause the ungrounding of a pinning point and the acceleration of the entire ice shelf [7], meaning it is important to realize when this has happened.Based on the Landsat imagery time series in Figure 7, ungrounding might have taken place sometime after 2013 due to the disappearance of the surface features.It is unclear at this point whether progressive unpinning of the ice rumple caused an increase in ice velocities of the Thwaites ice shelf.However, we suspect that the ice rumple investigated in this study did not have a significant impact on the stability of the 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 1).However, there are still some discrepancies, such as the modelled and measured ice thicknesses.Uncertainties within Young's modulus and viscosity are possible reasons for such deviations because we only adopted the rheology of the ice shelf based on previous research.It is suspected that the used ice rheology might be different from the actual state of the Thwaites ice shelf or is spatially varying, especially at the exact location of local grounding.Nevertheless, it is still 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 [31], and disks [32]) could be used to fit the observed deformation.
In order to analyze the dominant factor in the optimization result, we conducted simulations of the deformation model by varying the estimated parameters around a nominal value (Figure 9).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 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 deformations, 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.

Conclusions
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 monitoring such features in other ice shelves.
From 2011-2013, our deformation maps showed the recent fading of a small ice rumple in the surface of the Thwaites ice shelf, West Antarctica.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.The deformation pattern we found, 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 is likely to have since been unpinned.We then used the deformation maps 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.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.
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 km 2 ) compared to other research in this area.Thus, our results may not represent an overall thinning process, but do suggest that significantly 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.Our results confirm that DEMs derived from TanDEM-X data have sufficient spatial and vertical accuracy to allow accurate monitoring of fast-moving ice bodies; this approach should be utilized in other regions.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 still await to be answered regarding this process, more research using an increased number of datasets at different regions will be required.

Figure 1 .Figure 1 .
Figure 1.Grounding lines of Thwaites Glacier in 1996 (green) and 2011 (red) estimated using the DDInSAR method with European remote sensing (ERS) satellites[13,15].The orange dotted rectangle Figure 1.Grounding lines of Thwaites Glacier in 1996 (green) and 2011 (red) estimated using the DDInSAR method with European remote sensing (ERS) satellites [13,15].The orange dotted rectangle in the eastern shelf indicates a larger ice rumple previously discussed by Tinto and Bell [7].A newly-generated digital elevation model derived from TanDEM-X data on 10 June 2011 is shown within the yellow rectangle.The small red feature inside the yellow dotted square indicates the smaller ice rumple considered in this study.The background image is the MODIS Mosaic of Antarctica (MOA) image map [16].The overlaid ice velocity map was extracted from Rignot et al. [17].

Figure 2 .
Figure 2. Schematic drawings of (a) an ice rumple (with the extent of grounding represented by an idealized spherical pressure source) and (b) the elastic deformation model with a finite spherical pressure source.The initial surface is not flat, thus the radial distance from the center of a sphere to the surface varies.

Figure 2 .
Figure 2. Schematic drawings of (a) an ice rumple (with the extent of grounding represented by an idealized spherical pressure source) and (b) the elastic deformation model with a finite spherical pressure source.The initial surface is not flat, thus the radial distance from the center of a sphere to the surface varies.

Figure 3 .
Figure 3. Generated high-resolution surface topography maps of the Thwaites ice shelf on (a) 10 June 2011; (b) 29 June 2012; and (c) 24 June 2013.Dark grey areas denote the locations of all the available CryoSat-2 SARIn altimeter elevation data, and red circles denote selected data.

Figure 3 .
Figure 3. Generated high-resolution surface topography maps of the Thwaites ice shelf on (a) 10 June 2011; (b) 29 June 2012; and (c) 24 June 2013.Dark grey areas denote the locations of all the available CryoSat-2 SARIn altimeter elevation data, and red circles denote selected data.
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 periods of 2011-2012 and 2012-2013, respectively.Both deformation maps show surface depressions with a circular main shape and a weak elliptical signal in the cross-flow direction (Figure 4d,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 periods of 2011-2012 and 2012-2013, respectively.Remote Sens. 2018, 9, x FOR PEER REVIEW 6 of 14 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 periods of 2011-2012 and 2012-2013, respectively.Both deformation maps show surface depressions with a circular main shape and a weak elliptical signal in the cross-flow direction (Figure 4d,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 periods of 2011-2012 and 2012-2013, respectively.

Figure 4 .
Figure 4. Detailed surface digital elevation models (DEMs) of the ice rumple from 2011-2013, showing 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.The InSAR grounding line in 2011 is plotted by the yellow dotted line [13,15].

Figure 4 .
Figure 4. Detailed surface digital elevation models (DEMs) of the ice rumple from 2011-2013, showing 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.The InSAR grounding line in 2011 is plotted by the yellow dotted line [13,15].

Figure 5 .
Figure 5.Comparison between the observed and estimated deformation from the modified deformation model after the global optimization process during 2011-2012, with the setting time as the actual elapsed time.The minimum root-mean-square (RMS) difference between the two deformation maps was 1.20 m.Profiles in the longitudinal and along-flow directions are plotted along the red dashed and dot-dashed lines, respectively.

Figure 6 .
Figure 6.Comparison between the observed and estimated deformations from the modified deformation model after the global optimization process during 2012-2013, with the setting time as the actual elapsed time.The minimum RMS difference between the two deformation maps was 2.11 m.Profiles in the longitudinal and along-flow directions are plotted along the red dashed and dotdashed lines, respectively.

Figure 5 .
Figure 5.Comparison between the observed and estimated deformation from the modified deformation model after the global optimization process during 2011-2012, with the setting time as the actual elapsed time.The minimum root-mean-square (RMS) difference between the two deformation maps was 1.20 m.Profiles in the longitudinal and along-flow directions are plotted along the red dashed and dot-dashed lines, respectively.

Figure 5 .
Figure 5.Comparison between the observed and estimated deformation from the modified deformation model after the global optimization process during 2011-2012, with the setting time as the actual elapsed time.The minimum root-mean-square (RMS) difference between the two deformation maps was 1.20 m.Profiles in the longitudinal and along-flow directions are plotted along the red dashed and dot-dashed lines, respectively.

Figure 6 .
Figure 6.Comparison between the observed and estimated deformations from the modified deformation model after the global optimization process during 2012-2013, with the setting time as the actual elapsed time.The minimum RMS difference between the two deformation maps was 2.11 m.Profiles in the longitudinal and along-flow directions are plotted along the red dashed and dotdashed lines, respectively.

Figure 6 .
Figure 6.Comparison between the observed and estimated deformations from the modified deformation model after the global optimization process during 2012-2013, with the setting time as the actual elapsed time.The minimum RMS difference between the two deformation maps was 2.11 m.Profiles in the longitudinal and along-flow directions are plotted along the red dashed and dot-dashed lines, respectively.

Figure 7 .
Figure 7. 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.The yellow dotted line was extracted from the grounding line of the MEaSUREs dataset [13,15].

Figure 7 .
Figure 7. 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.The yellow dotted line was extracted from the grounding line of the MEaSUREs dataset [13,15].

Figure 8 .
Figure 8. 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 8 .
Figure 8. 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 9 .
Figure 9.The contribution of variations in input parameters to the output of the surface deformation model, determined by changing the variables within 10% from the estimated values.Note that the values for X and Y are plotted in the same location.

Figure 9 .
Figure 9.The contribution of variations in input parameters to the output of the surface deformation model, determined by changing the variables within 10% from the estimated values.Note that the values for X and Y are plotted in the same location.

Table 1 .
Minimum, maximum, and optimized values of each unknown variable for constrained global optimization of the viscoelastic surface deformation.
* Each deformation period was run with five unknown parameters (setting time equal to the measured elapsed time).