Retrieve Ice Velocities and Invert Spatial Rigidity of the Larsen C Ice Shelf Based on Sentinel-1 Interferometric Data

: The Larsen C Ice Shelf (LCIS) is the largest ice shelf in the Antarctica Peninsula, and its state can be considered to be an indicator of local climate change. The goal of this paper is to invert the rigidity of the LCIS based on the interferometric synthetic aperture radar (InSAR) technique using Sentinel-1 images. A targeted processing chain is ﬁrst used to obtain reliable interferometric phase measurements under the circumstance of rapid ice ﬂow. Unfortunately, only the descending data are available, which disallows the corresponding 2-D velocity ﬁeld to be directly obtained from such measurements. A new approach is thus proposed to estimate the interferometric phase-based 2-D velocity ﬁeld with the assistance of speckle tracking offsets. This approach establishes an implicit relationship between range and azimuth displacements based on speckle tracking observations. By taking advantage of such a relationship, the equivalent interferometric signals in the azimuth direction are estimated, thereby recovering the interferometric phase-based 2-D ice velocity ﬁeld of the LCIS. To further investigate the state of the LCIS, the recovered 2-D velocity ﬁeld is utilized to invert the ice rigidity. The shallow-shelf approximation (SSA) is the core of the reverse model, which is closely dependent on boundary conditions, including kinematic and dynamic conditions. The experimental results demonstrate that the spatial distribution of the rigidity varies approximately from 70 MPa · s 1/3 to 300 MPa · s 1/3 . This rigidity distribution can reproduce a similar ice ﬂow pattern to the observations.


Introduction
Ice shelves account for 75% of the Antarctic coastline, covering an area over 1.561 million km 2 [1]. They retreat/thin and collapse when the climate warms [2], which reduces the buttressing effect regulating the contribution of the continental ice sheet to the sea level. These changes of ice shelves influence the ocean circulation as they are generally accompanied by a plethora of melting water [3]. To a certain degree, ice shelves can be considered to be a potential signal of the environmental change. The Larsen C Ice Shelf (LCIS), approximately 51,000 km 2 [4], is the largest ice shelf in the Antarctica Peninsula [5] and is sensitive to atmospheric and oceanic forcing [6]. In the past few years, as warm air/ocean water drives more melting, the LCIS has continuously thinned at a rate of 0.38 m/a [7]. Especially, it calved on 12 July 2017, along with a greater than 200 km crack, forming a huge tabular iceberg. Because of discharges, the corresponding resistance decreases, which leads to an increase in inland glacier flow velocity. Such events render the LCIS minimum observation size since the birth of the remote sensing technique.
In this paper, a velocity field of the LCIS is first mapped by the differential interferometric synthetic aperture radar (D-InSAR) method [8]. It retrieves the interferometric phase to map a glacier velocity field with a higher accuracy and a finer spatial resolution [9], which is beneficial to obtaining the spatial rigidity of the LCIS more reliably. However, the interferometric observations can be only derived from one direction (descending orbit) [10], which leads to a failure in resolving the complete horizontal velocity components. Theoretically, the correlation between the displacements in the range and the azimuth directions is fixed at a given time. That is to say, as long as a velocity map is retrieved from the same SAR images, it is possible to recover the 2-D horizontal velocity vectors from the available displacements in the line-of-sight (LOS) direction. It is well known that the speckle tracking technique can obtain a 2-D velocity field with a sub-pixel measurement accuracy [11], which can be used to explore the implicit interrelation between the azimuth and the range displacements. Due to the large size and the high nonlinearity of the displacement data, the relationship is extremely complicated. In other words, it is impossible to express it as an explicit mathematical function. It has been repeatedly proofed that the artificial neural network (ANN) is capable of accurately fitting a function with arbitrary forms [12]. Noting this, a targeted ANN is constructed to depict the implicit relationship between the displacements in the azimuth and the range directions. The speckle tracking observations are used to train this network in a supervised manner. With the use of the model, an implicit function correlating the range and the azimuth displacements is generated, thereby estimating the equivalent interferometric signals in the azimuth direction based on the single direction phase observations. The proposed method has two advantages: (1) Because the ice motion is associated with several parameters (e.g., ice thickness [13], surface slope [14], and ice bed conditions [15]), it is necessary to take these factors into account when computing the function. The proposed method can extend the input to multiple variables, allowing the introduction of a series of parameters. It enables the relationship to be depicted more comprehensively, leading to an improvement in estimation reliability.
(2) The proposed method is immune to outlying speckle tracking estimates to some extent as the network is defined in a global sense.
Next, the inverse control model of the ice sheet and sea-level system model (ISSM) [16] is applied to the interferometric phase-based ice velocity field, thereby deriving the spatial rigidity of the LCIS. The model consists of the forward and the inverse parts [17]. The foundation of the forward part is actually established by a set of diagnostic equations of ice shelf flow [18], which connects the surface velocity with the depth-averaged viscosity. Especially, the viscosity is a function of effective strain rate and ice rigidity [18]. In order to avoid the interpretation trouble, the rigidity is directly substituted to calculate rather than deduced through viscosity [19]. On the other hand, the inverse part is responsible for seeking an appropriate rigidity distribution to match the observed velocity field as precisely as possible. A cost function is constructed in this part to calculate the misfit between the observed field and the modeled counterpart. The steepest-descent algorithm is introduced to minimize the cost function by adjusting the ice rigidity along with the direction of the gradient of the cost function. The strengths of the inverse control model are enumerated as follows: (1) It is unpractical to find the exact solution of the inverse parameter. By iteratively optimizing the cost function, the model can achieve the best match between the observations and the model outputs, thereby providing an optimal solution of the rigidity.
(2) The model is easily automated with a little help of the programming technique [20].

Input Data
The following datasets are employed in this paper.

1.
The ice velocity map of the LCIS is retrieved from six standard Sentinel 1A/B terrain observation with progressive scan (TOPS) Level-1 single look complex (SLC) images acquired in interferometric wide (IW) swath mode. The images were acquired on 15 May 2020 (orbit number: 21589) and 21 May 2020 (orbit number: 32660). These SLC products preserving phase information have been referenced in radar range observation coordinate using the orbit and attitude data from the satellite [21]. The heading angle of the acquired images is around −61.27 degrees, whereas the incidence angle varies approximately from 29 degrees to 42 degrees. Note that the horizontal transmit and horizontal receive (HH) polarization possesses higher signal-to-noise ra-tio and supplies better amplitude characteristics compared to the horizontal transmit and vertical receive (HV) polarization, implying that the former is more suitable for mapping the ice motion [22]. Therefore, the HH channel is chosen in this paper.

2.
The MEaSUREs BedMachine Antarctica (Version 2) [23] dataset provides the necessary ice thickness ( Figure 1a) and the surface topography ( Figure 1b) for the modeling manipulation. The bed parameter is derived from the difference between the surface elevation and the ice thickness. The ice thickness and the surface elevation are presented in ice equivalent, which avoids the trouble of correcting for the firn density [24]. These data have a spatial resolution of 500 m and are geolocated in a polar stereographic projection with the standard longitude of 0 • E and latitude of 71 • S.

InSAR Processing
Due to the rapid flow of the glacier, the traditional procedure of the D-InSAR technique cannot be directly utilized to generate an interferogram with respect to the glacier velocity [25]. A targeted InSAR processing chain is used to deal with the input data (see Figure 2), which is roughly described as follows.
(1) A standard Sentinel-1 SAR image in IW mode generally consists of three subswaths, and each sub-swath usually comprises nine bursts. Fractional areas of neighboring bursts overlap in azimuth. Similarly, neighboring sub-swaths cover a certain common area in range. Consequently, burst alignment operations are first applied to each sub-swath of the TOPS acquisitions according to the time parameters. Following this, the corresponding sub-swaths are merged to generate the consecutive master and slave magnitude images.
(2) With the use of the external digital elevation model (DEM) and the precise orbit information generated by Copernicus precise orbit determination (POD) service, the relative offset of each measurement point is estimated using the geometric coregistration. As a result, an initial rough look-up table is generated.
(3) According to the look-up table, the slave magnitude image is resampled with respect to the master magnitude image.
(4) The offsets between the master and the resampled slave images are obtained by the speckle tracking technique with a window size of 32 × 128, which utilizes the classic cross-correlation function to determine the range and the azimuth offsets within a sub-pixel measurement accuracy. The speckle tracking result also plays another crucial role in this paper. It is used to supervise the training process of the artificial neural network depicting the relationship between the range and the azimuth displacements (see Section 3.1).
(5) The look-up table is updated based on the speckle tracking estimates, and the slave image is resampled again according to it. (6) As the accuracy of the speckle tracking based look-up table cannot meet the requirement of TOPS image registration, inevitable phase jumps will be presented at burst/swath edges. Following the idea proposed in [26], phase jumps are estimated and compensated.
(7) Eventually, the interferogram is generated based on the result in (6). Note that the interferometric phase components induced by topography are estimated based on the external DEM and removed during the process of interferogram generation.

Merging Burst
Merging Burst

Merging Swath Merging Swath
Precise Orbit External DEM

Geometric Coregistration
Resampling SLC1 SCL2 Look-up Table   Speckle Tracking Update Look-up Table   Resampling Phase Jump Mitigation

Master Image
Slave Image Interferogram Figure 2. The interferometric synthetic aperture radar (InSAR) processing chain used in this paper.

Ice Velocity Field Recovering
Based on the speckle tracking technique, the 2-D ice surface velocity field of the LCIS can be produced with the use of the offsets in the range and the azimuth directions. Theoretically, there exists an interrelation between the displacements in the two directions. However, it cannot be expressed as a determinate function (e.g., a sine function or an exponential function). This study describes this implicit function as follows: where v azi and v ran refer to the azimuth and the range displacements, respectively. f (·) is the function delineating the correlation between v azi and v ran . The vector x is a series of physical parameters associated with v azi and v ran , including spatial coordinates, surface slopes, ice thickness, etc. Frankly, it is hard to find a function f (·) by the traditional fitting method, as the large dimensions of the input-output data bring about severe nonlinearity and unacceptable computational time. In recent years, the artificial neural network widely attracts the focus of researchers [27,28], because it has made several breakthroughs in science and society [29,30]. It is an information processing system which imitates the human's brain. It plays an important role in highly complex, nonlinear, and parallel computing. The architecture of the network used in this study is illustrated in Figure 3. In this network, the range displacements obtained from the speckle tracking technique, the coordinates of the points in the ice displacement field, the surface slopes, and the ice thickness are regarded as the inputs. As the input layer is desired to have as many independent data dimensions as possible, surface slopes are decomposed into the partial derivatives of the surface elevations. The azimuth displacements gained from the speckle tracking technique are assigned as the labels associated with the output layer. The depth of the hidden layers is set to 4, and the size of each hidden layer is uniformly set to 100. Tanh is set as the activation function of each hidden layer. After the construction of the network, we commence on the hyperparameters related to the training (including the learning rate, the batch number) and the optimization strategies (including the optimization algorithm and the regulation). The learning rate and the batch number are set to 0.001 and 32, respectively. Adam is chosen to be the optimization algorithm. It combines the advantages of the AdaGrad and RMSProp, which is helpful for handling sparse gradients on noisy problems. In order to avoid model overfitting, L2 regulation is carried out.

Inversion for Rigidity
As ice can be characterized as fluid, the Stokes equations is widely used to depict the glacier ice flow [16]: Equations (2) and (3) represent the stress balance (neglecting the acceleration and the Coriolis forces) and the incompressibility, respectively, where ρ is the ice density, u is the ice velocity vector, T is the Cauchy stress tensor, p = − 1 3 ∑ i T ii is the pressure, g is the gravity acceleration, and ∇· denotes the divergence operator. The glacier ice is considered as viscous density-preserving fluid, and its constitutive relation can be expressed as where T D refers to the stress deviator. It is split from the Cauchy stress tensor: where I is the identity matrix in mathematics. D is the strain rate tensor, which can be further expressed as Note that the trace of the D is equal to the divergence of u, thus the value of tr(D) is also zero. µ is the ice effective viscosity. Due to the glacier ice belonging to the non-Newtonian fluid, µ follows the Glen's flow law [31]: where B is the ice rigidity (i.e., the objective of this inversion), n the Glen's flow law exponent (typically chosen as 3). D e is the effective strain rate and defined as Let x and y be horizontal coordinates, z be the vertical coordinate (positive upward). With the use of Equations (4)-(6), Equations (2) and (3) can be rewritten in the form of the velocity components in the Cartesian coordinates [16]: ∂u ∂x + ∂v ∂y where u, v, and w are the components of the velocity vector u in x, y, and z directions, respectively. These equations are the so-called full-Stokes model, which contains four unknowns (u, v, w, and p). It is computationally time-consuming to find the solutions. In practice, ice shelves are relatively flat and thin. In addition, the viscosity of sea water is negligible compared to that of glacier ice. Therefore, it is reasonable to assume that the horizontal velocities and strain rate are independent of z. Based on the assumption, the full-Stokes model can be rearranged. The shallow-shelf approximation (SSA) is thus obtained [18]: Additionally, Equation (7) where z s and H are the elevation of an ice shelf and the local ice thickness, respectively.
Obviously, the SSA model has only two unknowns (u and v), the computational complexity is therefore reduced. It must be noted that the viscosity in the Equations (13)- (15) refers to the depth-averaged effective viscosity due to the temperature dependence of ice rigidity.
To solve the SSA model that consists of two partial differential equations, the appropriate boundary conditions are required. In this paper, both the kinematic and the dynamic conditions are introduced. The kinematic conditions specify the velocities (usually the observed surface velocities) at locations neighboring the grounding line of an ice shelf. On the other hand, the dynamic conditions apply the stress constraints to the ice front, which is the depth-integrated balance: where n is the outward-pointing normal vector of the ice front, ρ w is the density of sea water.
Inversely, the SSA model can be also considered as a direct approach to solve the depth-averaged effective viscosity when the velocities, the local ice thickness, and the ice surface elevation are known. However, it may not be applicable in practice mainly due to the errors in the observed velocities [20]. To deal with the problem, a control method is introduced to estimate the distribution of viscosity [20]. The method treats the SSA equations as the "forward model" to compute the modeled velocity (u m , v m ) beginning with an initial viscosity value. Moreover, it adjusts the viscosity iteratively by gradient descent algorithm to achieve the best match between the modeled velocities (u m , v m ) and the observed velocities (u d , v d ) under the constraints given by Equations (13) and (14). The misfit between the model outputs and the observations is evaluated by a least-squares cost function with Lagrange-multipliers: where Ω is the study domain and (α, β) is the Lagrange-multiplier vector. In mathematics, obtaining an extremum of the cost function necessitates that all the variations of the function with respect to the free variables are equal to zero. Note that the variations of α and β are naturally equal to zero due to (u m , v m ) being obtained from the "forward model". The variation of J with respect to µ is written as [ To invert the ice rigidity directly, Equation (15) is used to deduce the variation of J with respect to B. Due to the noise presented in the observed velocities, a Tikhonov regularization term is introduced to the cost function to stabilize the inversion [16].
The partial differential equations are discretized by the continuous Galerkin finite element method [33], and implemented with the use of triangular Lagrange P1 elements in the ISSM, which helps us improve the computational efficiency by constructing meshes with different resolution on the whole study region. In addition, it is beneficial to preserving small-scale features in the area where the ice dynamics are of great concern.

Results
Note that the results obtained from both the speckle tracking and the D-InSAR techniques are projected to the polar stereographic coordinate from the radar coordinate. As the observations in the non-glacier regions are not associated with glacier activities, they may impose negative effect on the ANN performance and the inversion of the ice rigidity. Therefore, the displacement estimates outside the LCIS were masked out. Figure 4a,b represents the speckle tracking displacement components of the LCIS in the range and the azimuth directions, respectively. It can be easily observed that there is a primary signal leap extending from the top left to the bottom right. Moreover, a few small distinct inconsecutive signals are presented in Figure 4b. This phenomenon should stem from the influence of the inaccurate radar timing information. Figure 4c demonstrates the LOS displacement map transformed from the interferometric phase. Essentially, Figure 4a reflects an identical motion information with Figure 4c. They both exhibit that the range motion varies from −6 m to −2 m. The negative values indicate that the ice moves along the opposite direction of the LOS. The minimum motion is observed around the promontories, whereas the maximum motion is located in the neighboring region of the ice front. Compared to Figure 4a, Figure 4c is much smoother (without evident jumps), primarily due to the higher spatial resolution and measurement accuracy of the interferometric phase based estimates. With the use of the ANN, Figure 4d maintains an almost same distribution pattern of the displacement as that in Figure 4b. Moreover, it gets rid of abrupt signals, which mainly because of two reasons. First, the introduction of the surface slope information and the ice thickness to the ANN forces the model to match the mechanism of ice flow, which mitigates the influence of outlying estimates. Second, the ANN is trained with regulations, the high-frequency noise is therefore eliminated. In this study, the results in Figure 4c,d are used as the observations for the inverse model. Figure 4b,d reveals that the maximum azimuth displacement is located at the northern part of the shelf. Further, both the positive and the negative azimuth displacements are observed, which implies the significant change of the ice flow direction. Figure 4e,f shows the displacement magnitudes of Figure 4a,b and Figure 4c,d, respectively, which demonstrate that the ice displacement speeds up from the inlet areas of the shelf to the ice front and reaches a maximum value (nearly 14 m) at the middle of the ice front.
Due to the phenomenon of interferometric decorrelation, there are several no data points in the ice displacement fields represented by Figure 4c,d, which may exert negative influences on the model inversion process. Noting this, these no data points were recovered using the method proposed in [34]. Next, the complete ice displacement fields were decomposed and transformed into two mean annual velocity components in the Cartesian coordinate (Figure 5a,b) to satisfy the requirement of the model inversion. Evidently, the x-velocity component approaches zero in the northernmost and the southernmost regions. On the other hand, the y-velocity component is nearly zero in the northernmost areas. The maximum values of the two velocity components are both located in the middle of the ice front. The study domain is adaptively discretized into 260,190 finite elements in the inverse model. Figure 5c,d represents the magnitudes of observed mean annual velocities and modeled counterparts, respectively. The black arrows indicate the ice velocity vectors. In general, the ice moves along a uniform direction from the inland to the sea, and it speeds up in the same direction. Figure 5e demonstrates the misfit between the observations and the model outputs. It indicates that the disagreement between them is small. Most discrepancies of the two velocity fields are distributed within 12 m/a.     Figure 6 exhibits the depth-averaged spatial distribution of the ice rigidity inverted from the ice velocity fields (Figure 5a,b). The rigidity varies in a huge span, approximately from 70 MPa·s 1/3 to 300 MPa·s 1/3 , and its distribution pattern is not random. For the convenience of description, five typical areas are circled in purple dash line, and they are marked as A-E, respectively. The ice in area A is relatively stiff (with a rigidity of 250 MPa·s 1/3 ), which is corresponding to an isothermal ice at approximately −28 • C (see Table 1). As the inverted rigidity is depth-averaged, the estimated temperature according to Table 1 is logically between the base and the surface temperatures. The soft ice is presented in areas B and E, which generally has a rigidity within 100 MPa·s 1/3 . Correspondingly, the temperature in these areas is no less than −7 • C. Note that the soft ice is observed along the three distinct rifts in area E. Area C is of great interest, as strong bipolar singularities of the rigidity are presented in a relatively small area.

Analysis of Ice Rigidity Distribution
Rigidity refers to the property of an object to resist deformation when forces are exerted on it. On the basis of Glen's law [31], the proportional factor between deviatoric stress and strain rate is the so-called glacier ice rigidity. Its value is dependent on temperature, water content, grain size, impurities, and preferred-orientation fabric [35]. Frankly, it is hard to obtain the variables like impurities for most researchers. However, temperature is relatively within reach, as glaciologists can obtain temperature data from laboratory experiments or public meteorology archives. Consequently, the relationship between temperature and rigidity is best understood. Due to the lack of adequate theoretical formulations about rigidity, researchers calculate it only through the temperature distribution over interested areas in most ice shelf modeling applications [32]. Obviously, the value of rigidity obtained by this method is not accurate.
In this paper, the rigidity distribution of the LCIS is inverted based on the 2-D velocity field, which provides us with an opportunity to investigate the state of the LCIS further. Several rifts and crevasses can be observed in the area B, D, and E from the base map of the topography, where the ice rigidity is low. To some extent, the effectiveness of the inversion model used in this paper is therefore confirmed. Figure 7 is the driving stress deduced from the modeling results. It can be observed that higher driving stress usually accompanies with the rifts and crevasses, which brings about higher strain rates and lower ice rigidity. Note that no information with respect to rifts and crevasses is involved in the model inversion process, because the field measurements have related the fractures to higher strain rates and larger velocities. The ice in the regions closed to the grounding line is relatively stiff, and its rigidity is between 150 MPa·s 1/3 and 250 MPa·s 1/3 . This is possibly due to three factors. First, the basal melting is dramatic around the grounding line, which causes a thinning of the ice shelf and a fall of the temperature, the rigidity thus arises. Second, the advection of cold ice from tributary glaciers reduces the temperature and enlarges the ice rigidity. Third, abutting the grounding line, the sea water beneath the ice shelf is comparatively shallow. Consequently, vigorous tidal mixing can easily transfer heat to the ice shelf, which triggers a great deal of melting, thereby increasing the rigidity. In area A, where the ice is relatively stiff, the accretion of marine ice might be extremely plentiful. As a result, it may isolate the ice shelf from warm sea water to a large extent, which might cool the temperature-depth profile. The existence of the strong bipolar singularities in area C is a little difficult to interpret. As there is no evident stagnant ice, the ice over this area should be very stiff. Such an abnormity is never mentioned in previous research works. It is suspected that this may result in a systematic bias, due to the non-local characteristic of the "forward model" [36].

Importance of Control Method
Traditionally, the ice rigidity parameter B is calculated based on the temperature. Unfortunately, it is hard to obtain the temperature data with a high spatial resolution. In most cases, B is usually inferred to be a uniform value over an extensive coverage, which is not capable of depicting the real state of the ice rheology and leads to inaccurate results in many glacier modeling applications. Owing to the use of the control method, a rigidity field with significant spatial variability can be generated. The velocities (Figure 5d) computed based on the inversion result reproduce the pattern of the observations accurately. To examine the effectiveness of the control method for inferring the ice rigidity parameter, the mean annual velocities were directly calculated by the "forward model" with a uniform B (Figure 8a). The value of B is set to 175 MPa·s 1/3 according to the recommended temperature (−19 • C) in the LCIS [37]. Figure 8b demonstrates that the misfit between the modeled velocities and the observations. Obviously, the disagreements are approximately 35 m/a in most areas. These areas are mainly covered by soft ice.

Conclusions
Because only descending data are available over the LCIS, it is directly impossible to retrieve the corresponding 2-D velocity field using the D-InSAR technique. In this paper, a new method is proposed to recover the azimuth displacement. This basic idea of this method is stimulated by the theory of ANN. A targeted network is constructed to fit the implicit relationship between range and azimuth displacements based on speckle tracking results. It introduces physical parameters related to the ice motion as the inputs, which constrains the inputs and outputs to the mechanism of glacier flow to some extent. It is beneficial in reducing the influence of the outlying estimates and noise. Next, the flow mechanism of ice shelves is discussed. The corresponding inversion implementation in the ISSM is used to infer the ice rigidity. It seeks the best rigidity distribution by minimizing the misfit between the interferometric phase-based 2-D velocities and the modeled counterparts. Compared to the uniform ice rigidity calculated from a fixed temperature, the inverted variable rigidity can reproduce the observed flow patterns more accurately.