一维水分运移模型水力参数敏感性分析

刘明明，段 磊，张 琛，宋 浩，吕婷婷，霍世璐

(1.长安大学环境科学与工程学院，西安 710054；2.旱区地下水文与生态效应教育部重点实验室，西安 710054)

近年来，随着计算的发展，利用数值模拟手段，研究包气带水分运移已经日益引起了人们的重视，并逐渐成为研究包气带水分运移规律的重要手段[1-3]，但是，由于模型参数的不确定性影响了模型精确度，给建立包气带一维水分运移模型的取值，造成了很大的困难,所以为了提高模型精度，针对参数的不确定性，进行定量的分析与研究有着重大意义。

在数值模拟中，分析和测定模型参数对模型输出结果的影响，这就是参数的敏感性分析[4,5]。目前，针对参数的不确定性对包气带水分运移模型的研究比较多。王志涛等[6]利用HYDRUS-1D，采用单因素扰动法，得出n、θs和k的扰动对垂直湿润距离计算值和累积入渗量的影响最大；毕经纬[7]、孟立江[8]以底部通量作为敏感性分析的指标，得出饱和水力传导度是最敏感参数；张添等[9]利用路面模型，得出浅层土壤水分对土壤质地较为敏感，敏感性系数达到了0.45以上 ，并且砂土含量对土壤水分的影响更为显著 ；但是，针对参数不确定性对非均质剖面含水率空间分布的研究比较少。因此，本文基于HYDRUS[10]模型利用局部敏感性分析的方法，对模型参数进行敏感性分析，旨在找出影响剖面含水率空间分布的显著参数，给借助模型研究土壤水分动态提供参考。

1 实验区概况

实验区位于新疆霍城县伊犁河流域中下游的细土平原区，该地区属温带半干旱气候，平均气温8.2～9.4 ℃，年平均蒸发量约850～1 000 mm，年平均降水量约140～460 mm。

为观测灌水之后剖面含水率的变化，利用含水率观测仪器TDR，分别对剖面在0、10、20、40、60、90、140、190、200和240 cm深度的含水率进行监测，监测频率为次/3 h。

2 数值模拟

2.1 概念模型

根据实验场概况，建立一维水分运移模型，模拟裸土试验场在灌水之后，水分在剖面上的空间分布。层状非均匀包气带岩性如表1所示，其中主要岩性为粉质黏土，并含有中砂以及粉土夹层(图1)。因此，基于包气带岩层的考虑，本文将模型概化为四层。上边界为开放边界，接受降水、灌溉水补给及土壤蒸发排泄，蒸发量根据中国气象网提供的气象数据，利用Peman-Monteith计算得到[11,14]。潜水埋深为4.5 m，大于潜水面毛细上升高度，因此，本文将下边界设置为自由排边水边界。

表1 土样颗粒分析实验表Tab.1 Soil sample particle analysis experiment table

图1 岩性剖面图Fig.1 Lithology profile

2.2 数学模型

(1)

式中：h为土壤负压，cm；c(h)为容水度，1/cm；h0为初始土壤负压值，cm；P为降雨量，cm/d；I为灌溉量，cm/d；E为蒸发量，cm/d；t为时间，d；B为模型下边界；

HYDRUS采用van-Genuchten模型描述土壤水分特征曲线和土壤导水率曲线两类参数[15]，模型表达公式如下：

(2)

(3)

(4)

式中：θ(h)为土壤的特征曲线；θ为土壤含水率，%；θr为土壤的残余含水率，%；θs为土壤的饱和含水率，%；h为压力水头；K(h)为压力水头下的导水率，%；ks为饱和导水率，%；α、n和m为van-Genuchten模型的形状参数；Se为有效饱和度。

2.3 模型的离散化

采用等间隔的方法，以1 cm为间隔将包气带剖分为241个单元格，模拟期从2016年10月1日-2017年9月31日，从中分析参数对剖面含水率空间分布的影响，模型采用变时间间隔的方法，初始时间间隔为0.001 d，最小时间间隔为0.001 d，最大时间间隔为2 d。

2.4 土壤参数的识别

利用筛选法将用环刀在实验场取的原状土样进行定名，再利用TST-70饱和渗透仪，测出饱和导水率，最后利用德国生产的ku-pF非饱和导水率测定系统[16]，得到负压和含水率的关系，再利用MATLAB最小二乘法对实验数据进行参数拟合[17-19]，得到土壤水力参数如表2。

基于残差平方和均方根误差进行误差分析，并从表3的结果可知，土壤水分特征曲线拟合比较好，符合要求，因此本文将表2作为模型的初始值。

表2 土壤水分特征曲线拟合参数以及实测参数Tab.2 Soil moisture characteristic curve fitting parameters and measured parameters

表3 土壤水分特征曲线拟合误差表Tab.3 Soil moisture characteristic curve fitting error table

2.5 模型识别与验证

模拟期内经过两次灌水试验，因此本文分别将两次灌水后的第8 d剖面含水率作为识别与验证期，在经过对比研究后，发现拟合效果良好(表4、图2)。

表4 识别后土壤水力参数Tab.4 Soil hydraulic parameters after identification

图2 土壤剖面含水率拟合图Fig.2

3 参数敏感性分析

3.1 参数敏感性分析理论依据

局部分析的方法，具有高效以及操作简单的优点。因此，本文将利用局部分析的方法对参数空间中的一点进行敏感性分析[20]。

本文使用相对敏感性值将参数进行归一化处理，计算出敏感指数 ，并根据敏感数进行分类[21,22]。

I=(ΔO/ΔFi) (Fi/O)

(5)

式中：O为模型输出结果；Fi为影响O的因子(参数)；ΔO为模型输出结果的改变值；ΔFi为影响O的因子(参数)的改变量。敏感性分类见表5。

表5 敏感性分类Tab.5 Sensitivity classification

3.2 参数取值

根据上述岩层残余含水率与室外实测初始含水率的对比发现，当残余含水率增加的值超过20%时，残余含水率将高于岩层最低初始含水率，因此为了便于比较参数对剖面含水率的敏感性，本文将水力参数上下变幅统一变为20%，以此来计算水力参数的敏感指数，表6为水力参数取值范围。

表6 水力参数的取值Tab.6 Value of hydraulic parameters

3.3 土壤含水量模型参数的敏感性分析

本文利用参数局部敏感性分析的手段，通过改变水力参数的单一变量，来研究灌水15 cm之后，水力参数对剖面含水率影响的敏感指数，结果如表7所示。

表7 裸地情况下影响土壤含水量模型水力参数的敏感性分析Tab.7 Sensitivity analysis of hydraulic parameters affecting soil water content model under bare soil conditions

本文认为当敏感指数绝对值大于0.05时，即可认为输出结果对参数变化敏感(分级按照水力参数对剖面含水率影响的最大敏感指数进行归类)。

从表7参数的敏感指数可知，粉质黏土层θs1上浮20%时，敏感指数属于Ⅲ级，最大敏感指数最大达到0.60，对土壤上部粉质黏土剖面含水率，影响比较大，但是对下部其他岩层的影响比较小，当θs1下浮20%时，敏感指数属于Ⅲ级，最大敏感指数达到0.91，但影响范围只局限于本岩层。中砂岩层θr2上下浮动20%时，敏感指数几乎相等，同时达到Ⅲ级，最大敏感指数达到0.98，对本层的剖面含水率影响最大，对其他岩层含水率并没有造成太大的影响。粉土层 上浮和下浮20%时，敏感指数差不多，最大达到0.72，同属Ⅲ级，对第三层剖面含水率影响比较大，当n3上下浮动时，发现当n3向下浮动时，其敏感指数达到了Ⅳ级，说明了n3是对第三层剖面含水率影响最大的参数。当θs4上浮和下浮20%时，敏感指数都达到Ⅲ级，但n4对剖面含水率影响指数达到了Ⅳ级，表明了在下部粉质黏土层中，n4对该层剖面含水率影响最大，而且对上部粉土层也产生了很大的影响。

为了验证上述敏感性分析结果，本次模拟将灌水变为20 cm，来研究水力参数对剖面含水率的影响(表8)。

表8 裸地情况下影响土壤含水量模型水力参数的敏感性分析Tab.8 Sensitivity analysis of hydraulic parameters affecting soil water content model under bare soil conditions

从表8参数的敏感指数可知，上部粉质黏土n1和θs1对剖面含水率的影响最大，不过n1对该层剖面含水率最大，而对于中砂层来说，θr2对该层的剖面含水率，起着决定性作用，而其他水力参数，对剖面含水率几乎没有造成影响，而粉土和下部粉质黏土层的n和θs共同影响着剖面含水率，但是n对该层剖面含水率最大，由此不难看出，灌水15 cm和灌水20 cm时，得出的结论一致。

3.4 分析与讨论

通过对表7和表8的对比发现，在上部粉质黏土层中，n1和θs1对第一层剖面含水率敏感指数最高，并且敏感指数随着深度的增加而减少。在中砂层中，只有θr2对本层剖面含水率敏感性最大，其他水力参数上下浮动，难以对剖面含水率产生很大的影响，粉土岩层以及下部粉质黏土层中n和θs对本层的含水率，有着重大的影响，由此可以看出，在细颗粒中n和θs对剖面含水率影响比较大，而且n的敏感指数要大于θs的敏感指数，这主要是VG模型参数中θs对水分特征曲线低吸引力段影响显著，而n变化会影响土壤水分特征曲线的弯曲程度以及整体形状，但其参数主要在土壤释水的后期，发挥着主导性作用[23,24]，因而在灌水后期，细颗粒岩层n对剖面含水率的影响程度大于θs。VG模型中θs在灌水后期开始发挥着重要作用[25]，而中砂层由于颗粒比较粗，所以土壤水在前期释水的过程中，θs2和n2发挥着重要作用，但是随着释水过程的推进，延迟与滞缓的时间，对剖面含水率影响逐渐减弱，当剖面含水率降低接近该岩层θr2时，θr2对该剖面含水率影响占据了主导作用。

4 结 语

水分经过300多天的水分运移之后， 在非均质岩层细颗粒以及粗颗粒岩层中，α和k参数的变化难以对剖面含水率产生太大的影响，而在细颗粒岩层中，θs和n对剖面含水率发挥着重要着作用，在粗颗粒岩层中，θr对剖面含水率的影响最大。

数值模型在实际应用的过程中，由于岩层复杂多变，给模型参数的取值造成了很大困难，因此本次研究结论，可以减少在非均质中构建模型时，需要调整的参数，对于提高模型精确度有着重大意义，所以，在进行数值模拟时，应对细颗粒岩层的θs进行加密测量，而对n要进行重点调参，在粗颗粒的岩层中，要提高对θr的重视程度。