1. 引言
地下水是非常重要的自然资源之一，地下水位的年、季变化将引起灌溉、饮水和种植结构的变化。区域地下水埋深是了解区域地下水资源状况的一个重要参数，它与区域生态环境息息相关。而潜流带是河流河床内水分饱和的沉积物层，是河流和地下水存在物质和能量交换相互作用的区域 [1] ，占据着水文循环过程中的重要位置。研究潜流交换过程及其特点对准确评价和开发利用流域水资源，维持和修复地表水及地下水生态系统健康具有重要意义。
作为研究潜流通量方法的一种，温度示踪法对于研究潜流通量具有其独特的优势，相对于液体比重法和地下水溶质示踪剂法，温度示踪法更加方便快捷；相对于化学示踪法，温度示踪法不会对水环境造成污染且易于观测。由于温度数据示踪潜流交换受河床空隙多少、大小、形状、方向性、联通程度及其空间变化的影响较小 [2] ，所以具有非常广泛的适用性。Peter C. Young等人于1999年，在动态谐波回归模型的基础上，提出了处理非静态时间序列温度数据的更简便的模型，使得热运移方程解析解的计算变得更加容易 [3] 。在2006年，Hatch等人提出了根据温度曲线相位及振幅的变化计算河床渗透率的方法 [4] ，并开发了适用于饱和多孔介质中一维垂向流的程序VFLUX [5] 。同时，John Keery等人在2006年也给出了已有热运移方程数值解的扩展解，该方法可以在较小的实验面积内计算潜流带不同深度处潜流交换的垂直流速 [6] 。Mc Callum等人和Luce等人分别在2012年与2013年对热运移方程进行了近一步的解算和数值分析，使得VFLUX模型不仅可以解算潜流通量，也可以解算热混合率、最优传感器间距等 [7] [8] 。尽管热示踪法在潜流通量的研究上得到了极大的发展，但是受到主观认知能力和客观高度复杂性的限制，河流潜流带潜流量计算的结果仍然充满了不确定性。温度示踪无法克服温度数据估计潜流交换的部分缺点，比如：河床的空间异质性以及河床一些热学参数的选取对潜流通量计算结果影响比较大。为了提高计算结果的准确性，更精准地说明各热力学参数对潜流通量的敏感程度，反映河流与地下水的交换机理，进行地下水潜流通量的敏感性分析变得尤为重要。
本研究旨在填补基于温度示踪法的河流潜流带潜流量计算参数的敏感性这方面的空白，将采取渭河陕西段实测数据作为原始数据，通过研究热弥散度、热传导能力、水的比热容、河床体积热容等热力学参数以及孔隙度对结果的敏感性来更精确地分析进而对变化情况和规律进行模拟，提高潜流通量计算结果的准确性，为当地的环境治理和经济开发作参考。
2. 研究区概况
如图1所示，本次试验的研究区位于渭河流域陕西省西咸新区渭河咸阳西安段。黄河流域中上游地区处于我国西北干旱地区，降水少，自然生态系统脆弱，渭河位于中国西部干旱半干旱地区，是黄河最
Figure 1. The map of Weihe River and the location of the test site
图1. 渭河水系和现场试验场地示意图
大的一级支流，发源于甘肃省渭源县鸟鼠山，流经陇、宁、陕三省，于陕西省潼关注入黄河，全长818 km，总流域面积1.35 × 104 km2。渭河陕西段全长502 km，流域面积6.76 × 104 km2，流域南北跨越黄土高原、关中盆地、秦岭北麓山地3个地貌单元，气候类型属典型的大陆性季风气候，冬季寒冷干燥，夏季炎热多雨，降水多集中于6~9月份且多暴雨，加之渭河流域大部分被黄土覆盖，质地疏松，易被水蚀致使渭河成为典型的冲淤性河流。
渭河下游的咸阳–草滩段位于西安凹陷的东北部，河道呈直线状的游荡型，东西长约40 km，东抵临潼隆起，南有沣河、灞河等注入，北有泾河相汇，所处地貌单元主要为关中平原，南与秦岭山前洪积扇及黄土台塬相邻，北与渭北台塬相接，南北宽50~80 km，平均海拔360~450 m。研究区的渭河河道宽浅，枯水期河宽50~100 m，洪水期河宽100~200 m。河心沙洲发育，主槽位置摆动频繁。河漫滩沿渭河两岸分布，宽度为2~6 km，低漫滩高出河床1~3 m，高漫滩高出河床5~6 m。受两岸河堤约束绝大部分河漫滩已不再受洪水淹没。草滩断面海拔高度360 m左右。
渭河及支流河漫滩与各级阶地广泛分布河谷平原潜水含水岩组。含水层岩性以全新统、上更新统和中更新统冲积砂、砂砾卵石为主，其次为亚粘土由西到东岩相变细，厚度变薄。河漫滩的含水层厚度一般为10~30 m；渭河一级阶地为45~60 m；三级阶地为4~8 m。潜水位埋深随地形升高而增大，一般为1~40 m [9] 。
3. 材料与方法
3.1. 热示踪法计算潜流通量原理
潜流带的温度具有典型的日变化和季节性特征，可以作为反映潜流交换的有效信息。理论上，垂向潜流交换通量的计算可以利用任意两组不同深度的温度监测数据实现 [10] 。一维条件下，考虑热传导，对流和扩散作用的热量运移方程如公式(1)所示：
(1)
公式中：
是温度(℃)，
是时间(s)，
是有效热弥散系数(m2/s)，
是距离河床表面的深度(m)，
是垂向潜流交换通量(m/s)，
是饱和介质热容与水热容之比。
Hatch等 [11] 给出了公式(1)的解析，可以分别采用振幅比和滞后时间两种方法来计算不同位置的潜流交换通量。河水温度可以作为埋深为0处的河床温度 [12] 。垂向潜流交换通量(
)分别采用振幅比法(
)和相位差法(
)可以表示为：
(2)
(3)
公式(2)和(3)中：
为河床内不同深度处两组温度信号的振幅之比(无量纲)，
为两个深度处温度信号的最大值或最小值出现的时间差(s)，
和
分别为介质骨架与水的热容(J/M3/℃)，
为河床中两测点的高程差(m)，
是热峰的运移速度(m/s)，
是温度信号的周期(s)，
是定义如下的参数：
(4)
饱和介质的热容按照下式进行计算：
(5)
公式(5)中：
为饱和介质的孔隙度(无量纲)。
采用Vukovic和Soro提出的经验公式进行估计 [13] [14] ：
(6)
公式(6)中：
是介质颗粒级配的不均匀系数 [15] 。公式(2)，(3)和(4)中的热弥散系数
采用公式(7)进行估计：
(7)
公式(7)中：
为有效热传导系数(J/s/m/℃)，
为热弥散度(m)。潜流通量计算值代表每两个温度监测位置之间的等效通量值，这里采用每对温度测点的中点位置来代表该潜流通量的空间位置。更多的计算方法细节可以参见相关文献 [16] 。
本次计算采用Gordon等 [17] 开发的VFLUX工具箱进行潜流通量计算，由于相位法和振幅法均可表示结果，但相比之下振幅比法更适宜作为小通量条件下(q < 1000 L/m2/d)的潜流通量计算 [18] 。鉴于本次计算结果多数小于1000 L/m2/d，因此这里仅对振幅比法的结果进行分析。我们采用的具体程序为MATLAB 7.10.0 (R2010a)的VFLUX2.0.0版，如果信号处理工具不可用，可前往http://www.es.lancs.ac.uk/cres/captain地址进行下载。
3.2. 敏感性参数及赋值范围
热传导性可以量化介质对水流热量传递的阻力，但热传导性随河床的变化却非常小。河床周围的介质粒径大小的变化会直接影响潜流带渗透性的变化，但热传导性与沉积层的结构无关，热传导性可以量化介质对水流热量传递的阻力，但热传导性随河床的变化却非常小。David通过试验发现粒径大小的变化对潜流带的热力学参数几乎不影响或只有很小的影响 [19] [20] ，见图2，所以本文需要分析的敏感性参数即为上述VFLUX计算所采用的热力学参数，即为孔隙度(
)、河床体积热容(
)、水热容(
)、热弥散度(
)、有效热传导系数(
)。
参数调试幅度也是要解决的难题之一，参数调试幅度太大，则可能引起模型失真，从而得出错误的敏感性分析结论。若参数调试幅度太小，则可能看不出变化，或者有变化，但可能是由模型本身计算误差所引起；参数值浮动范围的选取应根据实际情况自行设定，尽量减小误差 [21] [22] [23] [24] [25] 。因为VFLUX程序经过了前期的重采样、数据过滤等程序，计算误差往往处于一个很低的数量级，远远小于敏感性参数赋值变化所引起的结果变化，所以计算误差可忽略。
根据前人经验和文献查阅 [20] ，针对实际研究区，确定本次研究各参数的典型值，上边界值和下边界值，如表1所示。
Table 1. Estimated thermal properties of the saturated streambed used for 1-D modeling of vertical hyporheic flux
表1. 垂向潜流通量计算中采用的热力学参数
Figure 2. Sensitivity of permeability coefficient and thermal conductivity to dielectric size [4]
图2. 渗透系数与热传导率对介质粒径的敏感性 [4]
3.3. 野外实验采集数据
为收集数据，进行野外原位实验。经现场勘探，暂定横断面具体地点为渭河西段左岸。同时，由于上下游流速不同，本实验需沿河岸在不同的横断面上设置温度记录仪。温度记录仪集测温、记录存储和数据传送为一体。可自动记录和存储传感器测得的温度，连接计算机即可将温度数据传送与计算机处理，得到原始数据。
温度传感器使用具有记录功能的温度记录仪iBUTTONDS1922L-F5，呈纽扣形状，直径17.35 mm，厚6 mm，分辨率采用0.0625℃，可以高精度实时记录温度细微变化，并通过专用读卡器与计算机连接，导出时间和温度系列实测温度数据集。
实验具体措施为将温度记录仪和电脑连接，对记录仪进行参数设定，采样速率设置为10分钟，利用冲击钻将温度记录仪钻孔，放入自封袋中，并将自封袋口利用硅胶进行密封，用铁丝把自封袋固定于木棒上。设计每根木棒上安置1个至5个温度记录仪。将1号温度传感器设置为河床下深度0.01 m，用来监测河水的温度，其它4个传感器的间隔依次为0.05，0.15，0.40，0.75 m。同时将对气温、水温和河床介质温度的监测频率均设置为10分钟；利用橡胶锤将固定了温度传感器的木棒打入河床，连续记录河床浅层沉积物的温度变化。
本次温度监测时段为2015年9月30日至10月3日，测得地表水温度在13℃~20℃之间变化，而地下水温度约稳定在21℃左右。
3.4. 敏感性研究方法
目前的敏感性分析中默认各参数间是相互独立的，本文采用单个控制变量法进行敏感性分析，即保持其他参数不变，只改动一个参数的取值，在深度为0.1 m，0.275 m，0.575 m的图形中寻找不同，分析不同深度同一参数的敏感性差异和同一深度不同参数的敏感性强弱，并分析其中原因，同时将数据进行汇总和整理，计算潜流通量受到各参数影响的变化率。
4. 结果与讨论
4.1. 同参数不同取值比较结果
根据参数的不同取值，通过VFULX程序得到了不同深度的垂向潜流通量结果(单位为m/s) (如表2)。
根据表2和表7，可以看出在0.05~0.70 m的深度时，孔隙度越大，垂向潜流通量越大，孔隙度越大，说明河床介质的渗透性越好，所以交换量越大。在0.15~0.40 m的深度孔隙度由0.28到0.35潜流通量变化最明显，说明孔隙度的敏感性在此阶段最高。在0.40~0.70 m的深度潜流通量随孔隙度变化最均匀，没有明显的急剧变化。
根据表3和表7，可以看出在0.05~0.70 m的深度时，热弥散如何变化，潜流通量都不会随之改变，说明潜流通量的大小在取值范围在0~0.1之间不受热弥散度影响，因为热弥散不论如何变化，热传感器接收到的温度信号的时间会变，但接受热信号反映的温度的大小不会发生改变，也决定了由热示踪法所得到的潜流带交换量不会发生改变。所以在热示踪法求潜流通量时，热弥散的敏感性几乎为零。
根据表4和表7，可以看出在0.05~0.15 m的深度时，热传导系数越大，潜流通量越大。在0.15~0.40 m的深度时，潜流通量不随热传导系数变化而保持不变。在0.40~0.70 m的深度，热传导系数越大，潜流通量越小。在0.05~0.15 m和0.40m~0.70 m的深度时，热传导系数在0.0025~0.003时变化量比0.002~0.0025变化量大。热传导系数随深度不同对潜流通量的敏感性差异比较大，较为复杂。
根据表5和表7，可以看出在0.05~0.70 m的深度时，河床体积热容的值越大，潜流通量的值越大，在0.05~0.15 m的深度时，河床体积热容在0.5~0.6时潜流通量的变化量是河床体积热容在0.4~0.5时潜流通量的变化值的三倍。在0.15~0.40 m的深度时，潜流通量随河床体积热容的变化均匀变化。在0.40~0.70 m的深度时，河床体积热容在0.4~0.5时潜流通量的变化量比河床体积热容在0.5~0.6时潜流通量的变化值大。河床体积热容在0.4~0.5时，0.15~0.70 m深度的潜流通量变化量大于0.05~0.15 m深度的潜流通量变化量。河床体积热容在0.5~0.6时，0.05~0.15 m深度的潜流通量变化量大于0.15~0.40 m深度的潜流通量变化量且大于0.40~0.70 m深度的潜流通量变化量，即深度越深敏感性越差。
根据表6和表7，可以看出在0.05~0.70 m的深度时，水体积热容越大，潜流通量越小。在0.05~0.15 m深度时，水体积热容在0.9~1之间潜流通量变化量大于1~1.1之间潜流通量变化量，在0.15~0.70 m深度均匀变化。
通过上表可得，孔隙度在0.15~0.40 m敏感性比0.40~0.70 m强，热弥散无敏感性，热传导系数的敏感性强弱：0.40~0.70 m > 0.05~0.15 m > 0.15~0.40 m；河床体积热容敏感性强弱：0.05~0.15 m > 0.40~0.70 m > 0.15~0.40 m；水体积热容敏感性强弱：0.05~0.15 m > 0.15~0.40 m > 0.40~0.70 m。
4.2. 不同参数比较结果
根据表8可以得知，0.05~0.15 m深度不同参数敏感性强弱：热传导系数 > 河床体积热容 > 水体积热容 > 孔隙度 > 热弥散；0.15~0.40 m深度不同参数敏感性强弱：孔隙度 > 河床体积热容 > 水体积热容 > 热传导系数 = 热弥散 = 0；0.40~0.70 m深度不同参数敏感性强弱：热传导系数 > 河床体积热容 > 孔隙度 > 水体积热容 > 热弥散。平均各变化率在0.05~0.70 m深度不同参数敏感性强弱：热传导系数 > 河床体积热容 > 孔隙度 > 水体积热容 > 热弥散。
表2. 孔隙度敏感性数据统计表
Table 3. Thermal diffusion sensitivity data table
表3. 热弥散敏感性数据统计表
Table 4. Thermal conductivity coefficient sensitivity data table
表4. 热传导系数敏感性数据统计表
Table 5. Bed volume heat capacity sensitivity data table
表5. 河床体积热容敏感性数据统计表
Table 6. Water volume heat capacity sensitivity data table
表6. 水体积热容敏感性数据统计表
Table 7. Different parameters of different values of different depth of the hyporheic flux change scale
表7. 各参数不同取值不同深度潜流通量变化量表
Table 8. Different parameters of different depth of the hyporheic flux rate of change table
表8. 不同参数不同深度潜流通量变化率表
5. 结论
根据以上，可得到以下结论：
1) 通过VFLUX软件在利用热示踪法求解潜流通量时，各参数取值的不同会影响潜流通量的大小：在0.05~0.70 ms深度，潜流通量随孔隙度、河床体积热容增大而增大，随水体积热容增大而减小。热弥散在0~0.01取值时不影响潜流通量结果。
2) 与其他参数不同，热传导系数对潜流通量的影响随深度不同而发生变化：深度在0.05~0.15 m时，潜流通量随热传导系数增大而增大。深度在0.15~0.40 m时，潜流通量随热传导系数增大而不变。深度在0.40~0.70 m时，潜流通量随热传导系数增大而变小。
3) 同参数在不同深度的敏感性强弱也有所不同：水体积热容的敏感性随深度增大而变弱。孔隙度在0.15~0.40 m时敏感性最强，热传导系数在0.40~0.70 m时敏感性最强，河床体积热容在0.05~0.15 m时敏感性最强。
4) 不同参数在同深度的敏感性强弱：0.05~0.15 m深度不同参数敏感性强弱：热传导系数 > 河床体积热容 > 水体积热容 > 孔隙度 > 热弥散；0.15~0.40 m深度不同参数敏感性强弱：孔隙度 > 河床体积热容 > 水体积热容 > 热传导系数 = 热弥散 = 0；0.40~0.70 m深度不同参数敏感性强弱：热传导系数 > 河床体积热容 > 孔隙度 > 水体积热容 > 热弥散。平均各变化率在0.05~0.70 m深度不同参数敏感性强弱：热传导系数 > 河床体积热容 > 孔隙度 > 水体积热容 > 热弥散。
基金项目
本课题为长安大学大学生创新训练计划(201610710225)和宁波市科技富民项目(2015c10043)资助。