回复审稿人#1：

我们感谢审稿人提出的详细意见，并已根据所有建议进行了修改。在本回复中，审稿人的意见以黑色标准字体显示。我们的回复以蓝色标准字体显示，对稿件的修改以蓝色粗体显示。

主要评论：

1）在第3.2节中，无人机覆盖范围似乎仅限于研究区域内单个冰川的一小部分末端区域。请阐明该训练区域在空间和时间（季节性/全年条件）上的代表性程度，并论证将无人机校准的参数/权重外推至整个区域和完整年度周期时的适用性。

感谢您的评论。我们已在3.2节中阐明，尽管无人机覆盖范围仅限于亚农冰川末端附近约30平方公里的区域，但无人机获取的数据涵盖了以下几个方面：(i) 从消融季到早期积累期（2023年6月至11月）的季节性条件；(ii) 从快速流动的冰川干流到流动较慢的边缘区域，呈现出清晰的空间速度梯度；(iii) 包括碎屑覆盖冰、洁净冰和裂缝带在内的多种表面相。这些特征提供了异质但相关的条件，使得无人机数据子集具有潜在的通用性，可用于评估Landsat、Sentinel-1和Sentinel-2速度反演的性能，以及校准本研究中使用的融合权重。

此外，我们使用配备M6 Pro (M6P)测绘相机的DJI Matrice 300 RTK无人机获取了用于参考和评估的无人机摄影测量数据。2023年6月至11月期间，我们采集了六幅正射影像，构成五组影像对。虽然覆盖范围仅限于亚农冰川末端附近约30平方公里的区域，但无人机测量涵盖了多种不同的冰川条件：(i) 测量时间段涵盖了消融期到早期积累期的过渡阶段；(ii) 测绘区域既包括快速流动的冰川干流区域，也包括流速较慢的冰川边缘；(iii) 表面类型包括覆盖碎屑的冰、洁净冰和裂缝区域。因此，这组无人机数据在空间和时间上具有潜在的可迁移性，可用于评估Landsat、Sentinel-1和Sentinel-2冰川速度反演的性能，以及校准区域融合中的融合权重。

2) 虽然稿件中指出来自三个数据源的速度数据已统一调整至 30 米的空间分辨率，但不同传感器/产品的底层原始网格未必已配准。为了确保逐像素加权最小二乘法 (WLS) 拟合及后续融合的有效性和可重复性，应明确描述融合前的配准/统一步骤，包括参考网格的选择、重投影和重采样方法（以及插值方案）、网格对齐策略、无数据/掩码传播规则，以及这些操作是否在 WLS 拟合之前执行。

Thank you for this helpful suggestion. We agree that explicit description of the pre-fusion grid harmonization is necessary for ensuring the validity and reproducibility of the pixel-wise WLS fitting and subsequent fusion. We have therefore added a brief statement in Sect. 3.3 clarifying (i) the reference grid used for harmonization, (ii) the resampling method adopted for co-registration, and (iii) that these operations are performed prior to the WLS fitting.

“To support pixel-wise WLS fitting and the subsequent fusion, all velocity products were harmonized to a common 30 m reference grid prior to WLS. We used the Sentinel-2 COSI-Corr velocity image as the reference geometry, and co-registered the Landsat-, Sentinel-1-, and UAV-derived velocity image to this grid using nearest-neighbour resampling, so that the original velocity values are preserved without interpolation smoothing.”

3) For the WLS-based fusion weights, please specify the exact objective function, the interpretation ofthe weights within the WLS formulation, and any constraints imposed during optimization. In particular, clarify whether non-negativity and/or a normalization constraint are enforced. If constraints are used, briefly describe the corresponding solution strategy/implementation to ensure reproducibility.

Thank you for this comment. We have clarified the WLS formulation by explicitly stating the constraints imposed on the fusion weights. In our implementation, the weights are constrained to be non-negative and to sum to one.

“Let  denote the UAV-derived velocity at pixel (, and ,, denote the Landsat, Sentinel-1, and Sentinel-2 velocities, with corresponding fusion weights ,,. The WLS objective is:

After defining the objective in Eq. (1), we estimate the fusion weights by minimizing  subject to the constraints , ,  and . This simplex constraint ensures that the weights are directly interpretable as non-negative mixing fractions. The resulting optimal weight triplet is then applied to fuse the remaining monthly velocity maps.”

4) In Sect. 3.4, (𝑖,𝑗)is defined as the pixel location in the image. However, in Eq. (1) the weight estimation would normally iterate over all pixels of the 2-D grid, whereas the summation appears to run only over 𝑖. Sentinel-1–derived velocities are used as the information source for filling gaps in the fused product (via the enhancement-coefficient field), but the rationale for selecting Sentinel-1 as the gap-filling reference is not sufficiently articulated. Also, the first two paragraph of Section 3.4 is a bit repetitive with Section 3.1.

Thank you for these helpful comments. We have revised Sect. 3.4 accordingly.

(1) We adjusted Eq. (1) and its notation to clearly indicate that the WLS objective is summed over the set of overlapping valid pixels on the 2-D grid.

(2) Before introducing the enhancement-coefficient infilling, we added a short rationale explaining why Sentinel-1 is used to guide gap repair: optical (Landsat/Sentinel-2) velocities are often missing due to cloud/snow/illumination limitations, whereas Sentinel-1 provides a more complete auxiliary field, helping preserve glacier-motion patterns while improving continuity.

(3) We rewrote the first two paragraphs of Sect. 3.4 to reduce redundancy with Sect. 3.1 and to focus the section on the weight-estimation formulation and subsequent steps.

“In this section, we describe how fusion weights are estimated using co-temporal UAV-derived velocities. Monthly velocity maps from Landsat, Sentinel-1, and Sentinel-2 are generated independently, and UAV orthomosaics are processed to provide a high-accuracy reference for the same periods. A WLS formulation is then used to estimate the optimal weights for the three satellite products, which are applied in the subsequent fusion.

Let  denote the UAV-derived velocity at pixel (, and ,, denote the Landsat, Sentinel-1, and Sentinel-2 velocities, with corresponding fusion weights ,,. The WLS objective is:

After defining the objective in Eq. (1), we estimate the fusion weights by minimizing  subject to the constraints , ,  and . This simplex constraint ensures that the weights are directly interpretable as non-negative mixing fractions. The resulting optimal weight triplet is then applied to fuse the remaining monthly velocity maps.

In the fusion stage, to address pixels with missing values (NoData), we introduce a binary mask  to locally renormalize the weights, so that the weighted average at ( is computed only over sources that are valid there (Landsat or Sentinel-2). The preliminary fused velocity is:

Although weighted fusion can effectively integrate multi-source information, some areas may still contain NoData. To further fill these gaps and enhance data continuity, this study introduces a sliding-window enhancement-coefficient infilling method. Because Landsat- and Sentinel-2–derived velocity maps often exhibit large gaps in this cloud- and snow-prone region and simple interpolation can be unreliable, we use Sentinel-1 SAR velocities—less sensitive to clouds and illumination—as a more complete auxiliary field to guide gap repair. Specifically, Sentinel-1 is used to infer the local spatial variation of the fused field and fill NoData accordingly: for each pixel , we consider a 10×10 sliding window  (stride = 1 pixel) and estimate a local enhancement coefficient  that describes the response of the fused value to Sentinel-1,:

”

5) In Eqs. (3) and (4), the enhancement coefficient is denoted as aΩ, which naturally reads as a single constant/parameter associated with a domain Ω. However, based on the stated computation and definition (window-based estimation followed by Gaussian smoothing to form a spatial field), the enhancement coefficient should vary spatially and thus be a gridded quantity.

Thank you for pointing this out. We agree that the enhancement coefficient is a spatially varying gridded field rather than a single constant associated with a window domain. We have therefore revised the notation in Sect. 3.4.

“Although weighted fusion can effectively integrate multi-source information, some areas may still contain NoData. To further fill these gaps and enhance data continuity, this study introduces a sliding-window enhancement-coefficient infilling method. Because Landsat- and Sentinel-2–derived velocity maps often exhibit large gaps in this cloud- and snow-prone region and simple interpolation can be unreliable, we use Sentinel-1 SAR velocities—less sensitive to clouds and illumination—as a more complete auxiliary field to guide gap repair. Specifically, Sentinel-1 is used to infer the local spatial variation of the fused field and fill NoData accordingly: for each pixel , we consider a 10×10 sliding window  (stride = 1 pixel) and estimate a local enhancement coefficient  that describes the response of the fused value to Sentinel-1,:

Only pixels where both and  are valid are included in the summation. Next,  is smoothed with a Gaussian filter to form a spatially varying enhancement-coefficient field  over the entire image. Finally, for missing pixels  in the fused map, infilling is performed as

Through this enhancement-based infilling, the spatial completeness of the fused image is markedly improved and gaps are reasonably reconstructed, providing more stable and continuous data support for subsequent glacier-change analyses and time-series modeling.”

6) In Sect. 4.1, the manuscript reports an “~70% / ~30% improvement” in valid pixels, but it is unclear whether this refers to an absolute increase (percentage points) or a relative increase with respect to a baseline. To avoid ambiguity, please provide the explicit definition/formula used to compute the improvement.

Thank you for this helpful comment. We agree that the phrase “~70% / ~30% improvement” is ambiguous. To avoid confusion, we have revised Sect. 4.1 to report the improvement in percentage points (pp) instead of relative percentages.

“Figure 3 presents the monthly fraction of valid pixels in velocity images from 2015–2024 for the two representative glaciers, Xirinongpu and Yanong. The results show that, in most months, the fused product attains a substantially higher valid-pixel ratio than either single-source dataset. Over Xirinongpu, the mean valid-pixel ratios for Landsat and Sentinel-2 are 57.4% and 54.7%, respectively, whereas the fused result reaches 97.9%, corresponding to absolute increases of 40.5 and 43.2 percentage points, respectively. Over Yanong, Landsat and Sentinel-2 achieve 74.7% and 77.1%, while the fused image further increases to 99.2%, corresponding to absolute increases of 24.5 and 22.1 percentage points, respectively. These findings indicate that the proposed fusion method effectively overcomes the limitations of individual sources and markedly enhances data availability and continuity.”

7) In Sect. 3.5, “image smoothness” is defined as the pixelwise standard deviation within the mask for each month (“…we compute the pixelwise standard deviation…”). However, later (Fig. 5 and associated text) the manuscript refers to using variance to assess image smoothness.

Thank you for pointing out this inconsistency. We agree that the terminology should be consistent. In our analysis of image smoothness, the metric used is variance, and the “standard deviation” wording in Sect. 3.5 was imprecise. We have revised Sect. 3.5 to consistently use variance when describing the image-smoothness assessment.

“(3) Image smoothness: For the three velocity products and for GoLIVE and ITS_LIVE, we compute the pixelwise variance within the mask for each month. Changes in variance indicate noise suppression and smoothing, enabling a comprehensive assessment of fusion-driven gains in coverage, temporal continuity, and smoothness.”

8) In Sect. 4.4, the manuscript states that the observed bimodal intra-annual velocity pattern “accords with” the subglacial hydrology-dynamics evolution framework for maritime glaciers. While this interpretation is plausible, the current discussion remains largely qualitative and lacks a clear evidential link to the data presented in this study.

Thank you for this important comment. We agree that the original wording could be interpreted as an overly strong attribution. In the revised manuscript, we have moderated the interpretation and reframed it as a plausible process-based explanation that is broadly consistent with the seasonal pattern observed in our data and with previous studies. We also added a citation to Nanni et al. (2023) to support that similar seasonal acceleration behaviour, including autumn acceleration, has been reported in mountain glaciers in other regions.

“To characterize intra-annual variability in the Kangri Karpo region, we extracted monthly velocities from the 2015–2024 fused product and averaged each calendar month across years to derive a climatological annual cycle (Fig. 11). The resulting series shows a clear seasonal signal: velocities are relatively low in January-March, increase from March onward, reach a first peak in May, remain comparatively stable during summer, rise again to a second peak in October, and then decline. This indicates a bimodal intra-annual velocity pattern in the study region.

This bimodal pattern is broadly consistent with process-based interpretations proposed for temperate/maritime glaciers, in which seasonal changes in subglacial drainage efficiency modulate basal sliding. In particular, early-season meltwater and rainfall inputs may enhance sliding under a relatively inefficient distributed drainage system, whereas subsequent drainage-system evolution toward more efficient channelized flow may reduce mean basal water pressure and limit velocities later in the melt season; partial re-closure or renewed pressurization may contribute to a later-season acceleration. Similar seasonal acceleration behaviour, including autumn acceleration, has also been reported for mountain glaciers in other regions (e.g., Nanni et al., 2023). However, we emphasize that the present study does not include direct observations of subglacial hydrology, and therefore this interpretation should be viewed as plausible rather than diagnostic.

The standard-deviation envelope in Fig. 11 further suggests month-dependent interannual variability. Standard deviations are larger in April-May, when seasonal acceleration initiates, indicating stronger year-to-year differences in the timing and magnitude of acceleration. By contrast, variability is comparatively smaller from June to the following March. Taken together, these results indicate a clear seasonal velocity response in the Kangri Karpo region, with a recurrent bimodal intra-annual structure.”

Nanni, U., Scherler, D., Ayoub, F., Millan, R., Herman, F., and Avouac, J.-P.: Climatic control on seasonal variations in mountain glacier surface velocity, The Cryosphere, 17, 1567–1583, https://doi.org/10.5194/tc-17-1567-2023, 2023.

9) Based on the velocity formulation described in Sect. 3.3, the primary unit of the derived glacier surface velocity should be m d-1. While, in several places (including Sect. 4.3) the manuscript reports velocities in m yr-1 or m/year. Please standardize the unit system throughout the manuscript.

Thank you for this helpful comment. We agree that the velocity unit should be standardized throughout the manuscript. The primary unit directly derived from the velocity formulation in Sect. 3.3 is m d⁻¹, and in the revised manuscript we have unified the velocity unit accordingly. Specifically, we revised the text, figure labels, and related descriptions (including Sect. 4.3 and associated figures) to consistently report glacier surface velocity in m d⁻¹ and removed the previous mixed usage of m yr⁻¹ / m/year.

Figure 7: Mean glacier surface velocity in the Kangri Karpo region (2015–2024).

Figure 8: Correlation between glacier area and velocity. Blue points denote glaciers in the study area; the red solid line is the fitted regression; axes are logarithmic.

Figure 9: (a), (c) Variations of glacier surface velocity and slope along the centerlines of the Yanong and Xirinongpu glaciers. Blue curve: velocity; red curve: slope. (b), (d) Scatterplots of velocity versus slope for all pixels within the Yanong and Xirinongpu glacier masks; red dashed line: linear fit of the slope-velocity relationship.

Specific comments

Title: “High resolution” could mean high spatial resolution and temporal resolution. As you indicated “monthly glacier surface velocity”, how about using “30-meter” to replace “highresolution”?

Thank you for this helpful suggestion. We agree that the term “high resolution” may be ambiguous because it can refer to either spatial or temporal resolution. To make the title more specific, we have revised it by replacing “high-resolution” with “30-meter”.

Revised title

30-meter monthly glacier surface velocity mapping in the Kangri Karpo region (2015–2024) using multi-source remote sensing data fusion

Line 17: This sentence confused me the fusion method is for whether high spatial resolution or high spatial and temporal resolution?

Thank you for this comment. We agree that the original wording was ambiguous and could be interpreted as referring to both spatial and temporal resolution. In the revised manuscript, we clarified this point by replacing “high-resolution” with “high-spatial-resolution” in the relevant sentence (Line 17) and throughout the manuscript where appropriate.

Line 25: It would be better to state as “with similar area and slope, south-facing glaciers are slightly faster than northern ones”.

Modified accordingly.

“…,within individual glaciers, velocity responds more strongly to slope; and, with similar area and slope, south-facing glaciers are slightly faster than north-facing ones.”

Line 96: Please add detail precipitation of this area to illustrate this area is “among the wettest sectors”. You also mentioned in Line 99 “high annual precipitation and humidity”, which may be repetition in meaning.

Thank you for this comment. We agree that the expression “among the Plateau’s wettest sectors” would be stronger if supported by explicit precipitation values. As the cited reference (Wu et al., 2021) provides a qualitative description rather than site-specific precipitation numbers for our study area, we avoided introducing unsupported quantitative values. Instead, we revised the sentence to a more cautious qualitative wording (“humid, monsoon-influenced sector”) and removed the potentially repetitive phrasing elsewhere in the manuscript. And by replacing the repeated climate description with a more specific statement on monsoon influence, seasonal hydroclimatic variability, and frequent cloud cover.

“Located at the eastern terminus of the Nyainqêntanglha Range on the southeastern Qinghai–Tibet Plateau, the Kangri Karpo region is a relatively wet sector of the Plateau (Wu et al., 2021).”

“Climatically, the area lies within the monsoonal temperate glacier region and is strongly influenced by the Indian monsoon (Yang et al., 2008), with pronounced seasonal hydroclimatic variability and frequent cloud cover.”

Figure 1: “KM” or “km”, but not “Km”.

Modified accordingly.

Figure 2: “Weight” not “Weigth”. Also, the use of blue and white in the figure is a bit misleading. UAV velocity image and Landsat velocity image are in white, which means they are processing and analysis steps as you indicate in the caption. What are the input datasets of them? Why the final outputs of the velocity results (2015-2024) are not in blue?

Modified accordingly.

Line 159: Please provide the detail acquisition time of the six UAV images, rather than listing

a time span of June to November in 2023.

Thank you for this suggestion. We have revised the manuscript to report the exact acquisition dates of the six UAV orthomosaics (8 June, 4 July, 10 August, 1 September, 3 October, and 20 November 2023), instead of only giving the June–November 2023 time span.

“In addition, we acquired UAV photogrammetry for reference and evaluation using a DJI Matrice 300 RTK equipped with an M6 Pro (M6P) metric mapping camera. Six orthomosaics were collected on 8 June, 4 July, 10 August, 1 September, 3 October, and 20 November 2023, forming five image pairs.”

Line 166: “32×32 pixels” not “32*32 pixels”. The results of “2 pixels for Landsat OLI” and “3 pixels for Sentinel-2 MSI” have the same spatial resolution. Why did you not set the sept size as 1 pixel for both Landsat and Sentinel-2, and resample the displacement result to 30m afterwards?

Thank you for this helpful comment. We have corrected the formatting from “32*32 pixels” to “32×32 pixels.”

关于步长设置，我们特意为Landsat OLI使用了2像素步长，为Sentinel-2 MSI使用了3像素步长，以便使位移输出的间距与目标约30米分辨率直接匹配。如果两个传感器都使用1像素步长，然后再重采样到30米，则会显著增加计算成本和处理时间。此外，与先生成密度更高的10/15米位移场，然后再重采样到30米相比，直接获取约30米分辨率的结果可以避免额外的重采样/插值步骤，并减少潜在的平滑效应，这更有利于后续的融合。

第 170 行：我建议在本段中更多地介绍时间序列的稳定性，包括 α 截尾均值滤波器的使用描述、选择 0.33 的阈值标准以及“*”的使用（这可能不是正文中的正式用法）。

感谢您的建议。我们在修改后的稿件中阐明了α = 0.33的设定依据。该参数是根据我们研究区域冰川速度时间序列的观测特征，通过经验方法选定的。在融雪季节（通常为6月至9月），光学遥感观测（Landsat和Sentinel-2）经常受到云雾污染（以及季节性积雪/云影效应）的影响，这会导致推导出的月平均速度出现异常高低的偏差。如果不加以抑制，这些异常值可能会严重干扰后续的融合结果。基于研究区域大部分年份的观测数据，年度月平均速度序列的中间部分提供了更为稳定的参考值，因此我们采用了α = 0.33的对称截尾方法（即截去下三分之一和上三分之一，并对中间约34%的数据取平均值），作为一种实用且稳健的设定。

对于Landsat和Sentinel-2影像，速度是通过COSI-Corr软件，利用频域互相关算法计算得到的。搜索窗口为32×32像素；Landsat OLI的步长为2像素，Sentinel-2 MSI的步长为3像素；相关阈值为0.95。东西向（EW）和南北向（NS）位移合并为总位移，然后除以成对的时间基线，得到30米分辨率的月度冰川表面速度。

为了提高可靠性，我们首先剔除了速度大于 3 md⁻¹ 的不匹配像素以及信噪比 (SNR) 小于 0.9 的低置信度像素。对于光学影像衍生的速度（Landsat 和 Sentinel-2），由于光学匹配更容易受到云/云影污染、季节性积雪、光照变化和地表纹理变化的影响，因此需要进行额外的稳定处理。这些因素可能导致某些月份出现异常大的或小的不匹配值。为此，我们应用了年内 α 截尾均值滤波器来稳定每个像素的月时间序列，并降低残余异常值的影响。具体而言，对于每个像素，我们对给定年份内的 12 个月速度值进行排序，并应用 α = 0.33 的对称截尾，即舍弃大约下 33% 和上 33% 的值，并将中间约 34% 的值取平均值，从而得到代表该年份的参考值。我们选择 α = 0.33 作为 12 个月年度序列的实用稳健设置，因为它既能抑制极端值，又能保留中心子集以实现稳定估计。然后，我们将月度值与该参考值逐像素进行比较；大于参考值 1.5 倍或小于参考值 0.5 倍的值被标记为异常值并移除。

第 311 行：“0.10，低于”而不是“0.10-低于”

已据此修改。

第 321 行：“product” 而不是 “produc”

已据此修改。

第 396 行：请使用连接号 (en dash)。“January–March” 而不是 “January-March”。

已据此修改。

Response to Reviewer #2,

We thank reviewer for the detailed reviews, and we made all suggested corrections. In this response, the reviewer’s comments are in black standard font. Our response is in standard blue font and the modifications to the manuscript are in blue bold font.

Major comments:

1) I strongly suggest switching velocities to meters per year rather than meters per day, for three reasons. Firstly, the velocities are generated from monthly pairs, yet their values and uncertainties are reported in meters per day. Although the conversion from one unit to the other is a simple operation, it is confusing to use meters per day because the temporal resolution (monthly) does not allow for the retrieval of results in a unit so small compared to the real measurement (meters per month). Reporting meters per day assumes that the velocity observed can be discretized in days, while it is not. The second reason is that in the manuscript, the authors also refer to meters per year in some cases. This introduces more confusion, and they should adopt a unique, coherent unit throughout their study. Finally, it is my personal opinion that meters per year are a more relatable unit as they allow for spotting outliers more instinctively in the dataset. For example, an uncertainty on stable terrain of 0.05 m per day appears unrealistically small, but in meters per year (0.05 m/d = 18.25 m/yr) this metric appears more reasonable.

Thank you for this important suggestion. We agree that expressing velocities derived from monthly image pairs in m d⁻¹ can be confusing, because the retrieved quantity represents the mean velocity over the full image-pair interval rather than a discretely observed daily motion. We also agree that a single and coherent unit should be used throughout the manuscript.

In response, we have revised the manuscript to report glacier surface velocity consistently in m yr⁻¹ throughout the text, figures, tables, and uncertainty descriptions. We have also revised Sect. 3.3 to clarify that displacement is first divided by the temporal separation of each image pair to obtain an interval-mean velocity, which is then expressed in m yr⁻¹ for reporting consistency and easier interpretation. In addition, all trend units have been converted accordingly (e.g., from m d⁻¹ decade⁻¹ to m yr⁻¹ decade⁻¹).

“For Landsat and Sentinel-2 imagery, velocities were derived with COSI-Corr using a frequency-domain cross-correlation algorithm. The search window was 32×32 pixels; the step size was 2 pixels for Landsat OLI and 3 pixels for Sentinel-2 MSI; the correlation threshold was 0.95. East-west (E-W) and north-south (N-S) displacements were combined to form total displacement, which was then divided by the pairwise time baseline to derive the interval-mean glacier surface velocity for each image pair. Because the retrieved quantity represents the mean velocity over the full image-pair interval rather than discretely observed daily motion, all velocities reported in this study are expressed in m yr⁻¹. The resulting monthly velocity products were generated at 30 m spatial resolution.”

These revisions remove the previous unit inconsistency and improve the interpretability of the reported velocities and uncertainties.

2) I don't know if I missed it in the text, but it seems that the authors do not convert the different remote sensing datasets to the same reference grid. This can yield important impacts on the results as the workflow only accounts for individual pixels. If velocities are generated on different grids, the data fusion will include overlaps and inconsistencies between the grids, potentially skewing the results especially for slow (stable) areas, which are the backbone of the uncertainty analysis and of the coregistration. I suggest that the authors resample each dataset to a common grid in order to homogenize the spatial location of each value.

Thank you for this helpful suggestion. We agree that explicit description of the pre-fusion grid harmonization is necessary for ensuring the validity and reproducibility of the pixel-wise WLS fitting and subsequent fusion. We have therefore added a brief statement in Sect. 3.3 clarifying (i) the reference grid used for harmonization, (ii) the resampling method adopted for co-registration, and (iii) that these operations are performed prior to the WLS fitting.

“To support pixel-wise WLS fitting and the subsequent fusion, all velocity products were harmonized to a common 30 m reference grid prior to WLS. We used the Sentinel-2 COSI-Corr velocity image as the reference geometry, and co-registered the Landsat-, Sentinel-1-, and UAV-derived velocity image to this grid using nearest-neighbour resampling, so that the original velocity values are preserved without interpolation smoothing.”

3) I am not entirely convinced by the selection of Sentinel-1 to fill the data gaps. I do not understand why Sentinel-1 is not included more in the analysis (L.221 provides some information), but it is used to fill the gaps. I understood that it is considered "unreliable" on steep terrain, but we can make the argument that glaciers usually have a much flatter topography than stable terrain. I would like to see more explanations concerning the choice of Sentinel-1 for the gap-filling part.

Thank you for this valuable comment. We agree that the rationale for using Sentinel-1 in the gap-filling step needed to be explained more clearly. In the revised manuscript, we have strengthened this explanation by explicitly noting that, although Sentinel-1 offset tracking is generally less reliable in steep and deeply incised mountain terrain, glacier surfaces are usually smoother and less topographically extreme than the surrounding stable slopes. Therefore, Sentinel-1 remains suitable as an auxiliary source for filling gaps within glacierized areas. We have incorporated this reasoning into Sect. 3.4 and clarified that Sentinel-1 is used as an auxiliary constraint for gap repair rather than as the primary source for direct interpretation of the final fused velocities.

“Although weighted fusion can effectively integrate multi-source information, some areas may still contain NoData. To further fill these gaps and enhance data continuity, this study introduces a sliding-window enhancement-coefficient infilling method. Because Landsat- and Sentinel-2-derived velocity maps often exhibit large gaps in this cloud- and snow-prone region, and simple interpolation can be unreliable in areas with spatially variable glacier motion, we use Sentinel-1 SAR velocities as a more spatially complete auxiliary field to guide gap repair. This choice is also supported by the topographic characteristics of glacierized terrain: although Sentinel-1 offset tracking is generally less reliable in steep and deeply incised mountain areas, glacier surfaces are usually smoother and less topographically extreme than the surrounding stable slopes. Therefore, within glacierized areas, Sentinel-1 can still provide useful motion-pattern information for gap filling, while avoiding an over-reliance on simple interpolation. In this framework, Sentinel-1 is used as an auxiliary constraint for reconstructing missing values rather than as the primary source for direct interpretation of the final fused velocities. Specifically, Sentinel-1 is used to infer the local spatial variation of the fused field and fill NoData accordingly: for each pixel…”

4) I suggest that the authors develop further the discussion about the limitations of the "small" UAV footprint compared to the study area. While their analysis is comprehensive and interesting, its foundations can seem weak based on the fact that the data-fusion workflow relies on the weights evaluated over the UAV's area of interest. This can have important repercussions as it assumes that the relationships (and weights) established for a very narrow area are applicable to the rest of the study area, knowing that this UAV's area of interest does not encompass the same variety of surface parameters (aspect, slope, surface type, glacier type) as the whole study. I am aware that this is a strong limitation, but it needs to be discussed further and stressed more in the discussion.

Thank you for this important comment. We agree that the limited UAV footprint relative to the full study area is a significant constraint of our fusion framework and should be discussed more explicitly. In the previous revision, we clarified in Sect. 3.2 that the UAV surveys, although spatially limited, cover heterogeneous conditions including fast-flowing trunk ice, slower marginal zones, debris-covered and clean-ice surfaces, crevassed areas, and a seasonal window from the ablation season to the early accumulation period. However, we agree that this does not fully resolve the broader issue of representativeness.

In the revised manuscript, we have therefore strengthened the Discussion to explicitly acknowledge that the UAV-calibrated weights are derived from a relatively small reference area near the Yanong Glacier terminus and are then extrapolated to glaciers with potentially different aspect, slope, surface characteristics, glacier types, and seasonal states across the Kangri Karpo region. We now state clearly that this transferability assumption is a methodological limitation and should be kept in mind when interpreting the regional applicability of the fused product. We have also expanded the future-work discussion to note that this limitation could be reduced by using spatiotemporally adaptive weighting strategies and by incorporating additional high-accuracy reference data from multiple glaciers and seasons.

In Sect 3.2：“In addition, we acquired UAV photogrammetry for reference and evaluation using a DJI Matrice 300 RTK equipped with an M6 Pro (M6P) metric mapping camera. Six orthomosaics were collected on 8 June, 4 July, 10 August, 1 September, 3 October, and 20 November 2023, forming five image pairs. Although coverage is limited to ~30 km² near the Yanong Glacier terminus, the UAV surveys span heterogeneous conditions: (i) the period covers the ablation season through the transition to the early accumulation period, (ii) the mapped area includes both fast-flowing trunk regions and slower glacier margins, and (iii) surface types include debris-covered ice, clean ice, and crevassed areas. Therefore, this UAV subset provides a heterogeneous and useful reference for benchmarking the satellite-derived velocity products and for calibrating fusion weights in this study, although its limited spatial footprint means that full regional and year-round representativeness cannot be guaranteed.”

In Sect 4.6：“A key limitation of this study is that the UAV reference data used to calibrate the fusion weights cover only a relatively small area near the Yanong Glacier terminus, rather than the full range of glacier and terrain conditions across the Kangri Karpo region. Although this UAV subset includes heterogeneous conditions, such as fast-flowing trunk ice, slower marginal zones, debris-covered and clean-ice surfaces, and observations spanning the ablation season to the early accumulation period, it does not fully encompass the broader diversity of aspect, slope, surface characteristics, glacier type, and seasonal states represented in the entire study area. Therefore, the present weighting scheme implicitly assumes that the relative performance of Landsat, Sentinel-1, and Sentinel-2 derived from this limited reference area is transferable to other glaciers and time periods, which may not always hold. This limitation should be kept in mind when interpreting the regional applicability of the fused product. Future work should therefore aim to develop spatiotemporally adaptive weighting strategies, for example by calibrating weights separately for different topographic settings, surface conditions, glacier types, and seasons, and by incorporating additional high-accuracy reference data from multiple glaciers and periods.”

5) I would like to address the uncertainty analysis, which is a delicate topic when applied to glacier velocities from remote-sensing. In the ITS_LIVE documentation, the authors acknowledge the strong limitations of uncertainty estimation based on stable ground pixels. They do mention that their propagation to on-ice pixels is not trivial, and vastly underestimates the uncertainties over the ice. This is an issue inherent to any remote-sensing-derived ice surface velocity dataset, and I do not expect this manuscript to come up with a solution that the entire community still fails to provide. However, I would suggest that a more precise uncertainty analysis should be provided, especially considering that the authors have in their possession highly precise UAV data. My suggestion is to focus part of the uncertainty discussion on the specific area covered by the UAV, and compare the uncertainties derived from UAV-only velocities to fused velocities from Sentinel 2 and Landsat. I would also like the authors to acknowledge that the uncertainties from the data-fusion method are likely to be underestimated.

Thank you for this important and constructive comment. We agree that uncertainty estimates derived from stable-ground pixels cannot fully represent the true uncertainty over glacier surfaces, and that their propagation to on-ice pixels is inherently limited. We also agree that, within a data-fusion framework such as ours, the resulting uncertainty may be underestimated if it is evaluated only from stable terrain.

Following your suggestion, we have added an additional uncertainty assessment within the UAV survey area, where high-precision UAV-derived velocities are available as a local on-ice reference. In the revised manuscript, we present a new Figure 7 showing the uncertainty results of the UAV-, Landsat-, Sentinel-2-, and fused velocity products within the UAV-covered area from June to November 2023. The mean 𝐸𝑜𝑓𝑓 values are 2.9 for UAV, 25.9 for Landsat, 21.7 for Sentinel-2, and 18.6 for the fused product. These results show that the fusion method effectively reduces uncertainty relative to the individual optical products, but that a substantial gap still remains between the fused result and the UAV-derived reference. This indicates that data fusion improves the velocity estimates, but does not reach the accuracy level of the high-precision UAV observations.

Based on this comparison, we have revised the manuscript to explicitly acknowledge that uncertainty estimates based on stable-ground pixels, although useful as a regional and long-term reference metric, cannot fully capture the additional error sources present over glacier surfaces, such as surface deformation, crevassing, seasonal snow cover, illumination changes, and texture evolution. We now clearly state that the uncertainty estimated from the present data-fusion framework is likely to be underestimated. We also emphasize that the UAV-footprint analysis serves as an important local on-ice complement to the stable-area uncertainty assessment, rather than a replacement for it.

“Figure 7 presents the uncertainty analysis results for the Landsat-derived velocities, Sentinel-2-derived velocities, UAV-derived velocities, and fused velocities within the UAV survey area. The results show that the fusion method effectively reduces uncertainty relative to the individual optical products. However, a substantial gap still remains between the fused results and the UAV-derived reference. This indicates that, although data fusion can reduce errors to some extent, it still does not achieve the accuracy level of high-precision UAV observations. Moreover, uncertainty estimates based on stable-ground pixels and subsequently propagated to glacier surfaces have inherent limitations. Therefore, the uncertainty obtained from the present data-fusion framework is likely to be underestimated.

Figure 7: Uncertainty analysis of the four velocity products within the UAV survey area from June to November 2023: (a) UAV-derived velocity uncertainty; (b) Landsat-derived velocity uncertainty; (c) Sentinel-2-derived velocity uncertainty; and (d) fused velocity uncertainty.

”

6) My last major comment focuses on the reconciliation of the different temporal baselines of the image pairs used for generating velocities. There is an obvious limitation to generating monthly velocities (or assumed monthly) when the temporal baselines of the image pairs fluctuate by several weeks and rarely correspond to a month exactly. This has the potential for skewing the analysis, knowing that an important part of it focuses on the monthly velocities. Recent work from Charrier et al. (2021 & 2025, complete reference in the pdf) provides tools to correct the overlap between different image-pair baselines and produce regularized velocity time-series. I am aware that the data analysis of the manuscript might have occurred before Charrier et al. published their Python package and manuscript, as I am aware that applying their method to the output of this study might represent a tremendous amount of work. I suggest trying to apply their workflow to the final output of this manuscript, or at least mentioning their method in the discussion of the manuscript (specifically the discussion part about "limitations and future directions"). I do not have a conflict of interest concerning this suggestion. I do however want to mention that I have worked with L.Charrier in the past which is why I am aware of her publication.

Thank you for this very important comment. We agree that reconciling the different temporal baselines of the image pairs is a key issue for nominal monthly glacier-velocity products. In our current workflow, glacier displacement is always divided by the exact temporal separation of each image pair, so the retrieved quantity is the interval-mean velocity for that pair rather than a velocity obtained by assuming an exact one-month baseline. However, we fully agree that differences in revisit cycles among Landsat, Sentinel-1, and Sentinel-2 lead to baseline-length mismatches, temporal overlap, and temporal smoothing, which may affect the strict comparability of the monthly time series and introduce uncertainty into the interpretation of intra-annual variability.

We appreciate your suggestion to consider the regularized time-series framework proposed by Charrier et al. We agree that this is a very valuable direction, as it can better reconcile irregular and overlapping image-pair baselines and generate temporally regularized velocity time series. However, applying that workflow to the full multi-sensor dataset used in this study would require substantial reprocessing and restructuring of the present pipeline, and is beyond the scope of the current revision. Following your suggestion, we have therefore strengthened the discussion in the “Limitations and future directions” section by explicitly acknowledging this limitation and by citing the approach of Charrier et al. (2021, 2025) as an important avenue for future improvement. We now also clarify that the present products should be interpreted as nominal monthly interval-mean velocities rather than strictly calendar-month velocities.

In sect 4.6:” In addition, differences in sensor revisit cycles (Landsat ≈ 16 d; Sentinel-1 ≈ 12 d; Sentinel-2 ≈ 10 d) preclude strict month-to-month alignment, producing baseline-length mismatches, temporal overlap, and temporal smoothing that inevitably introduce additional uncertainty into the nominal monthly velocity series. Although glacier displacement was always divided by the exact temporal separation of each image pair, the resulting products should therefore be interpreted as nominal monthly interval-mean velocities rather than strictly calendar-month velocities. This limitation is particularly relevant for the analysis of intra-annual variability, because irregular and overlapping baselines may blur short-term temporal signals. Future work should explore regularized velocity time-series approaches that explicitly reconcile overlapping image-pair baselines, such as the framework proposed by Charrier et al. (2021, 2025), in order to generate temporally more rigorous and directly comparable glacier-velocity series.”

Detailed comments:

L.65: Since the uncertainty of the reduction is stated and different from the uncertainty of 2015, state the uncertainty of the 1980s result.

Thank you for this comment. We agree that, since the uncertainty of the area reduction is reported, the uncertainty of the 1980s glacier-area estimate should also be stated explicitly. In the revised manuscript, we have corrected this sentence according to the original values reported in Wu et al. (2018). The glacierized area in 1980 is now given as 2728.00 ± 34.24 km², rather than the previously cited inventory value without uncertainty, and the corresponding reduction is consistently reported as 679.51 ± 59.49 km² relative to the 2015 area of 2048.50 ± 48.65 km². The text has been revised accordingly.

L.82: Are topographic shadows also an issue for optical datasets in this region? Especially if we consider small glaciers in steep valleys.

Thank you for this helpful comment. We agree that topographic shadowing is also an important limitation for optical remote sensing in this region, especially for small glaciers located in steep and deeply incised valleys. In the revised manuscript, we have therefore expanded the description of the optical-data limitations to explicitly include terrain-induced illumination effects (e.g., topographic shadowing), in addition to cloud cover and seasonal snow.

“Optical sensors offer high spatial resolution, but their applicability is strongly constrained by persistent cloud cover, seasonal snow, and illumination effects such as topographic shadowing, especially for small glaciers in steep valleys, thereby limiting temporal coverage and data usability (Scherler et al., 2008; Berthier et al., 2005)”

L.124: What is the reference grid? Since the sensors all have their own grid, their derived velocities need to be resampled on a reference grid. Can you detail how you do this step?

Thank you for this comment. We agree that the reference grid and the resampling procedure should be described explicitly, because pixel-wise fusion requires all input velocity products to be aligned to a common grid. In the revised manuscript, we have clarified that, prior to WLS fitting and fusion, all velocity products were harmonized to a common 30 m reference grid. Specifically, we used the Sentinel-2 COSI-Corr velocity field as the reference geometry, and the Landsat-, Sentinel-1-, and UAV-derived velocity fields were co-registered and resampled to this grid using nearest-neighbour resampling, so that the original velocity values were preserved without interpolation smoothing. We have added this explanation explicitly in Sect. 3.3.

“To support pixel-wise WLS fitting and the subsequent fusion, all velocity products were harmonized to a common 30 m reference grid prior to WLS. We used the Sentinel-2 COSI-Corr velocity image as the reference geometry, and co-registered the Landsat-, Sentinel-1-, and UAV-derived velocity image to this grid using nearest-neighbour resampling, so that the original velocity values are preserved without interpolation smoothing.”

L.125: Can you give a quick example in parentheses?

Thank you for this helpful suggestion. We agree that a brief example improves clarity. In the revised manuscript, we have added a short parenthetical example after “basic corrections” to clarify the preprocessing step, namely applying orbit files to Sentinel-1 images prior to offset tracking.

“(1) Image pairing and preprocessing: Assemble monthly image pairs from Landsat-8/9 OLI and Sentinel-2 MSI (optical) and Sentinel-1 IW GRD (radar); perform geometric co-registration and basic corrections(e.g., applying orbit files to Sentinel-1 images).(2)…”

L.130: Can you specify which robust filtering you are using? Or say directly if it's detailed later in the methods.

Thank you for this helpful comment. We agree that the wording here was too brief. In the revised manuscript, we have clarified that the robust filtering method is described in detail later in Sect. 3.3.

“.(4) Temporal robustification and quality control: Remove mismatches using the signal-to-noise ratio (SNR) threshold and a robust time-series filtering method (described in detail later in Sect. 3.3) to improve cross-sensor consistency and reliability.(5)…”

L.131: Can you explain in a few words what it means? I do not understand what an "enhancement-factor field" is in this case. Especially considering "enhancement-factor field" does not appear later in the text.

Thank you for this helpful comment. We agree that the term “enhancement-factor field” was unclear in this summary description. In the revised manuscript, we have replaced it with a more explicit and consistent explanation. We now clarify that this step estimates a local enhancement coefficient from overlapping Sentinel-1 and preliminary fused velocities, smooths this coefficient spatially, and then uses it to reconstruct missing pixels. This wording is now aligned with the terminology used later in the Methods section.

“(6) Sentinel-1-guided gap filling: For residual voids, construct a spatially smoothed enhancement-coefficient field from the relationship between Sentinel-1 and the preliminary fused velocities, and use it to reconstruct missing pixels and improve spatiotemporal continuity.(7)…”

L.147: I would suggest separating this supplementary table into 2: one for Landsat 8 and one for Landsat 9, or at least specify whether you used a mix of the two satellites for the pairs or if you generate pairs only with the same satellite.

Thank you for this helpful suggestion. We agree that the Landsat pairing strategy should be clarified more explicitly. In this study, Landsat velocity pairs were allowed to be generated from both same-satellite and cross-satellite acquisitions, i.e., Landsat-8/Landsat-8, Landsat-9/Landsat-9, and Landsat-8/Landsat-9 pairs, depending on image availability and cloud conditions. In the revised manuscript, we have added an explicit clarification in the main text to state this pairing rule. We retained the existing supplementary table format, but the pairing strategy is now clearly described in the manuscript.

“Detailed information on the Landsat imagery used in this study is provided in Table S1 of the Supplementary Material. Landsat-8 OLI and Landsat-9 OLI-2 were treated jointly rather than as independent sensor streams because of their high inter-sensor consistency (Xu et al., 2024), and Landsat velocity pairs were therefore allowed to be formed from both same-satellite and cross-satellite acquisitions (Landsat-8/Landsat-8, Landsat-9/Landsat-9, and Landsat-8/Landsat-9).”

L.159: Could you add a supplementary table of the exact dates at which the UAV mosaics were collected? It would be consistent with the fact that you provide tables for Sentinel 1 & 2 and Landsat.

Thank you for this helpful suggestion. We agree that the exact acquisition dates of the UAV mosaics should be stated explicitly. Because the number of UAV mosaics used in this study is limited, instead of creating a separate supplementary table, we have added their exact acquisition dates directly in the main text. The revised manuscript now states that six UAV orthomosaics were collected on 8 June, 4 July, 10 August, 1 September, 3 October, and 20 November 2023, forming five image pairs.

L.190: This echoes my comment from L.147: I suggest you write a short sentence about why Landsat 8 and 9 are not considered as independent sensors. Maybe you can cite: Xu, Hanqiu, Mengjie Ren, and Mengjing Lin. 2024. “Cross-Comparison of Landsat-8 and Landsat-9 Data: A Three-Level Approach Based on Underfly Images.” (they seem to show that the two sensors are similar). I am asking because in comparison, ITS_LIVE separates the two missions when using image pairs (although I am not certain why)

Thank you for this helpful suggestion. We agree that the rationale for not treating Landsat-8 and Landsat-9 as independent sensors should be stated more explicitly. In the revised manuscript, we now clarify that Landsat-8 OLI and Landsat-9 OLI-2 were treated jointly in this study because of their high inter-sensor consistency, and because combining them helps increase the number of usable image pairs in this cloud-prone mountain region. We have also added the suggested reference (Xu et al., 2024), which shows a high level of consistency between Landsat-8 and Landsat-9 data.

L.199: See specific comment section.

Modified accordingly.

L.205: Is this the "enhancement-factor field"? If yes, they have different names and this can be confusing.

Thank you for this comment. We agree that using different names for the same concept could be confusing. In the revised manuscript, the terminology has been made consistent throughout, and this term is now referred to uniformly in the text.

L.206: I suggest specifying the standard deviation used for the Gaussian filter, as this variable can have huge impacts on the results of the filter (for reproducibility purposes).

Thank you for this helpful comment. We agree that the Gaussian filter parameter should be reported explicitly for reproducibility. In the revised manuscript, we now specify that the enhancement-coefficient field was smoothed using a Gaussian filter with a standard deviation of 𝜎=2 pixels. We also clarify that nearest boundary handling was used during filtering.

“Next,  is smoothed with a Gaussian filter(σ=2 pixels) to form a spatially varying enhancement-coefficient field  over the entire image. Finally, for missing pixels  in the fused map, infilling is performed as…”

L.226: I suggest using a column : for coherence with the rest of the paragraph.

Modified accordingly.

L.227, 229, 230: Needs a space after the dot.

Modified accordingly.

L.232: This comment might apply to the introduction. It would be worth mentioning that the glaciers are not surge-type, for two reasons: 1) Other glaciers in the area are surge-type (Guillet et al., 2022). 2) A surge-type glacier could cause a great change in the standard deviation, thus invalidating these claims. You could clear any doubt concerning these claims by specifying that the glaciers are not surge-type. This would also help justify that you consider velocities above 1.5* or below 0.5* their reference as outliers.

Thank you for this important comment. We agree that this assumption should be stated explicitly. In the revised manuscript, we now cite the published surge-type glacier inventory for High Mountain Asia (Lv et al., 2022). After checking this dataset against our study area, we found that no glaciers within the study area are identified as surge-type glaciers. This clarification has been added in Sect. 2 and also supports the interpretation of the smoothness analysis and the outlier-screening strategy in Sect. 3.3.

In sect 2:” Glaciers are widely distributed across the massif. According to RGI 7.0, the region contains 218 glaciers (Figure 1): 162 smaller than 1 km², 30 between 1 and 5 km², and 26 larger than 5 km². Among the largest, the Yanong and Xirinongpu glaciers have surface areas of approximately 165 km² and 94 km², respectively (RGI Consortium, 2023). According to the published surge-type glacier inventory for High Mountain Asia (Lv et al., 2022), no glaciers within our study area are identified as surge-type glaciers. Accordingly, surge-type behaviour is not considered in the subsequent analysis.”

In sect. 3.3“To improve reliability, we first removed mismatches with speeds > 1095 m yr⁻¹ and low-confidence pixels with SNR < 0.9. For the optical-image-derived velocities (Landsat and Sentinel-2), additional stabilization is necessary because optical matching is more susceptible to cloud/cloud-shadow contamination, seasonal snow cover, illumination changes, and surface-texture variations, which can produce anomalously large or small mismatches in some months. Because no glaciers within our study area are identified as surge-type glaciers in the published High Mountain Asia inventory of Lv et al. (2022), such abrupt excursions are unlikely to reflect true surge-related dynamics. We therefore applied….”

L.239: How do you select these areas? Using RGI outlines or by manually filtering snow and ice-free pixels?

Thank you for this comment. In this study, the stable areas were selected manually rather than automatically extracted from the RGI outlines. Specifically, we manually identified topographically stable, snow-/ice-free bare-rock areas surrounding the glaciers and excluded zones that showed obvious snow cover, shadow, water, or potential surface instability. We have clarified this point in the revised manuscript.

“Considering the scarcity of long, continuous ground observations in high-mountain regions, we employ a stable-area velocity-fluctuation method, constructing an error model in manually selected topographically stable, snow-/ice-free bare-rock zones surrounding the glaciers to indirectly evaluate the velocity uncertainty of different products.”

L.272: I notice in the supplementary materials that the southern and southwester areas are consistently "flaring up", especially from 2017 to 2022. 2022 also has really high velocities that look like noise and not a coherent pattern. Could you comment on this? It's a non-negligible part of the study area that seems impacted.

Thank you for this careful observation. We agree that the southern and southwestern parts of the study area show persistent local “flaring” patterns in some years, especially between 2017 and 2022, and that the particularly high velocities in 2022 do not always form a fully coherent glaciological pattern.

We interpret this behaviour mainly as a consequence of the challenging observation conditions in that sector. The southern and southwestern parts of the study area are located near mountain passes, where glaciers occupy narrow, steeply incised valleys with very large local elevation differences. In addition, this area is a major corridor for moisture transport, and is therefore frequently affected by persistent cloud and fog cover. Under these combined conditions, all three satellite data sources used in this study are more prone to mismatches, decorrelation, and residual retrieval noise than in other parts of the region.

Although we applied relatively strict outlier-screening and filtering procedures, these measures cannot completely eliminate retrieval artefacts in such particularly unfavorable terrain and atmospheric conditions.

L.310: How would you explain the jumps in variance in the Landsat and Sentinel products that are not replicated by the fusion? Are these outliers or actual ice acceleration/deceleration?

Thank you for this important question. We interpret the abrupt jumps in variance seen in some monthly Landsat and Sentinel-2 products mainly as retrieval artefacts rather than as true glacier acceleration or deceleration signals. In this cloud- and snow-prone mountain region, optical matching is highly susceptible to cloud/fog contamination, shadow, seasonal snow, and local illumination changes, which can introduce anomalous mismatches in certain months and cause abnormally large pixelwise variance.

These variance jumps are not reproduced in the fused product because such anomalous optical mismatches are suppressed during the fusion workflow. Specifically, obvious low-quality matches are first screened using quality-control criteria, and the optical time series is then further stabilized using the robust temporal filtering procedure described in Sect. 3.3. As a result, anomalous values that artificially inflate the variance in individual Landsat or Sentinel-2 products are largely removed before or during fusion, and therefore do not propagate into the final fused result.

L.321: Needs a "t".

Modified accordingly.

L.325: If I understood well, you take the GoLIVE uncertainties directly from the dataset. I see a potential issue for direct comparison, because GoLIVE uncertainties probably do not consider the same off-ice pixels as you do, and might therefore include more outliers and impact the sigma values. It would be worth mentioning this, otherwise the datasets are not directly comparable for their uncertainties.

Thank you for this important comment. We would like to clarify that the uncertainty values compared in this study were calculated within the same type of manually selected off-ice stable areas, rather than directly adopting the uncertainty values provided in the original datasets. In other words, the uncertainty estimates for GoLIVE, ITS_LIVE, and the products generated in this study were all evaluated over manually delineated ice-free stable zones surrounding the glaciers, using a consistent framework.

L.345, 356, 358: I understand that it could be the custom to just write p<0.05 rather than its value. I thought maybe it could add more details about the analysis by allowing for the comparison of the significance between each test. Up to you if you want to change it.

Thank you for this helpful suggestion. We agree that exact p values may provide additional detail for comparing significance levels across tests. However, to maintain a concise presentation and a consistent reporting style throughout the manuscript, we prefer to retain the notation p < 0.05 in the main text.

L.351: It might be worth showing the correlation just for the centerline points as well. It seems that they correlate much better, which is a fascinating result: the correlation between slope and velocity could be a function of the distance to the centerline.

Thank you for this insightful suggestion. We agree that examining the slope–velocity relationship specifically along the glacier centerline could provide additional information, and that the correlation may indeed vary with distance from the centerline. This is an interesting direction for further analysis. However, the present study focuses on the regional-scale relationship between slope and velocity based on the full glacier area, and a dedicated centerline-based analysis would require an additional methodological framework beyond the current scope of the manuscript. We have therefore not added this extra analysis here, but we agree that it is a valuable topic for future work.

L.402: I suggest switching to parentheses because the "-" sign impacts readability.

Modified accordingly.

L.471: You did mention precipitations earlier in the text. Why not include an ERA-5 plot of the average precipitation on the same plot (secondary Y-axis)?

Thank you for this helpful suggestion. We agree that precipitation is also an important climatic factor potentially related to glacier velocity. However, we chose not to add ERA5 precipitation to this figure for two main reasons. First, compared with air temperature, the influence of precipitation on glacier surface velocity is more indirect and complex. In particular, precipitation falling on a glacier does not necessarily or immediately reach the glacier bed; instead, it may be temporarily stored within the snow/firn/ice system, and the associated hydrological response can vary substantially among glaciers. Therefore, placing precipitation on the same plot could easily suggest a simple direct relationship that is difficult to justify in the present context.

Second, precipitation reanalysis products are generally subject to relatively large uncertainties in this high-mountain region. A rigorous precipitation-based interpretation would therefore require additional processing and evaluation (e.g., bias correction and/or downscaling), which is beyond the scope of the present study. For these reasons, we retained the current figure design and focused this part of the discussion on temperature only.

L.485: I think this section could benefit from two points: 1) You should expand on the possible impacts of calibrating the weights on the glacier tongue only, and what it might imply for areas at higher altitude with different surface cover (snow, less debris, different aspects, slope). It is good that you are acknowledging this limitation, but in my opinion this potentially yields stronger limitations on the applicability of those weights (as shown by the consistently noisy results south of the domain, as shown by the mosaics in the supplements for years 2017 to 2022). 2) The second point concerns the adaptive weighting. This is a note I will put in the general comments: some work has been done that would complement yours greatly, by Charrier, L., Dehecq, A., Guo, L., Brun, F., Millan, R., Lioret, N., Copland, L., Maier, N., Dow, C., and Halas, P.: TICOI: an operational Python package to generate regular glacier velocity time series, The Cryosphere, 19, 4555–4583, https://doi.org/10.5194/tc-19-4555-2025, 2025.

I am aware that most of your analysis might have been done before this paper was released. They provide a straightforward and ready-to-use code on their GitHub. I think it addresses some of your future direction points because your results are inherently impacted by the issue of non-alignment of the dates for the monthly pairs, and thus have some overlap between the velocities when you fuse them.

Modified accordingly.

L.504: Inconsistent use of units, I suggest either using only m/d or m/yr. I personally have a preference for m/yr, as detailed in the specific comments.

Modified accordingly.

Response to Reviewer #3,

We thank reviewer for the detailed reviews, and we made all suggested corrections. In this response, the reviewer’s comments are in black standard font. Our response is in standard blue font and the modifications to the manuscript are in blue bold font.

Discussing (or at least mentioning in the introduction) current state-of-the-art fusion algorithms

Introduction : Here you need to mention current state-of-the-art fusion algorithms

Charrier et al., 2025 : h ps://doi.org/10.5194/tc-19-4555-2025, 2025.

Greene et al., 2021 : h ps://doi.org/10.5194/tc-14-4365-2020

Derkacheva et al 2020 : h ps://doi.org/10.3390/rs12121935

Thank you for this helpful suggestion. We agree that the Introduction should better reflect recent state-of-the-art developments in multi-sensor and multi-temporal glacier-velocity integration and regularization. In the revised manuscript, we have expanded the Introduction to include the suggested studies, namely Derkacheva et al. (2020), Greene et al. (2020), and Charrier et al. (2025). These references help position our work more clearly with respect to recent advances in multi-sensor velocity processing, seasonal velocity detection, and temporally regularized glacier-velocity time-series reconstruction.

“Recent studies have further advanced the integration and regularization of multi-sensor glacier-velocity observations. For example, Derkacheva et al. (2020) explored statistical and regression-based approaches for reducing and combining ice-velocity information derived from Landsat-8, Sentinel-1, and Sentinel-2. Greene et al. (2020) developed an approach for detecting seasonal ice dynamics from satellite image pairs, highlighting the importance of resolving sub-annual velocity variability. More recently, Charrier et al. (2025) proposed the TICOI framework to generate temporally regularized glacier-velocity time series from heterogeneous and overlapping multi-sensor observations. These studies demonstrate the growing importance of multi-source integration and time-series regularization for improving the continuity and interpretability of glacier surface velocity products.“

Explaining why velocities were computed between images, spaced by one month only. This is surprising, since it drastically reduces the number of image-pair velocities.

Thank you for this important comment. Our intention was not to maximize the number of available image-pair velocities, but to construct a nominal monthly glacier-velocity product with a consistent temporal meaning. For this reason, we preferentially selected image pairs with an interval close to one month. Using substantially longer baselines would indeed increase the number of valid pairs, but it would also smooth short-term variability and mix velocity signals from different parts of the melt season, thereby reducing the interpretability of the monthly time series. In other words, the one-month strategy was adopted as a compromise between image-pair availability and the temporal resolution required for monthly monitoring.

At the same time, we did not impose a rigid 30-day rule. In practice, we preferentially used low-cloud scenes near the beginning and end of each month, and in months with poor optical conditions the pairwise interval was extended when necessary to maintain data usability.

Give a more detailed explanation of your methodology, including the name of the null hypothesis used, the reason for resampling the images and the presence of gaps after enhancement-based infilling.

Thank you for this helpful suggestion. We agree that the methodological description in the previous version was too concise in several places. In the revised manuscript, we have therefore provided a more detailed explanation of the workflow. Specifically, we now clarify:

(1) the null hypothesis used in the statistical analysis, namely the default two-sided test in scipy.stats.pearsonr, for which the null hypothesis is that the two variables are uncorrelated;

(2) the reason for resampling the images/products, which is to harmonize the different native grids of the multi-sensor datasets onto a common reference grid prior to pixel-wise comparison, coefficient estimation, and fusion; and

(3) the reason why some gaps remain after enhancement-based infilling, namely that the enhancement coefficient is computed within a 10 × 10 sliding window, and when all pixels within a local window are NoData, no valid coefficient can be estimated for that location.

These points have now been clarified explicitly in the revised manuscript to improve reproducibility and reduce ambiguity.

Provide a comparison of your fused velocities with the UAV velocities, even if it’s only covering a part of the glacier, and have been used for validation, this is still a way to validate your results.

Thank you for this helpful suggestion. We agree that, even though the UAV data cover only part of the glacier area, they still provide an important local reference for validating the fused velocity product. Following your suggestion, we have added an additional comparison within the UAV survey area in the revised manuscript (Figure 7).

Figure 7 presents the uncertainty analysis results for the UAV-, Landsat-, Sentinel-2-, and fused velocity products within the UAV-covered area from June to November 2023. The results show that the fused product has lower uncertainty than the individual optical products, indicating that the fusion improves local velocity retrieval quality. At the same time, the fused result still remains notably more uncertain than the UAV-derived reference, which further demonstrates that the UAV comparison provides a meaningful local validation of the fused product. We have incorporated this new analysis and its discussion into the revised manuscript.

In sect 4.2 “Figure 7 presents the uncertainty analysis results for the Landsat-derived velocities, Sentinel-2-derived velocities, UAV-derived velocities, and fused velocities within the UAV survey area. The results show that the fusion method effectively reduces uncertainty relative to the individual optical products. However, a substantial gap still remains between the fused results and the UAV-derived reference. This indicates that, although data fusion can reduce errors to some extent, it still does not achieve the accuracy level of high-precision UAV observations. Moreover, uncertainty estimates based on stable-ground pixels and subsequently propagated to glacier surfaces have inherent limitations. Therefore, the uncertainty obtained from the present data-fusion framework is likely to be underestimated.

Figure 7: Uncertainty analysis of the four velocity products within the UAV survey area from June to November 2023: (a) UAV-derived velocity uncertainty; (b) Landsat-derived velocity uncertainty; (c) Sentinel-2-derived velocity uncertainty; and (d) fused velocity uncertainty.”

1. 128 : do you convert range/azimuth to E-W / N-S component ? if yes, it needs to be specified

Thank you for this helpful comment. In our Sentinel-1 workflow, the offsets were first estimated in SAR image geometry (range/azimuth) using SNAP offset tracking, and the results were then geocoded through terrain correction to the map-projected ground geometry used in this study. We have clarified this processing sequence in the revised manuscript.

“(3) Radar velocity retrieval: Use intensity-based pixel-offset tracking on Sentinel-1 to estimate range/azimuth offsets in SAR image geometry, then apply terrain correction to geocode the results to map-projected ground geometry and convert them to velocity, providing stable coverage under cloudy conditions. “

Fig. 2 : why are you resampling the images ?

Thank you for this comment. The images/products were resampled because the different sensors have different native grids, and a common reference grid is required for the subsequent pixel-wise WLS fitting and fusion. In the revised manuscript, we have clarified this point and added the detailed resampling procedure in Sect. 3.3.

“To support pixel-wise WLS fitting and the subsequent fusion, all velocity products were harmonized to a common 30 m reference grid prior to WLS. We used the Sentinel-2 COSI-Corr velocity image as the reference geometry, and co-registered the Landsat-, Sentinel-1-, and UAV-derived velocity image to this grid using nearest-neighbour resampling, so that the original velocity values are preserved without interpolation smoothing.“

l.174 why 1.5 or 0.5 ? why not  (annual mean + 2 std mean) which is the zscore ?

Thank you for this important question. We did not use a z-score criterion based on the annual mean ± 2 standard deviations because our purpose here was not to detect statistical outliers relative to an assumed approximately normal annual distribution, but to remove obviously anomalous mismatches in monthly glacier-velocity time series derived from optical imagery. In this study, the annual sample size for each pixel is only 12 monthly values, and the distribution may be affected by seasonality as well as by occasional cloud-/snow-related mismatches. Under these conditions, the mean and standard deviation can themselves be sensitive to the outliers that we aim to remove.

We therefore adopted a ratio-based threshold relative to a robust annual reference value derived from the intra-annual α-trimmed mean. The 1.5× and 0.5× thresholds were selected as practical empirical bounds to identify values that deviate strongly from the central, year-representative velocity level while retaining normal seasonal variability. In other words, this procedure is intended as a robust quality-control filter rather than a formal distribution-based anomaly test.

eq 1 is not a weighed least square because you did not weight the residuals. There is also a problem in the summa on : it should be over i and j.

Thank you for this careful comment. We agree that Eq. (1), as originally written, should not have been referred to as weighted least squares. In the revised manuscript, we have therefore renamed the method as a constrained least-squares-based weighted fusion framework and revised the text accordingly. We have also clarified the summation notation explicitly over both 𝑖 and 𝑗.

“In this section, we describe a constrained least-squares-based weighted fusion framework in which fusion weights are estimated using co-temporal UAV-derived velocities. Monthly velocity maps from Landsat, Sentinel-1, and Sentinel-2 are generated independently, and UAV orthomosaics are processed to provide a high-accuracy reference for the same periods. A constrained least-squares formulation is then used to estimate the optimal mixing coefficients for the three satellite products, which are subsequently applied in the weighted fusion of the monthly velocity maps.

Let  denote the UAV-derived velocity at pixel (, and ,, denote the Landsat, Sentinel-1, and Sentinel-2 velocities, with corresponding fusion weights ,,. The constrained least-squares objective is:

After defining the objective in Eq. (1), we estimate the fusion weights by minimizing  subject to the constraints , ,  and . This simplex constraint ensures that the weights are directly interpretable as non-negative mixing fractions. The resulting optimal weight triplet is then applied to fuse the remaining monthly velocity maps.”

Eq 1 do you have one weigth for the entire me-series of each sensor ?

Thank you for this question. Yes, in the present implementation we estimate one set of sensor-level mixing coefficients from the co-temporal UAV-overlap periods, and this same weight triplet is then applied to the monthly velocity maps throughout the time series. In other words, the weights are treated as sensor-specific but temporally invariant in the current framework, rather than being re-estimated separately for each month.

1. 205 Why are you filling the gaps based on Sen nel-1 velocities mainly (even with the enhancement factor) ?

Thank you for this important question. Sentinel-1 is used as the main auxiliary source for gap filling primarily because it provides the most spatially continuous and weather-independent observations in this cloud- and snow-prone mountain region. In contrast, the Landsat- and Sentinel-2-derived velocity products, although generally more reliable under clear conditions, often contain large gaps caused by persistent cloud cover, seasonal snow, and illumination effects.

In our framework, Sentinel-1 is not used to directly replace missing values in the fused product. Instead, it is used to derive a local enhancement coefficient from the relationship between Sentinel-1 and the preliminary fused velocities in overlapping valid areas. This coefficient field is then smoothed and used to guide the reconstruction of missing pixels. In this way, Sentinel-1 contributes spatial continuity and local flow-pattern information, while limiting the direct propagation of its own retrieval errors into the final product.

In addition, although Sentinel-1 offset tracking is generally less reliable in steep and deeply incised mountain terrain, glacier surfaces are usually smoother and less topographically extreme than the surrounding stable slopes. This makes Sentinel-1 suitable as an auxiliary source for gap filling within glacierized areas. We have clarified this rationale in the revised manuscript.

In sect 3.4 “Although weighted fusion can effectively integrate multi-source information, some areas may still contain NoData. To further fill these gaps and enhance data continuity, this study introduces a sliding-window enhancement-coefficient infilling method. Because Landsat- and Sentinel-2-derived velocity maps often exhibit large gaps in this cloud- and snow-prone region, and simple interpolation can be unreliable in areas with spatially variable glacier motion, we use Sentinel-1 SAR velocities as a more spatially complete auxiliary field to guide gap repair. This choice is also supported by the topographic characteristics of glacierized terrain: although Sentinel-1 offset tracking is generally less reliable in steep and deeply incised mountain areas, glacier surfaces are usually smoother and less topographically extreme than the surrounding stable slopes. Therefore, within glacierized areas, Sentinel-1 can still provide useful motion-pattern information for gap filling, while avoiding an over-reliance on simple interpolation. In this framework, Sentinel-1 is used as an auxiliary constraint for reconstructing missing values rather than as the primary source for direct interpretation of the final fused velocities. Specifically, Sentinel-1 is used to…”

1. 220 clouds are also present in Landsat images

Thank you for this helpful comment. We agree that cloud contamination also affects Landsat imagery, and the original wording was not sufficiently precise. In the revised manuscript, we have corrected this sentence to clarify that both Landsat and Sentinel-2 are affected by cloud interference in the study region.

“To visualize quality changes before and after fusion, we compared the fused results separately with velocities from two mainstream optical sources, Landsat and Sentinel-2. Both are widely used for glacier-velocity retrieval, but both are affected by cloud contamination in this region. Such contrasts help verify whether fusion effectively integrates complementary strengths and mitigates the weaknesses of the individual optical products. Note that although Sentinel-1 SAR is included, it is not contrasted separately:”

1. 224 the enhancement factor is used only for Nodata no ? why do you say that sentinel-1 contribution on is further modulated by the enhancement-coefficient ?

Thank you for this comment. We agree that our previous wording was not sufficiently precise. In the revised manuscript, we now clarify that the enhancement-coefficient field is used only in the subsequent gap-filling stage for residual NoData pixels, rather than as a modulation of the standard weighted fusion for all valid pixels. The text has been revised accordingly.

“Note that although Sentinel-1 SAR is included in the framework, it is not contrasted separately: (i) it mainly supplements optical observations during cloudy periods and has limited standalone reliability in steep, deeply incised terrain; and (ii) its contribution to the main weighted-fusion result is limited by its relatively low sensor-level fusion weight (≈ 0.08). The enhancement-coefficient field is applied only afterward, in the gap-filling stage for residual NoData pixels. Accordingly, this section focuses on before–after comparisons for Landsat and Sentinel-2 to highlight systematic improvements to the primary optical observations.”

Fig4 . why do you still have gaps after the enhancement-based infilling ?

Thank you for this question. The enhancement-based infilling in our framework is not able to fill all gaps unconditionally. Specifically, the enhancement coefficient is estimated within a 10 × 10 sliding window. When all pixels within a given window are NoData, no valid local relationship can be established and the enhancement coefficient cannot be estimated for that location. As a result, some residual gaps may remain even after the enhancement-based infilling step.

Fig 5. What’s the unit of the axes ? You need to explain why ITS_LIVE have variance than your fuse products : because you used the annual product probably.

Thank you for this helpful comment. We agree that the axis units should be stated explicitly, and we have revised Figure 5 accordingly.

We also agree that the lower variance of ITS_LIVE relative to our fused product requires clarification. In this comparison, ITS_LIVE represents an annual-scale product, whereas our fused velocities are monthly products. Because the ITS_LIVE annual product effectively integrates velocity information over a much longer temporal window, short-term fluctuations are smoothed, resulting in lower apparent variance. By contrast, our fused product preserves month-to-month variability and seasonal signals more explicitly, which naturally leads to higher variance. We have added this explanation in the revised manuscript.

“Overall, the fused product exhibits a more stable, continuous, and low-fluctuation variance curve in both areas, indicating effective suppression of local noise and improved image consistency. At Yanong, the fused mean variance is 13.2×103 m²·yr⁻²-lower than GoLIVE (30.6×103 m²·yr⁻²), comparable to Sentinel-2 (13.3×103 m²·yr⁻²), and lower than Landsat (14.6×103 m²·yr⁻²). At Xirinongpu, the fused mean variance is 16×103 m²·yr⁻², below Landsat (18.7×103 m²·yr⁻²) and Sentinel-2 (19.9×103 m²·yr⁻²), and also better than GoLIVE (37.3×103 m²·yr⁻²). Moreover, relative to the two pre-fusion optical products, the fused series shows markedly reduced variance fluctuations, indicating greater smoothness and stability. The lower variance of ITS_LIVE compared with the fused product because ITS_LIVE is used here as an annual-scale product, whereas the fused results are monthly velocities. The longer temporal aggregation of ITS_LIVE smooths short-term fluctuations and therefore tends to yield lower apparent variance.

 “

Fig 6. Same what’s the unit of the axis ?

Modified accordingly.

1. 345 could you precise the null hypothesis that you have tested here ? is it the one for scipy ? h ps://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.pearsonr.html

Thank you for this helpful comment. Yes, the significance test reported here follows the default hypothesis test implemented in scipy.stats.pearsonr. Specifically, the null hypothesis is that the two variables are uncorrelated (i.e., the population Pearson correlation coefficient is r=0), and we used the default two-sided test.

Fig 12 : here there is no way of deciding if the trend is significant or not. For that, you need to perform a null hypothesis test, as done in Halas et al :

https://doi.org/10.1016/j.rse.2022.113419

In the section just after you are providing a p value but it’s not clear : what test did you used ?

Thank you for this important comment. We would like to clarify that Figure 12 does not represent a significance-tested trend map. Instead, it shows the direction and magnitude of the pixel-wise linear trends in glacier surface velocity during 2015–2024. The statistical significance test and the reported p values apply only to the glacier-level trend analysis presented in Figure 13, where linear regression was performed for each glacier (> 1 km²) and the resulting trends were classified as significantly accelerating, significantly decelerating, or non-significant.

We agree that the original wording could lead to confusion, especially where we used the term “significant” in the discussion of Figure 12. In the revised manuscript, we have therefore clarified this distinction explicitly and replaced the relevant wording in the Figure 12 discussion so that it refers only to the spatial pattern and magnitude of trends, not to statistical significance. The explanation of the p-value test has also been made explicit in the Figure 13 section.

Fig 14 : why is it written slope for the y-axis and not decadal trend in velocity ?

Modified accordingly.

In the conclusion  « Pearson correlations indicate significant relationships between velocity and glacier area, slope, and aspect: » is a strong statement recording the pearson value you got.

Thank you for this helpful comment. We agree that the original wording in the Conclusion was too strong given the Pearson correlation coefficients obtained in this study. In the revised manuscript, we have softened this statement and now refer more cautiously to statistically significant correlations/associations between glacier surface velocity and glacier area, slope, and aspect, rather than implying a stronger relationship than is supported by the results.

“In spatial terms, glacier surface velocity during 2015–2024 exhibits the canonical “fast center, slow margins” pattern: multi-year mean maxima exceed 700 m·yr⁻¹, whereas values in lower reaches and most tributaries are generally <100 m·yr⁻¹. Regarding attribute controls, Pearson correlation analysis indicates statistically significant but moderate associations of velocity with glacier area, slope, and aspect: larger glaciers tend to flow faster; within individual glaciers, velocity is more strongly correlated with slope; and after accounting for area and slope, south-facing glaciers are slightly faster than glaciers of other aspects.”