remote Assessment of Recent Flow, and Calving Rate of the Perito Moreno Glacier Using LANDSAT and SENTINEL2 Images

: We mapped ﬂow velocity and calving rates of the iconic Perito Moreno Glacier (PMG), belonging to the Southern Patagonian Iceﬁeld (SPI) in the Argentinian Patagonia. We tracked PMG from 2001 to 2017, focusing mostly upon the latest images from 2016–2017. PMG delivers about ca. 10 6 m 3 day − 1 of ice in the Lago Argentino, and its front periodically reaches the Peninsula Magallanes. Therein, the PMG causes an ice-dam, clogging Brazo Rico channel, and lifting water level by about 10 m, until ice-dam failure, normally occurring in March. Here, we used 36 pairs of satellite images with a resolution of 10 m (SENTINEL2, visible, 9 pairs of images) and 15 m (LANDSAT imagery, panchromatic, 27 pairs of images) to calculate surface velocity ( V S ). We used Orientation Correlation technique, implemented via the ImGRAFT ® TemplateMatch tool. Calving rates were then calculated with two methods, namely, (i) M1 , by ice ﬂow through the glacier front, and (ii) M2 , by ice ﬂow at 7.5 km upstream of the front minus ablation losses. Surface velocity ranged from about 4 m day − 1 in the accumulation area to about 2 m day − 1 in the calving front, but it is variable seasonally with maxima in the summer (December–January–February). Calving rate ( CR M ) ranges from 7.72 × 10 5 ± 32% to 8.76 × 10 5 ± 31% m 3 day − 1 , in line with recent studies, also with maxima in the summer. We found slightly lower ﬂow velocity and calving rates than previously published values, but our estimates cover a different period, and a generally large uncertainty in ﬂow assessment suggests a recent overall stability of the glacier.


Introduction
The Perito Moreno Glacier (hereon PMG) is located in the Campo de Hielo Patagonico Sur (Southern Patagonian Ice field, hereon SPI), the largest temperate ice masses in the Southern Hemisphere and the largest continental (not polar) reserve of fresh water therein.
Recent climate change in the central Andes caused measurable reduction of snow cover [1]. Moreover, the ice cover was largely reduced recently, and this is expected to accelerate in the future [2,3]. The question arises whether South Patagonia will also undergo such phenomena.
Patagonia is a main landmass at its latitude, so SPI is of great significance for studies of climate history and ice-climate interactions [4][5][6][7]. Several outlet glaciers flow radially from SPI, mostly calving into fjords on the western side, and into lakes on the eastern side [8,9]. During the last 50 years, out of the 22 major calving glaciers in Patagonia, nine have been fluctuating within ±1 km, and 12 have been retreating considerably in a range from −1 to −13 km [10,11]. PMG instead experienced lower fluctuation during the last half-century, albeit the glacier snout oscillated frequently, causing damming of

Input Data
To investigate calving rate of PMG from 2001 to 2017, we used (i) open-source satellite imagery to derive motion vector fields of the glacier, (ii) the perimeter of the glacier, and (iii) ablation measurements from 1995 to 2003 [27].
We used images from the LANDSAT and SENTINEL2 collections. The images covering the PMG area are positioned at 50 • 17 11 S, 73 • 23 12 W (path: 231, row: 95). Both collections are sensible to cloud cover, due to the optical nature of the sensors. Here, we used images in pairs, with a maximum time window of 100 days. In Table 1, we report the pairs used.
We used LANDSAT4-5 images (from 1999 to 2017, with 30 m of spatial resolution, which was found too low to be used here after a preliminary screening) and LANDSAT7 (from 1999 to 2003) and LANDSAT8 (from 2003 to 2017) images with a spatial resolution of 15 m (panchromatic band).
Eventually, the processed cloud-free images from 1999 to 2017 were 62 for LANDSAT7 and 24 for LANDSAT8. Some images were discarded due to a large acquisition interval or the low quality of the feature tracking.
In addition, we considered nine images in natural color from SENTINEL2 (bands 2, 3 and 4 corresponding to BGR, at 10 m of spatial resolution) to map the glacier's front variation. Most of our processed images eventually covered the years 2014-2017, providing a most recent estimation of GPM dynamics.
The outline of the glacier was defined by the GLIMS glacier database relative to the year 2012, published by the National Snow and Ice Data Center (http://nsidc.org/, (accessed on 21 December 2021). To investigate the topography of the glacier, the ASTER Global Digital Elevation Model (Version 2) was used. This DEM has a spatial resolution of 30 m, and it is delivered by NASA and freely available at https://earthexplorer.usgs.gov/ (accessed on 21 December 2021). We used topographic features, such as slopes and aspect, to qualitatively validate the estimated fields of velocity on the glacier. Table 1. Pairs of images used as input data. ∆t is the time window between the two images (days), and SNR AVE is the mean signal to noise ratio. Albeit meteorological data can be very useful in measuring and modeling glacier dynamics (i.e., snow and ice ablation and accumulation whenever occurring, e.g., [3]), here no such data were available. The nearest weather station is at El Calafate (50 • 20 24 S, 72 • 16 13 W, 199 m a.s.l.), about 60 km from the glacier front and thereby little usable for the purpose of linking climate variation to PMG dynamics. The unique available meteorological dataset at the glacier surface that we know of covers a very short time period, namely, Stuefer [29], which measured monthly temperatures at the glacier front from 1996 to 1998. While for the lower tongue it is possible to derive a (linear) relation that describes the temperatures on the glacier [30], this is not possible for the accumulation area, which is too high and far from any weather station. Therefore, it is not possible to distribute air temperature along the whole glacier, even for short periods (1996)(1997)(1998). Finally, grid data are available, such as VASClimO used by Stuefer et al. [27], but they are not as reliable as field measurements, especially in areas with such a complex topography. Concerning precipitation, Stuefer et al. [27] used monthly data from the closest grid point of VASClimO as a proxy.
Ablation measurements could be used in our study. Stuefer et al. [27] measured annual ablation in profiles A, B, D, and L ( Figure 1) from 1995 to 2003, including flow velocity at different stakes. From these measurements, a significant difference in flow velocity occurs between the central zone of the glacier and the marginal and crevassed zone, where ablation is higher, due to a larger surface exposed and particularly turbulent air fluxes. Thus, we divided the area in two zones, namely, zone 1 (central) and zone 2 (marginal and crevassed), and ablation was computed as the weighted (against area) mean of the ablation measurements in the two zones.

Surface Velocity
The surface velocity was computed here by comparing satellite images taken in different days. This comparison was made by means of ImGRAFT-Image Georectification And Feature Tracking (http://imgraft.glaciology.net/home, accessed on 21 December 2021), an open source toolbox, developed by A. Grinsted and A. Messerli in MATLAB code [24]. TemplateMatch is a function of ImGRAFT, used to investigate the shift between two images, based upon Orientation-Correlation techniques. The two images are defined as image A (the first taken) and image B (most recent one). This function allows to create a velocity field on the entire area of the glacier. This is different from classic feature tracking, which computes velocity values only on reference points. Calibration and validation of the method are specific for every glacier and in general for every spatial body. Thus, these steps are not yet automatic, and supervision has an important role.
Even if the image pixel size is in the order of magnitude of meters, the displacements computed by ImGRAFT TemplateMatch may be in the order of centimeters.
Highest mean correlation is normally used as a goodness of fit parameter. However, only validation against field measurements or supervision on image classification can define the accuracy of the estimated displacement and the magnitude of noise.
In addition to supervised classification, to avoid false matches or outliers (e.g., cloud cover), four indexes are computed in the code: (i) mean correlation index (C AVE , mean of the largest correlation matrix for each match), (ii) mean Signal to Noise Ratio (SNR AVE ), (iii) the percentage of valid cells (C % , cells with SNR < 2 are automatically discarded), and (iv) the SNR corrected with the number of valid cells (SNR COR , here we took SNR COR ≥ 10, see e.g., [31]). The SNR AVE is the ratio between signal and noise typical of the signal analyses, defined as: where C AVE is the mean correlation coefficient between template A and B (−1 to +1, cells with C ≤ 0.6 are automatically discarded). C NOISE is the module of the average of the correlation coefficients, across the whole search region, and it is an estimate of the noise level. Here, we applied a filter for velocity values based on thresholds, namely: (i) V S,max ≤ 5 m day −1 on the entire glaciers (whenever necessary, this was raised to 7 m day −1 ); (ii) V S,max ≤ 4 m day −1 and no West ward components in the ablation area; (iii) V S,max ≤ 2 m day −1 and flow direction only towards East in the front area, covering approximately the last 4 km.
Such thresholds, based upon observation of the velocity fields over the many images, were used to remove coarse noise and errors and (too) large movements detected, as sometimes observed over large time windows.
Across glacier tracks are used here. Some tracks overlap those of Stuefer et al. [27], and others cover the glacier area and front, with a denser scheme. The field measurements in [27] are used to benchmark our results.

Calving Rate
The calving rate (CR) depends on the glacier dynamic, meteorological conditions, and the features of the sea or lake where the glacier ends. First, one needs to know whether the PMG tongue is floating or resting upon ground. Critical buoyant thickness is given by [32].
where H M is ice thickness at the front margin, H w is water depth, ρ w and ρ i are the density of water and ice, respectively. Because H M is 50-80 m higher than H w (about 150 m, [27]), PMG tongue should rest on ground everywhere. To quantify the calving rate, we followed the approach reported in [27], based on the velocity at the glacier front. From surface velocity one can derive vertically averaged glacier velocity as: where V is the mean velocity in m day −1 , and Vs is the surface velocity representative for a certain point or cross-section area of the glacier. This correction factor is suggested by Cuffey and Paterson [32], as a widely used approximation in glaciology, and we used it here accordingly. The glacier flow rate is computed multiplying velocity V by the corresponding crosssection area (area-velocity method). If a cross-section is split into several sub-sections, and each subsection displays a velocity measurement or estimate, the sum of flow rates of each sub-section gives the total flow rate. Flow rates are computed at track B and at the front of the glacier. This allows two different approaches for the calving rate estimation.
The first method, M1 is based upon computation of glacier flow rate at the front. Given that all calving ice has to cross the front cross-section, and it has to collapse into the Lago Argentino, it is possible to calculate the net advance of the glacier's front as: where CR M1 is the volumetric total calving rate [m 3 day −1 ] using this method, Q F [m 3 day −1 ] is the flow rate (speed) at the front of the glacier, and G F [m 3 day −1 ] is the net advance of the glacier front position. Here, we embed into CR M1 horizontal ice ablation or melting rate at the front and net calving rate, consistently with current literature (e.g., [22]) given that separation of such components is complex. To approximate glacier's cross sections for flow rate assessment, we used the transects as reported by Stuefer et al. [27] (see Figure 4 in [27]). SENTINEL2 images in natural color (bands 4, 3, and 2) with 10 m of spatial resolution were used to map the glacier front variations.
The second method M2 is based on the glacier flow rate computed at the B profile. The rationale is similar to approach M1, but it is necessary to introduce some assumptions. The ice in Profile B needs on average about 4000 days to reach the front. Therefore, such analysis is only possible thanks to large stability of the PMG, as confirmed in many studies and by the analysis of the ice speed variations during seasons and years made in this study. Second, from Profile B to front, it is likely that mass loss due to ablation is not negligible. Ablation measurements are here available from Profiles A and B, quite uniform, but not from Stake D (at the front) and from Profile L (above Profile B) [27]. Moreover, at Stake D, ablation is likely to be very high due to the influence of the lake and the lower altitude there. In order to compute the ablation on the whole area, two subsections are used [29]: central zone and marginal and crevassed zone. Moreover, altitude and climate data are considered when computing the ablation data as function of the altitude, to not omit climatic and glaciological factors. Nonetheless, neither with altitude nor temperature is it possible to explain the increase of ablation from Profiles A, B, and L, to the Stake D. Accordingly, ablation is computed without considering the Stake D and producing a systematic underestimation of mass loss, and thus an underestimation of the calving rate. Inserting ablation loss in the approach, the calving rate can be computed as: where Q B [m 3 day −1 ] is the glacier flow rate at Profile B, and Q A [m 3 day −1 ] is the mass loss by ablation in the whole area between the Profile B and the glacier front. We tentatively implemented both such methods, to verify whether (scarce knowledge of) ablation rate and mass balance may affect calving rate assessment.

Model Calibration
Before actual model tuning, we pursued a preliminary assessment of the best crosscorrelation method in TemplateMatch for the PMG area, among the four available methods (i.e., NCC, OC, PC, and CCF) as reported above.
The algorithm takes in input two subsequent (in time) images covering the same area (A-B), and outputs the correlation of any cell in image A with all the cells in image B, belonging to a certain research region. The largest correlation of each cell in image A with a cell in image B provides a guess of the new position of the feature in image A within image B (i.e., of its displacement in time), and according to the difference in time, an estimate of the velocity. All methods use Fast Fourier Transformation, to allow detection of movement smaller than the image pixel resolution.
Here, a preliminary adoption of these four methods was pursued to assess flow velocity of the PMG. This is fully reported in Chirico [33], while here we report the main findings.
Preliminary application of the PC [34] and CCF [35] methods did not produce significant results in the PMG area ( Figure 2). The NCC [36] method provides acceptable results, but some edge noise is still present, and the method is not optimal, given that it cannot manage well image rotations and distortions ( Figure 2). The OC method [37] provides good results, with a proper signal to noise ration SNR ( Figure 2). This method works on the intensity gradient of the images and so the method is fast, statistically robust, exhaustive, and invariant to image brightness. In addition, working with correlation, some weak spots of optimization are avoided. For instance, the OC method does not need an initial state, and it does not need a priori knowledge (such as, e.g., NCC), and thanks to the exhaustive research, local minima do not lead to sub-optimal solutions. Accordingly, we decided to use OC methods for PMG tracking here.
Further, TemplateMatch behavior was tested against different spectral bands, given that the algorithm uses single bands (i.e., not RGB visible bands), which can have different resolution. We initially used visible bands (i.e., 2, 3, 4 of SENTINEL2, and 1, 2, 3 of LANDSAT). Subsequently, we tested panchromatic band of LANDSAT (band 8) which has high resolution (15 m). SENTINEL2 images have RGB bands with the same resolution as panchromatic (10 m), and TemplateMatch provided here similar results for all bands.
The time step (∆t) between two images A-B also affects the results. Given the lack of debris cover on the glacier usable for reference and the relatively high flow velocity (≈ 1 m day −1 on average), a too large value of ∆t would not be suitable, because the characteristic of the images would change too much. Crevasses can be used as target features, but crevassed zones on the glacier tend to remain the same in time, so not actual movement may be found. The largest value of ∆t is 80 days, as reported in Table 1. Moreover, given image availability and cloud cover in the period of study (2002-2017 with LANDSAT), the lowest ∆t here is 13 days (see Table 1).
A first parameter to be chosen for model tuning is the number of templates (i.e., spacing between template centers) to set into image A and then follow in image B. This parameter does not influence the results largely, but the higher the number of templates, the (exponentially) longer the computational time. A regular grid of template centers is used in this study, with density fixed, so as to investigate different areas and phenomena (e.g., a dense grid of 5 by 5 pixels is used for the analysis made on the sole front, more interesting to assess calving, and a sparser grid of 10 by 10 pixels is used to investigate the largest ablation area).
The template size (TS) is defined by width and height (that are equal for a square template as used here, with side length TS) over the image A (older image) where FFT is pursued. The size of the template is an important parameter. A large template is more efficient, so it is less likely that false matches are attained, i.e., finding the same (similar) traits in a wrong region of the image B. On the contrary, a small template results into low computational time. A trade-off between these two objectives provides fast and efficient calculation. The search region over the image B (more recent) can also be either square or rectangular, but here, we used a square shape with side length RS (search region size), which is more regular.
A better knowledge of the glacier motion allows to define the smallest search region as possible. Here, after a preliminary analysis, we found that a square template with TS = 70 pixels in size provides the best performance, without very large computational time. The default size of RS, also used here after a preliminary assessment, is RS = TS + 41.
The last tuning parameters dx and dy define the initial hypothesis of the 2-D directional displacement of the template TS from the initial point. An accurate initial hypothesis can improve the performances, allowing to use a smaller Search Region, already centered on the likely final destination of the template. Here, the dx and dy values were adjusted iteratively based upon preliminary assessment of flow velocity for each pair of images, and depending upon flow velocity and time lag ∆t.

Surface Velocity Validation
After tuning of the function TemplateMatch of ImGRAFT as reported above, we eventually derived the surface velocity for the whole glacier (see Figure 3, giving Vs during the austral fall/winter of 2016). The estimated values are in general agreement with topography. The Southernmost part of the glacier shows faster velocity (about 3-5 m day −1 ), likely due to narrowing of the flow section, leading to an acceleration of the glacier flow. Along the ablation tongue, the movement and direction are more regular because there are no obstacles or sudden changes in section area or in slope. Close to terminus flow is split in two parts. In front of the Peninsula Magallanes, the lake bottom is higher, and from there, the glacier follows two different directions (i.e., towards Lago Argentino and Brazo Sur), where the speed is about 2 m day −1 higher than in the central area.
The accumulation area is poorly studied because of harsh conditions and perennial snow cover. Snow cover makes it difficult to reckon reference points when using feature tracking methodologies (such as those used by TemplateMatch). In addition, this area is often covered by clouds, and accordingly, pairs of images with clear-sky conditions at a short distance over time are rare to find. For these reasons, we validated our results only along the ablation tongue.
To better assess the velocity fields so obtained, we defined some cross-sections of the glaciers, where velocity profiles were studied, and benchmarked against studies in the literature. We used 10 sections (transects), of which 5 (TA Transects, Figure 3b within the broad ablation area and 5 along the PMG front (last 4 km, TF transects, Figure 3c. In Figure 4, we report for TF transects the estimated surface flow velocity (average of all images) and a comparison between the two satellites LANDSAT and SENTINEL2 and their image resolution, from 30 to 10 m of pixel size. We performed a comparison using all images since there is not a common period between LANDSAT4-5 and SENTINEL2A. Therein, it is shown that images with 30 m resolution (LANDSAT4-5; see Table 1) and at 20 m resolution (SENTINEL2) provide systematic underestimation of flow velocity. Accordingly, we decided to discard the images from LANDSAT4-5 here.    Table 1). Moreover, we report the surface profile of the glacier along the transects, as from ASTER Global Digital Elevation Model (Version 2), and an interpolated bedrock profile, according to Stuefer et al. [27] and Stuefer [29].
In Figure 5, we report a more accurate comparison between SENTINEL2 and LAND-SAT8, by providing mean flow velocity (here, in transect TA3 in the ablation zones, with similar results in all transects, also in front zone; see [33]) during three broadly overlapping periods in 2016 (winter 2016, ID3, ID19; end of summer 2016, ID15, ID24; and start to mid summer 2016, ID16, ID5; see Table 1). Therein, the two summer periods differ slightly in the shear margins, possibly given the temporal distance of about 1 month.
Validation was pursued against estimates of flow velocity of the PMG as reported from [27]. Namely, we compared flow velocity at transect B therein, in Figure 1, by creating a new section TB62, slightly downstream of TA2 here, and subsequently assessing flow velocity thereby. Figure 6 shows the results so obtained. We report therein the different flow velocities estimated at Transect B in [27], namely, from (i) tracking of ablation stakes (1995-2003), (ii) assessment of velocity field from SIR-C/SAR-X (October 7, 9, and 10, 1994), (iii) SENTINEL2 (RGB, 10 m pixel size), and (iv) LANDSAT7-8 (panchromatic, 15 m pixel size), in this study. Despite the different periods of assessment, the velocity patterns along the glacier seem to be quite similar.  Table 1). Also; we report the surface profile of the glacier along the transect, as from ASTER Global Digital Elevation Model (Version 2) and an interpolated bed profile, according to Stuefer et al. [27] and Stuefer [29]. Notice that the estimated velocities from TemplateMatch were lumped as a mean value from those of the adopted satellite, while the velocity field extracted from the SIR-C/X-SAR refers to one short period only (7, 9, and 10 October 1994). Tracking of ablation stakes as shown in Figure 6 covered a period of 8 years or so. We estimated flow velocities within the areas where Stuefer et al. [27] mapped movement of their ablation stakes, and accordingly, we could compare our estimates against their stakes' velocity data. In Table 2, we report goodness of fit indicators of our estimated flow velocity of the PMG. Table 2. Goodness of fit of the estimated flow velocity of the PMG against ground-based assessment using tracking of ablation stakes by Stuefer et al. [27].

Satellite
Bias

Calving Rates and Validation
In Figure 7, we report our estimated calving rate (CR), at the front (M1) and at transect B (M2), against estimates reported in [27], together with the corresponding flow rate at front in m 3 day −1 , indicating relative advancement (+) or retreat (−). Methods M1 and M2 generally provide good agreement, with the exception of the paired images ID11 in 2017 (2 June-12 July 2017) and mostly CR M1 < CR M2 (unless for ID7, 20 October-2 November 2017). Both methods deliver lower CR than estimated before, also possibly given the approximated positioning of the transect TB62, ideally overlapping their transect B (in Figure 1b). On average, we estimated here for M1, CR M1 = 7.72 × 10 5 ± 32% m 3 day −1 , and for M2, CR M2 = 8.76 × 10 5 ± 31% m 3 day −1 , against 1.08 × 10 6 m 3 day −1 reported in [27], i.e., −29/−19% here. However, one has to take into account that (i) CR in [27] refers to one single estimate, based upon an ice velocity field from SIR-C/X-SAR during 3 days (7-10 October 1994), while we referred to several images in more years, and (ii) CR in [27] come from an assessment of surface flow velocity at the front of the PMG, while we referred to PMG front and section TB62 further upstream (including ablation).

Discussion and Benchmark against Recent Studies
Among others, Ciappa et al. [7] estimated ice velocity in two cross sections of PMG's front during 2009, using Cosmo-Skymed images at 1 m resolution, processed with the maximum correlation approach. They found flow velocity very similar to our values here, i.e., about 2 m day −1 (main flow line) at flow peak season (austral summer, December-January-February). Moreover, they reported seasonal dependence as we found here, i.e., the warmer the faster (see Figure 8 therein). More recently, Minowa et al. [22] provided a benchmark for analysis of PMG flow features. They used normalized cross correlation of spatial variations in image intensity (e.g., [38]), also validated using GPS positioning of an ice drilled pole (upstream section TA1 in Figure 3b here) during some periods in 2010-2014.
They provide surface flow velocity mapping ( Figure 5 in [22]) somewhat consistent during the (austral) winter and summer.
Their estimated velocities vary from about Vs = 4 m day −1 in the upper (Southern) part of the glacier (winter accumulation zone, about 120 m a.s.l.), to about Vs = 2 m day −1 at the front, along the main (central) flow line. Lower flow velocity is mapped in the ablation tongue along the borders (V S = 0-1 m day −1 ).
Such estimates are well consistent with those from our values here (see Figure 3b, covering the same areas of Figure 5 in [22]).
We found Vs = 4 m day −1 in the South (accumulation) area of the PMG, with about 2 m day −1 in the calving front instead. Moreover, here along the edges of the calving front, lower velocities are found.
In Figure 8, we report seasonal variability of surface flow velocity. We report averaged velocity Vs along transect TB62 (North to South) during the summer (December-January-February, average of 9 couples of images from SENTINEL2), fall (March-April-May, average of 9 couples of images), and winter (June-July-August, average of 5 couples of images). Spring (September-October-November) is also reported, but it is likely less accurate because we had only a couple of images available.
In general, higher velocities occur in the summer/fall with warmer temperature, and lower velocities occur in the winter/spring with colder weather. This is in agreement with presently available literature on the PMG and especially [22], where surface velocity at the front fluctuates by ±25% or so (about 1.5-2.5 m day −1 Vs about 2 m day −1 on average), between summer and winter (see Figure 8 therein).
Moreover, from benchmarking of the estimated velocity profiles against approximate estimation of glaciers depth (altitude of surface profile minus altitude of bedrock profile, see  the largest (surface) velocity occurs mostly in the deepest areas, as expected. Figure 9 demonstrates flow velocity change along the glaciers from 2015 to 2017 using SENTINEL2 images. We calculated mean velocity across the glacier on transects starting from the front and moving backward with steps of 230 m, until 7.5 km upstream, at about the TA3 section in Figure 3b. We report the mean flow velocity in each transect (16 sections), together with the local deviation in each section for each acquisition date, colored seasonally. The largest Vs is seen at the front, and it fluctuates largely between different images (overall µ VS = 1.10 m s −1 , with σ VS = 0.44 m s −1 ). Large variability in seen, but in general, Vs increases with the season, and at the front, low values are seen especially in winter and spring (i.e., with largely smaller negative fluctuations). However, more images would be possibly needed to fully assess seasonal variations.
Minowa et al. [22] also reported ( Figure 8 therein) estimated mean calving rate (1999-2013) in m day −1 (calculated in practice like in our Equation (4) here; they set CR = . c + . m, i.e., flow plus net advance at glacier's front, which they call total ablation rate). By taking their constant front section area as A m = 0.68 km 2 (Equation (3) therein), they found a maximum in summer (February) of about CR Max = (3.25 ± 1) × 0.68 × 10 6 = 2.21 × 10 6 ± 6.80 × 10 5 m 3 day −1 and a minimum in winter (August) of about CR Min = (0.75 ± 1) × 0.68 × 10 6 = 5.1 × 10 5 ± 6.80 × 10 5 m 3 day −1 , i.e., with a four times bigger rate in summer than in winter. Here, from Minowa et al. [22] estimated on average CR Av,m/d = 1.3 ± 0.7 m day −1 , or CR Av = 8.84 × 10 5 ± 4.75 × 10 5 m 3 day −1 , which is coherent with our estimates here (CR M1,Av = 7.72 × 10 5 ± 2.46 × 10 5 m 3 day −1 , and CR M2,Av = 8.76 × 10 5 ± 2.72 × 10 5 m 3 day −1 ) and with those from [27], CR Av = 1.08 × 10 6 m 3 day −1 . Here, calving rate based on method M2 is mostly larger than that with method M1. Given that section B used for ice flow assessment with method M2 is 7.5 km upstream of the front, the effects due to lag and ablation en route may affect calving rate assessment. The largest differences seem to occur in the coldest seasons. We used here seasonal ablation estimates from [27] from 1995 to 2003 (Table 1 thereby), which may not be representative for our period of study. PMG ablation obviously requires specific assessment, which we could not set up here in lack of ablation measurements, and supraglacial climate data. However, one may take here that assessment as close as possible to the front is possibly more dependable. Within our better measured period (i.e., from 2013 onwards), we were able to observe four events of damming (i.e., with the PMG glacier hitting the opposite shore and clogging the Lago Argentino lake with and ice dam between Canal de los Témpanos and Brazo Rico, Figure 1b), namely, in November 2013, January-February 2016, August-October 2016, and July-November 2017, i.e., slightly less than one event per year. Minowa et al. [22] report four events (with duration of about 6 months) from 2004 to 2012 and three from 2006 to 2012, i.e., about 1 event every two years. This circumstance requires further proper assessment, to evaluate a potential increase of the damming frequency, maybe under modified climate recently in the Patagonia region.
Bown et al. [39] recently used a set of satellite images, including SENTINEL2, and some weather data nearby to map ice flow and mass balance of the Jorge Montt Glacier, Southern Patagonia Icefield, about 230 km N of the PMG, with specific focus on the period 2011-2018. The authors found high ice velocities and strong calving, altogether resulting in high mass imbalance, mostly attributed to frontal ablation. However, according to the authors, such glacier behaves differently from other glaciers in the area, such as the PMG, and comparison may only be qualitative.
Calving rate assessment, which gives possibly the most interesting results, clearly suffers from (poor) estimation of its terms (namely, flow rate, front position and ablation), all subject to considerable noise. Analysis of noise in calving rates however would require assessment of flow velocity as estimated via satellite, against those assessed in the same periods from ground observations, which was not possible here, if not only qualitatively. Moreover, actual front positioning from local assessment would be probably necessary. Ablation would also need to be measured on the ground in some points and maybe modelled accurately glacier-wise. All such pieces of information are sparsely available.
Here, for the purpose of the paper, we could state that our estimates reasonably well match those from other authors in other periods and thereby seemingly indicate substantial stability of the glacier.
Overall, very little weather and glaciological local information is available for the PMG, and full understanding of all mechanisms leading to glacier's dynamics and stability is difficult without further field investigation. Notwithstanding, our findings are of interest given that we investigated for the first time PMG with SENTINEL2 data that we know of, and we extended the present knowledge by 4 years as given in the literature. Recently, ICESat-2 satellite, a NASA follow-up mission to ICE-Sat (2003-2010, Res. 60-70 m) was launched, with the goal to continue measuring and monitoring the impacts of the changing environment. The ICESat-2 satellite carries an improved laser altimeter called ATLAS (Advanced Topographic Laser Altimeter System), designed to measure ice-sheet topography, sea ice freeboard as well as cloud and atmospheric properties and global vegetation, with a resolution of about 20 m, i.e., slightly coarser than, but comparable to, Landsat7-8 and Sentinel 2 here. In the future, some attempts at using also data from this satellite may be pursued, including for the PMG area, to gain further insight.

Conclusions
The worldwide renowned Perito Moreno glacier is a paradigmatic one to assess stability of Patagonian ice bodies under climate change. However, meteorological information is sporadic, and especially ice dynamics (flow and ablation) are very little explored, given complex local conditions, preventing in situ investigation.
PMG displayed anomalous glaciological behavior, because its front position did not significantly vary for many years. Periodically, its front moves forward (by about 100 m) to reach the Peninsùla de Magallanes, damming the Lago Argentino with subsequent ice fall [23], a phenomenon followed worldwide. This normally happens in the austral summer during March, after lifting of the water level in Brazo Rico by +10 m or more, and it is supposed that changes in such dynamics may be in response to the locally modified climate.
Accordingly, one way to investigate changes in the PMG dynamics is to observe flow (velocity/calving) dynamics using remote sensing.
We used the most recent and remote sensing tools with the finest resolution, namely LANDSAT7-8 and SENTINEL2, which we processed using ImGRAFT Templatematch, an open source tool, to obtain recent (2001-2017) and especially current (until 2016-2017, unmapped hitherto), and finest scale (10 m from 2016 to 2017) PMG flow velocity and calving rate.
The highest altitude accumulation area of PMG is often cloud covered, and the ice surface is mostly smooth due to snow cover, making velocity assessment based upon feature pattern recognition complex therein. We could obtain best results here, originally using couples of images from SENTINEL2, somewhat confirming the fitness of this satellite constellation for ice flow tracking [25,26].
A comparison with ice velocity measurements on the ground [27] reveals acceptable accuracy of the method, especially when using fine resolution (≤15 m) satellites, and comparison with recent satellite mapping exercises [22] further gives ground to our results.
Seasonal dependence of flow velocity (i.e., the warmer the faster) is also confirmed, albeit requiring further investigation.
Overall, we could gather that calving rate is possibly better estimated by focusing upon the PMG front, rather than on upstream areas, given uncertainty coming from poorly known ablation rates en route to the front.
Eventually, our results, in comparison with recent literature, seemingly indicate the present constant dynamics of the PMG in terms of flow velocity and calving because no clearly faster or lower Vs or CR could be highlighted recently. However here, some uncertainty remains, which has the potential to hide possible present (or future) changes. We found, as reported, more events of damming than normally expected (i.e., four events from 2013 to 2017), which may be worthy of further investigation.
Eventually, our preliminary exploration here demonstrates that SENTINEL2 and LANDSAT8, maybe interpreted using ImGRAFT TemplateMatch as we did here or similar methods, could be used henceforward to monitor PMG glacier dynamics on the whole glacier body, but especially focusing on the calving front. Finer resolution (e.g., local, UAV based) images may help refining glacier's flow analysis at smaller spatial scales than here. In March 2017, we acquired drone-based images, covering a small patch on the South ablation tongue of PMG, to assess flow velocity, with spatial resolution in the order of few centimeters. Preliminary analysis therein indicated components of flow velocity towards glacier's (Southern) edge, i.e., a predominant crosswise flow direction nearby the glacier's borders. Such components cannot be highlighted using satellite images, which suggests that fine scale measurements of ice flow velocity/direction are valuable, whenever one wants to gather a clearer understanding of PMG's (as of any other glaciers') dynamics.
Little weather and glaciological local information is available for the PMG, and full understanding of all mechanisms leading to the glacier's dynamics and stability is therefore difficult nowadays. However, our findings here seem to confirm the present stability of the glacier and are of interest given that we investigated for the first time PMG with SENTINEL2 data that we know of, and we extended by 4 years glacier's flow and calving analysis, thus contributing to an increased knowledge of the PMG.