冰川均衡调整(glacial isostatic adjustment，GIA)效应是黏弹性地球对冰负荷变化的动力学响应，会引起地壳运动与地幔流动，可为研究地幔对流、板块构造和板块俯冲等提供基本的约束条件^[1]。地球表面经历过多次冰期，其冰进—冰退过程自然能够引发GIA效应。GIA效应引起的全球响应变化使得量化GIA信号成为地球科学研究中的一个关键问题，已有许多学者在此领域进行过研究与讨论。关于量化GIA信号的方法大致可以分为两类：一是基于数值模拟的方法构建GIA模型，如ICE系列模型、Paulson07模型、Caron系列模型等；二是借助现代大地测量技术使用实际观测数据对GIA信号进行量化^[2-15]。其中，重力恢复与气候实验(gravity recovery and climate experiment，GRACE)卫星能够实现全球时变重力场连续监测，在现代GIA信号量化研究中应用广泛，特别是在北美、北欧、南极等区域。比如，Steffen等^[16]基于GRACE重力数据得到斯堪的纳维亚的重力变化趋势，发现与前人对该地区的GIA研究结果非常一致。Tamisiea等^[17]利用GRACE重力数据确定一个横向均质地球模型，可以限制GIA效应对劳伦蒂亚静态重力异常场的贡献。van der Wal等^[18]基于GRACE和水文模型量化哈德逊湾的GIA信号。高春春等^[19]基于GRACE和遥感卫星数据，利用全球定位系统(global positioning system，GPS)球谐拟合改正等方法分离了南极的GIA信号。张腾宇和金双根^[20]利用GRACE、GPS和雷达卫星估计了青藏高原的冰后回弹效应。

在南美洲的巴塔哥尼亚高原上存在南半球最大的极地外冰原——巴塔哥尼亚冰原(Patagonia Icefield，PIF)。研究表明近些年来PIF内的冰川融化速率较高且在逐年加快^[21-23]，这种大型温带冰块快速变化引起的GIA效应不可忽视。但目前巴塔哥尼亚地区的GIA响应研究很少，其GIA响应信号需进一步明确，因此本文研究目标是利用GRACE卫星时变重力数据量化该区域的GIA响应信号。

基于GRACE卫星研究GIA效应本质上是信号分离的过程，从GRACE探测信号中扣除水文信号可以得到GIA响应变化，其中水文信号是估计重力变化率的重要组成部分^[18]。Jiao等^[24]为了降低水文信号的不确定性，使用多种水文模型的均值作为青藏高原的水文变化估计值，结合GRACE卫星重力信号分析青藏高原的莫霍面变化。Shen等^[25]基于GRACE卫星研究中国陆地非水文时变重力的空间分布特征及机理时，利用多源水文模型数据精细构建了一个水文模型，用以准确扣除水文信号。Sun等^[26]在研究智利三联点附近的南美大陆莫霍面变化特征时，对比不同水文模型的雪水变化趋势，发现有较大差异。这说明，现有水文模型之间存在一定的差异，使用单一水文模型不确定性较大。为减少不确定性而获取准确的水文信号，本文综合不同水文数据完善了南美巴塔哥尼亚的水文模型。

本研究利用GRACE卫星重力数据与水文数据量化南美巴塔哥尼亚地区的GIA响应，主要研究目标为：1)完善该高原地区的水文模型，了解陆地水文信号的变化特征；2)从GRACE卫星综合质量变化信号中扣除水文信号得到该高原地区GIA响应信号，分析其变化特征；3)基于GPS数据了解该地区地表抬升情况，讨论GIA效应对PIF区域地表抬升的贡献。

1 研究区域概况

巴塔哥尼亚高原位于南美洲板块、纳兹卡板块及南极洲板块的交界处——智利三联点以东，覆盖南美洲板块的最南端，有一部分位于斯克提亚微板块之上(图 1)。高原西起安第斯山脉，东临大西洋，地形从西向东呈阶梯状倾斜。高原西部受到西风带和太平洋湿空气团输送的影响，降水充沛^[27]。由于智利海隆扩张中心与南美大陆碰撞，纳兹卡板块和南极洲板块相对于南美洲板块的汇合速度存在差异，巴塔哥尼亚高原下方发育了一个很大的板片窗^[28]，其边界参考Russo等^[29]。PIF分为北巴塔哥尼亚冰原(Northern Patagonia Icefield，NPI)和南巴塔哥尼亚冰原(Southern Patagonia Icefield，SPI)两部分。其中NPI面积约4 200 km²，体积约1 200 km³；SPI面积约13 000 km²，体积约4 300 km^3[22]。GPS观测PIF存在异常隆升现象，其垂直速度超过4 cm/a。Bourgois等^[30]根据数值模拟方法分析该异常隆升过程与内部构造、外部气候因素密切相关。

2 数据与方法 2.1 重力数据与正向建模法

本文采用美国德克萨斯大学空间研究中心发布的2003年1月至2016年12月的RL06版本的60阶月地球重力场模型(http://icgem.gfz-potsdam.de/series)，进行低阶项补充与替换、高斯滤波、去相关滤波等处理。其中高斯滤波半径300 km，去相关滤波选择P4M6方式^[31]。得到以等效水高形式表示的质量变化

$\begin{gathered} \Delta H(\theta, \lambda)=a \frac{\rho_{\text {ave }}}{3 \rho_{\mathrm{w}}} \sum\limits_{l=0}^{\infty} \sum\limits_{m=0}^l \bar{P}_{l m}(\cos \theta) \frac{2 l+1}{1+k_l} \times \\ \left(\Delta C_{l m} \cos (m \lambda)+\Delta S_{l m} \sin (m \lambda)\right), \end{gathered}$

(1)

其中：a是地球半径；l、m分别表示阶和次数；θ、λ分别代表余纬和经度；ρ_ave是地球平均密度，约5 517 kg/m³；ρ_w是水的密度，为1 000 kg/m³；P_lm(cosθ)是完全正规化的勒让德函数；C_lm、S_lm是无量纲的球谐系数。

再用最小二乘法拟合质量变化长期趋势^[26]

$\begin{aligned} \Delta g= & A+B \Delta t+\sum\limits_{i=1}^3 C_i \cos \left(\omega_i \Delta t\right)+ \\ & D_i \sin \left(\omega_i \Delta t\right)+\varepsilon, \end{aligned}$

(2)

其中：t是时间；A是常数；B是年变化速率；C_i、D_i是振幅，ω_i是相位，i=1表示年周期变化，i=2表示半年周期变化，i=3是与S2潮汐波相关项；ε是残差。

在GRACE数据处理中，截断和滤波会导致信号泄漏问题^{[24, 26]}。本文采用正向建模方法对GRACE观测信号进行恢复^[32]。该方法可以用于校正GRACE在数据处理中引起的泄漏误差。主要步骤如下：首先将研究区域的初始质量变化m₀作为迭代方法的输入信号，记为f₀；其次将m₀扩展为球谐系数，截断到60阶并进行滤波，得到正向模型表观质量率f₁；然后将f₀与f₁的差值添加到初始信号m₀中，比例因子设置为1.2，用于加快迭代的收敛速度，即m₀=m₀+1.2×(f₀-f₁)。当f₀和f₁一致时，或者在一定次数的迭代之后，停止迭代，m₀即为恢复后的真实变化信号。

2.2 CryoSat-2卫星的冰川数据结果

巴塔哥尼亚高原的冰川变化信息主要体现为PIF的质量变化。本文根据Foresta等^[22]基于CryoSat-2条带雷达卫星观测得到的PIF质量损失结果，通过以下公式得到PIF各区域的质量变化ΔH(等效水高形式)

$\Delta H=\frac{\Delta M}{S},$

(3)

其中：ΔM是质量损失(质量单位)；S是对应的面积。1°×1°格网中存在多个冰川区域对应一个格网的现象，选择格网中冰川质量变化最大值赋予该格网，得到PIF冰川质量变化趋势(PIF区域外的格网值为0)，球谐系数展开至60阶，再滤波平滑处理。

2.3 水文模型 2.3.1 GLDAS水文模型

全球陆地数据同化系统(global land data assimilation system，GLDAS)输出1979年至今的陆地表面水文数据，但该模型没有模拟深层地下水以及地表湖泊水库的变化影响，本研究中可用的数据仅土壤水、冠层水和雪水3种成分。GLDAS水文模型包含不同子模型，如NOAH、CLSM、VIC模型。经对比(图 2)，不同子模型之间存在一定的差异，说明模型存在较大的不确定性，但目前还没有更好的数据能够应用于研究中。因此，为最大程度地减少单个模型带来的偏差进而获取到更为可靠的水文成分变化，本文选择不同子模型的平均值作为水文成分的最终估计值，空间分辨率为1°×1°，时间分辨率为1个月。为了与GRACE卫星数据共同计算，对GLDAS水文数据进行球谐截断、滤波平滑处理。

2.3.2 WaterGAP水文模型

WaterGAP全球水文模型量化了1901—2016年全球陆地水资源，不仅模拟了人类对地下水和地表水的利用，还单独模拟了每个水文成分，空间分辨率为0.5°×0.5°，时间分辨率为1个月。这里对2003年1月至2016年12月WaterGAP水文模型的湖泊、地下水、湿地、水库与河流数据重采样至1°×1°后，进行与GLDAS水文模型相同的处理。

2.4 GIA信号的提取

由GRACE重力卫星观测得到的质量变化(ΔGRACE)反映地球系统质量迁移的综合信号，在GIA近场区域，主要包含水文与GIA信息，存在以下关系：

$\Delta \mathrm{GIA}=\Delta \mathrm{GRACE}-\Delta \text { Hydrological },$

(4)

$\begin{gathered} \Delta \text { Hydrological }=\Delta \text { Glacier }+\Delta \text { Snow }+\Delta \text { Soil }+ \\ \Delta \text { Lake }+\Delta \text { Groundwater }+\Delta \text { Canopy }+ \\ \Delta \text { Wetland }+\Delta \text { Reservoir }+\Delta \text { River }, \end{gathered}$

(5)

式(4)中ΔHydrological指水文质量综合变化，ΔGIA指由GIA效应引起的质量变化；式(5)表示考虑的水文信息有冰川、雪水、土壤水、湖泊、地下水、冠层水、湿地、水库及河流，对应的质量变化由ΔGlacier、ΔSnow、ΔSoil、ΔLake、ΔGroundwater、ΔCanopy、ΔWetland、ΔReservoir及ΔRiver分别表示。

2.5 GPS数据

本文使用共计75个GPS站点的垂直速度数据，这些数据来自内华达大地测量实验室和前人研究^{[26, 33-34]}，用以分析GIA响应在PIF地表隆升中的贡献率。我们选择的GPS站点的工作时间主要集中在2003—2016年；分布范围主要集中在PIF区域。

为得到巴塔哥尼亚高原的垂直速度场，要选用合适的插值方法。这里选择克里金插值法，适用于分布不均匀的数据插值。Jiao等^[24]提出GPS站点垂直速度与等效水高质量变化存在以下关系

$\Delta H_{\text {water }}=\frac{\Delta h_{\mathrm{GPS}} \rho_{\mathrm{crust}}}{\rho_{\text {water }}},$

(6)

其中：ΔH_water是1°×1°网格质量变化对应的等效水高值；Δh_GPS是网格中GPS站点数据的平均垂直位移；ρ_crust是网格中的地壳密度，计算方法是根据CRUST 1.0全球地壳模型中的8个层的密度取加权平均值，权重是层的厚度；ρ_water=1 g/cm³。

3 结果 3.1 GRACE观测的巴塔哥尼亚高原综合质量变化特征

GRACE重力卫星观测到的巴塔哥尼亚高原质量变化趋势的空间分布如图 3所示。其中，图 3(a)是GRACE卫星观测到的巴塔哥尼亚高原的表观质量率，图 3(b)是通过正向建模方法恢复的真实质量率，图 3(c)是对真实质量率经截断、滤波计算的模拟表观质量率，图 3(d)是图 3(a)和3(c)的差异。

由图 3(b)可知，巴塔哥尼亚高原的质量变化趋势具有明显的空间差异性，呈现出高原中部质量亏损、高原南部和北部质量增加的空间分布特征。其中质量亏损最显著的区域位于SPI最北端((73°W，49°S)，格网点B)，变化趋势约为(-2.23±0.12) cm/a。从该区域向外质量变化速度逐渐增加，直至在高原西北缘((72°W，43°S)，格网点A)表现出质量增加的最大趋势，速度约为(0.78±0.09) cm/a。该区域冰体快速融化与质量显著亏损密切相关。

3.2 巴塔哥尼亚高原水文模型的空间分布特征

本研究主要参考Foresta等^[22]的PIF质量变化结果与GLDAS、WaterGAP水文模型数据，完善巴塔哥尼亚高原的水文模型。考虑2003—2016年巴塔哥尼亚高原内冰川、雪水、土壤水、湖泊、地下水、冠层水、湿地水、水库和河流的质量变化信息，提取陆地水文信息特征。

3.2.1 冰川变化的空间分布

本文根据Foresta等^[22]计算的PIF冰川质量变化结果，利用式(3)得到PIF各区域的冰川质量变化，如图 4(a)所示。可以看出，冰川质量变化的信号中心位于PIF区域，最大变化趋势约为(-2.83±0.42) cm/a，越远离信号变化中心，变化幅度越趋近于0。

3.2.2 土壤水、冠层水、雪水变化的空间分布

土壤水、冠层水、雪水的质量变化来自GLDAS的3种子模型的均值，如图 4(b)~4(d)所示。其中土壤水整体呈现质量减少的趋势，东北部亏损最严重；冠层水变化微弱；雪水亏损在PIF附近最为显著，其他区域亏损趋于零值。

3.2.3 地下水、湖泊等其他水文成分变化的空间分布

基于GLDAS水文模型获取土壤水、冠层水和雪水3种水文成分数据，还需要其他水文成分来构建巴塔哥尼亚高原的水文模型。这里选择WaterGAP水文模型，得到湖泊、地下水、湿地水、水库和河流的质量变化，其空间分布如图 4(e)~4(i)所示。其中湖泊、地下水和水库的质量变化相对显著，三者的质量亏损主要区域分别位于PIF及其附近区域、高原靠近太平洋的西海岸区域、高原最北端；而湿地和河流的质量变化趋势较弱。

3.2.4 巴塔哥尼亚高原水文模型的空间分布

本文联合不同水文数据完善巴塔哥尼亚高原的水文模型(图 4(j))，该模型反映了巴塔哥尼亚高原总的水文信息质量变化。其变化范围在(-3.36±0.51)~(-0.13±0.05) cm/a，整体呈现出质量亏损的状态，高原内最显著的负信号出现在PIF区域，形成以冰原为中心向外辐射、负信号逐渐减弱的空间分布特征。Ndehedehe和Ferreira^[35]用测高卫星数据计算2002—2014年南美洲的水文质量变化，其中巴塔哥尼亚高原的水文质量年变化约为-3.5~0.5 cm/a，该结果与笔者得到的结果基本一致，证明了本研究水文信息特征提取的可靠性。

3.3 巴塔哥尼亚高原GIA响应变化的空间分布特征

由式(4)可知，GRACE重力卫星探测信号中主要包含水文信息与GIA信息。将水文信息从GRACE重力卫星探测结果中扣除，就得到巴塔哥尼亚高原的GIA信息，再通过正向建模法恢复真实信号，如图 5所示。结果表明，GIA响应与高原冰川分布情况一致，南部高于北部、西部高于东部，在高原的北缘和SPI南部存在2个信号中心，引起的质量变化值如表 1所示。变化范围在0~2 cm/a，最高达(1.97±0.35) cm/a，使高原整体隆升。利用最小二乘法计算GRACE和水文信号的变化趋势与标准差，将不确定性估计为2个标准差，再根据误差传递公式计算得到GIA信号的不确定性，置信水平为95%。

GIA信号量化有2种方法——数值模拟或大地测量。为了解这2种GIA信号方法在南美巴塔哥尼亚高原的GIA信号差异、验证本文GIA结果的可靠性，将通过大地测量手段得到的结果与数值模拟方法给出的GIA模型进行对比分析。

首先选择近10年内发布的4种成熟的全球GIA模型，分别是Geruo13模型、ICE-6G_D模型、Caron18模型和Caron19模型，在图 6(a)~6(d)和表 2中分别给出4种模型的空间分布情况及其与本文结果的相关系数。整体上GIA模型与本文结果均表示巴塔哥尼亚高原的GIA响应使地表隆升；相关系数介于0.52~0.82，具有良好的相关性，说明其空间分布特征相似。

其次搜集2种以PIF为基准建立的区域GIA模型，分别为模型A与模型B，其冰川历史和地幔黏度参数如表 3所示。根据Lange等^[33]提供的模型A、B隆升速率数据，利用式(6)计算出区域GIA模型A、B表示的质量变化，其空间分布如图 6(e)~6(f)所示。结合图 5(b)，在PIF区域模型A、B最大变化值分别约为1.81和2.04 cm/a。这2种区域GIA模型冰川历史与本研究时间尺度接近，GIA信号特征相似，证明了本研究结果的可靠性。

3.4 GPS观测的巴塔哥尼亚高原地表隆升情况

通过克里金插值与式(6)得到巴塔哥尼亚高原地表垂直速度引起的质量变化空间分布，如图 7(a)所示，变化最显著的区域位于SPI北端。考虑到GIA信号是根据球谐系数截断、滤波等处理后的GRACE信号与水文信号综合计算得到的，我们对GPS结果进行了相同的处理，如图 7(b)所示。可以看出在巴塔哥尼亚高原中部、PIF及其附近区域，GPS观测的垂直速度引起的质量变化最为显著，最高达(2.79±0.36) cm/a。

4 讨论 4.1 GIA信号与全球GIA模型存在差异的原因

GIA效应是3个不同时间阶段的综合结果^[36]：第1个阶段是冰后回弹，主要由约2.1万年前的大冰川消退引起，该时段也称为末次冰盛期；第2个阶段是小冰期后退，发生在几百年以前；第3个阶段是现在正发生的冰川融化。由于大地测量技术应用时间的局限性，本研究得到的GIA信息主要是黏弹性地球对第3阶段现今冰川融化的响应信号，研究时间段为2003—2016年。而全球GIA模型是通过不同冰川历史与负荷的冰模型和不同地球物理参数的地球模型综合计算得到的^[1]，主要包含末次冰盛期以来的GIA响应。通过数值模拟得到的GIA模型与本文通过大地测量方法计算的GIA结果均为线性年变化率。随着冰川融化速率逐年加快，GIA响应会逐渐增强，因而出现本文短时间尺度内计算的GIA结果高于全球GIA模型的现象。

4.2 GIA效应在PIF地表抬升中的贡献率

巴塔哥尼亚高原位于大型温带冰块快速变化和特殊构造的地区，其垂直变形可以通过GIA效应解释。为准确了解GIA在PIF地表抬升中的贡献率，计算了PIF区域GIA响应(图 5(b))与GPS垂直速度(图 7(b))的比值，将NPI及其附近划分为区域N，SPI周围划分为区域S(如图 7(a)所示)。在N区域和S区域，GIA响应与GPS垂直速度的比值均值分别约为0.692 5和0.827 0，表明在这2个区域GPS探测信号中分别约有69.25%和82.70%的垂直信号可以通过GIA效应解释。N、S区域GIA信号对地壳隆升的贡献存在约10%的差异，主要原因可能在于NPI和SPI的气候不同导致冰川卸载速率不同。SPI与NPI中间被一条宽阔的非冰川河流分离，并不相连接。SPI的年降水量约为NPI的4倍^[37]，年冰川损失质量约为NPI的3~5倍^[28]，气候与冰川卸载速率存在显著差异，导致GIA响应对地壳隆升的贡献也不同。

在PIF区域的地壳隆升主要由GIA效应贡献，但考虑PIF区域的构造环境，除GIA效应以外，板片窗打开引起的地形动态变化与气候因素也有可能对地壳隆升造成影响。一方面，PIF恰好位于板片窗之上。Georgieva等^[28]认为在巴塔哥尼亚高原沿大西洋海岸阶地的隆起和倾斜的河流梯田反映了大陆岩石圈对板片窗打开的动态响应；另一方面，PIF位于南北向的安第斯山脉顶部，是天然的屏障，山脉两侧气候差异非常大。Hu等^[38]应用三维、全动态俯冲模型数值模拟安第斯山脉的生长机制，认为纳兹卡板块、南极洲板块向南美板块的俯冲汇聚时，气候因素使得地表沉积物减少，板块耦合进一步增强，安第斯山脉隆升。这说明气候与冰川侵蚀速率变化也有可能导致PIF山峰高度增加^[28]。

因此，引起PIF区域地壳隆升最主要的因素即GIA效应，贡献约70%~80%；此外还受到板片窗打开、气候等因素的综合影响，对地壳隆升的贡献率在20%~30%左右。

4.3 地壳构造运动的影响

地壳构造运动是引起地表隆升或沉降的长期过程。GRACE重力卫星能够探测地球表面与地球深部的综合质量变化，要想从中精确提取GIA响应信号，应当进一步考虑地壳构造运动的影响。但Klemann等^[39]研究GIA与板块运动关系时发现，GIA效应引起的地幔流动是地壳构造运动的主要驱动力之一。这说明二者相互耦合。此外，在时间尺度上，Klemann等^[39]提出，地壳构造运动的特征时间尺度为1~100 Ma，远大于本文研究时间尺度(14 a)。Thomson等^[40]在研究巴塔哥尼亚的造山运动时，发现该地区平均构造抬升速率约为0.07~0.2 mm/a。这表明，在本文研究时间尺度内，地壳构造运动引起的隆升与质量变化非常小。因此本研究在计算过程中忽略了地壳构造运动的影响。

5 结论

南美巴塔哥尼亚高原位于冰雪覆盖的构造地区，GIA响应不可忽视。本文综合GRACE卫星重力数据和水文数据，优化了对南美巴塔哥尼亚高原GIA响应信号的估算。主要结论如下：

1) 通过完善巴塔哥尼亚高原的水文模型，发现2003—2016年该高原地区水文质量整体亏损，形成以PIF为中心向外辐射、负信号逐渐减弱的空间分布，变化范围在(-3.36±0.51)~(-0.13±0.05) cm/a。这种亏损分布特征与冰体融化密切相关。

2) 从GRACE综合信号中剔除水文信息，认为剩余信号主要是GIA响应，使巴塔哥尼亚高原地表隆升，最高达(1.97±0.35) cm/a。位于PIF区域的GIA信号变化特征与区域GIA模型高度一致；与全球GIA模型相似且有良好相关性，但数值上有较大差异，这与GIA模型的构建原理和时间尺度有关。

3) 基于GPS观测数据，GIA效应在NPI和SPI区域分别能够解释约69.25%和82.70%的地表抬升信号，其余影响因素有气候变化和板片窗的打开等。

4) 为估计巴塔哥尼亚高原的GIA响应信号，完善了巴塔哥尼亚高原的水文模型，虽然该模型包含几乎所有的陆地水文信号，但仍存在精度不高、不确定性较大的问题。此外，GRACE卫星数据的处理过程与GPS站点分布、数据插值导致的精度误差，都会影响最终的研究结果。这些问题导致的不确定性有待进一步探讨与改进。

感谢ICGEM(http://icgem.gfz-potsdam.de/series)提供的GRACE重力卫星数据；感谢NASA(https://disc.gsfc.nasa.gov/datasets)提供的GLDAS水文模型数据；感谢Müller Schmied等学者提供的WaterGAP水文模型数据(http://www.watergap.de/ListOfPubl.html)；感谢内华达大地测量实验室提供的GPS数据(http://geodesy.unr.edu/index.php)；感谢Peltier、Geruo、Caron、Lange等学者分享的GIA模型数据。本文图件主要使用ArcGIS软件绘制。