青海贵德地区潜水含水层弥散系数研究

摘要 以青海贵德地区作为研究区域，以NaCl作为示踪剂，一次性瞬时投放，进行野外弥散试验，观测示踪剂浓度随时间的变化过程，使用标准曲线配线法和使用Matlab数据拟合法进行求解，得到了贵德地区潜水含水层弥散系数，并分析了试验结果，比较了两种计算方法。

关键词 弥散试验;弥散系数;标准曲线配线法;Matlab数据拟合

中图分类号 S181.3;P641.89 文献标识码 A 文章编号 0517-6611（2014）33-11819-03

Study of Unconfined Aquifer Dispersion Coefficient in Qinghai Guide Area

PI Kai-hong， KANG Wei-dong， LUO Qi-bin et al

（Department of Geology/State Key Laboratory of Continental Dynamics， Northwest University， Xi’an， Shaanxi 710069）

Abstract With Qinghai Guide as an object of study areas， and NaCl as a tracer， a one-time instantaneous delivery was conducted， field diffusion experiments were carried out， the change of tracer concentration with time was observed， using standard curve fitting method and the Matlab data fitting method to resolve， unconfined aquifer diffusion coefficient of Guide area was obtained， the experiment results were analyzed and the two calculation methods were compared.

Key words Diffusion experiments; Diffusion coefficients; Standard curve fitting method; Matlab data fitting

作者简介 皮锴鸿（1990- ），男，陕西眉县人，硕士研究生，研究方向：水文地质。

收稿日期 2014-10-10

地下水弥散系数是用来表征污染物在多孔介质中运移规律的重要参数，一般具体地区含水层的弥散度和弥散系数通过野外弥散试验来获得[1]。野外试验的流场一般在天然流场和天然人工叠加流场中进行，其中天然人工叠加流场使用较多，其试验过程为在投源孔注入示踪剂，在各个观测孔观测其浓度在时间、空间上的变化，再依据客观条件使用相应的计算方法求解参数。对于野外弥散试验弥散参数求解的方法有多种[2-3]，包括直线图解法[4]、改进直线解析法[5]、弥散晕面积法、标准曲线配线法[6]、逐点求参法等。

贵德县位于青海省东部，处于黄河上游龙羊峡与李家峡之间，全境沟壑纵横，山川相间，地势南北高、中间低，且呈现多级河流阶地，形成四山环抱的河谷盆地。境内有多条河流由南北两侧汇入黄河，水资源较为丰富，再加上近些年各种引水工程的迅速建设，致使工农业迅速发展，进而有可能引发污染问题。该研究即针对贵德地区地下水污染物溶质运移弥散系数通过野外试验进行具体研究。

1 弥散试验

试验场地位于河谷阶地之上，为河流堆积相，岩性为灰黄色粉土夹杂砂砾石，上部为腐殖质粉土层，地层结构松散，砂砾石粒径一般细，砾石成分以砂岩为主，次为花岗岩及脉石英，分选磨圆较好，均为次圆状。区域地层孔隙率约为0.2，钻孔资料显示含水层厚约30 m。

钻孔布置及钻孔间距如图1所示。在抽水状态下，自然流场相对于抽水时人工形成的流场非常小，可以忽略，故可认为地下水流方向沿X轴方向。

图1 钻孔布置示意

该试验选取氯化钠（NaCl）为示踪剂，20 ℃时NaCl溶解度为36 g[7]。用NaCl作示踪剂对全区域环境基本上无损害，非常适合用于小范围的简单连通示踪试验。在进行示踪试验时，直接用化学方法检测示踪剂浓度，比较麻烦而且不经济;通常测量含NaCl地下水的导电率，再利用NaCl溶液电导率与其浓度之间关系，间接换算得到NaCl浓度变化曲线。

这次试验抽水孔抽水速率为50 m3/h，同时观测地下水位，待水位稳定之后进行弥散试验，将175 kg NaCl溶解后，一次性瞬时投放入投源孔，并开始计时，观测g2孔和CH孔浓度变化。整个试验持续时间根据NaCl浓度随时间变化曲线具体确定[8]，当NaCl浓度从最低的背景值到达峰值，再降低到最低，并且之后持续一段时间NaCl浓度不再变化，终止试验。

这次弥散试验历时55 h，NaCl浓度随时间的变化过程见图2。经测定，试验场地地下水电导率背景值为0.504～0.606 mS/cm，即NaCl浓度背景值为170.33～204.80 mg/L。由图2可以看出试验数据具有相似的趋势：g2观测孔中NaCl浓度开始基本无变化，到23：00开始升高，在26.83 h NaCl浓度达到峰值，8 516.34 mg/L，之后逐渐下降到背景浓度附近;CH监测孔NaCl浓度，从试验开始到第52 h时基本无变化，之后开始升高，在52.83 h后抽水孔中NaCl浓度达到峰值831.36 mg/L，之后逐渐下降到背景浓度附近。其中在浓度接近背景值试验时间段，观测浓度在背景值上下游波动，其原因在于取样及测试的随机性等所引起[9]。

图2 示踪监测孔NaCl浓度随时间的变化曲线

2 标准曲线配线法

由于抽水孔附近天然流场中地下水流速与人工抽水所形成的流场中地下水流速相比可以忽略不计，可认为形成以抽水井为中心的径向流场，同时径向地下水流速很大，则忽略横向弥散作用不会导致较大误差，故径向渗流场中溶质运移的对流-弥散方程可以在柱坐标系下简化为：

ct=aL|u|2cr2-ucr

法国水文地质学家Sauty[10]采用有限差分的数值法进行参数计算求得以Peelet数p为参数，以无因次浓度Cr为纵坐标和无因次时间tr为横坐标的标准曲线，用以确定含水层的纵向弥散度。

依据现场试验实测数据，按照下式将示踪剂观测浓度C换算成无因次浓度Cr，试验观测时间t换算成无因次时间tr。

Cr=C-C0Cmax-C0

tr=tt0

式中，C0为示踪剂的背景浓度（mg/L），Cmax为示踪剂的峰值浓度（mg/L），t0为纯对流时间（h）。

纯对流时间（t0）计算公式为：

t0=π（r22-r21）mmQ

式中，r2为投源孔至抽水孔的距离（m），r1为溶质浓度检测孔至抽水孔的距离（m）。

对试验数据进行无因次浓度Cr、无因次时间tr计算处理后，绘制与标准曲线同模的Cr-tr实测曲线，并与标准曲线相匹配，找到实测曲线与某一P值的标准曲线匹配最好的情况，确定P值，再由公式αL=r/P计算得到纵向弥散度（径向弥散度）αL值。

该试验观测孔与投源孔间距r为5.06 m，纯对流时间t0为27.9 h，抽水状态下计算得到地下水流速为0.052 4 m/h，通过配线得到P值为800（图3），则纵向弥散度αL为0.006 325 m，纵向弥散系数为DL为0.000 331 4 m2/h。抽水孔与投源孔间距9.98 m，配线得到P值为1 000（图4），则纵向弥散度αL为0.009 98 m，纵向弥散系数DL为0.000 523 m2/h。

图3 g2孔标准曲线配线结果

图4 CH孔标准曲线配线结果

3 数学软件数据拟合法

对试验观测数据直接拟合求解参数是一种相对比较直观、容易理解的求参方法，但由于所涉及的未知参数较多，且方程比较复杂的时候求解起来就非常困难，所以常借助数学软件求解，常用拟合求解参数的软件有SPSS、SAS、EXCEL、Matlab、Mathematica、Maple等，而Matlab相比其他软件具有语言简单、功能强大、扩充能力强、移植和开放性好等特点，同时提供多种数学工具箱供用户使用。该试验使用Matlab编写程序对纵向弥散参数DL求解[11]。

首先，将对流-弥散方程变形简化。当观测孔位于x轴上时（即y=0），同时考虑示踪剂在含水层中的反应速率[2]，解的形式为：

c（x，t）=m4πntDLDTe[-（x-vt）24DLt-λt]（1）

根据试验场地下游前人弥散试验成果以及经验判断，场地横向弥散度大致为纵向弥散度的1/10，得到：

c（x，t）=m4πntDL0.1e[-（x-vt）24DLt-λt]（2）

式中，c（x，t）为观测浓度，可根据弥散试验观测得到，含有未知参数纵向弥散系数DL，地下水流速v，反应速率λ。

使用Matlab编写程序，将（2）式设置为拟合函数模型，时间t定义为自由变量，观测浓度c为非自由变量：

ft = fittype（... " 175000/（4*pi*0.2*30*dl*sqrt（0.1）*t）*exp（-（5.06-u*t）^2/（4*dl*t）-r*t） " ， " independent " ， " t " ， " dependent " ， " y " ）;

%设置模型

opts = fitoptions（ft）;

opts.Display = " Off " ;

%设置模型运行选项

opts.StartPoint = [ 0.2 0.1 0.1];

%设置拟合参数的起始点，DL的起始点设置为0.2，λ 的起始点为0.1，u的起始点为0.1。程序的起始参数是模型运行的起点，如果设置的不合适可能无法拟合出结果，所以在计算的时候需要多次尝试给出不同的值。

[fitresult，gof]=fit（xData，yData，ft，opts）

%拟合数据x，y列数据，得到拟合结果和好的拟合结果统计量。

将弥散试验观测结果带入程序中，得到对监测孔g2拟合结果（图5）：纵向弥散系数为0.000 452 1 m2/h，反应速率系数为0.158 8，地下水流速为0.185 4 m/h;抽水孔CH拟合结果（图6）：纵向弥散系数为3.778E-005 m2/h，反应速率系数0.167 9，地下水流速为0.189 m/h。

图5 监测孔g2试验数据与拟合曲线

图6 抽水孔CH试验数据与拟合曲线

4 结论

（1）通过对两种不同方法计算得到的结果进行比较，得到青海贵德地区纵向弥散系数约为0.000 331 4～0.000 452 1 m2/h之间，横向弥散系数约为0.000 033 1～0.000 045 2 m2/h之间。其中通过两种方法计算得到抽水孔CH的值相差比较大，故不作为最终结果考虑。

（2）两种计算方法各有优缺，配准曲线法需要通过其他试验求得地下水流速，但求解过程相对容易操作;通过软件拟合求参可将地下水流速设为参量直接拟合求得，但对于参数的起始点较难设定，需要多次尝试给予赋值，才能得到理想的拟合结果。

参考文献

[1] 谢荣焕，刘满.基于MATLAB的水质模型参数的确定方法[J].工业安全与环保，2006（2）：28-29.

[2] 郭建青，马致远.分析二维水动力弥散试验数据的线性回归法[J].西北水资源与水工程，1999（1）：14-18.

[3] 吴飞，曹玉清，李苍松.二维水动力弥散的随机数学模型[J].数学的实践与认识，2004（6）：106-111.

[4] 吴耀国，田春声，李云峰.判定地下水二维水动力弥散参数的直线图解法[J].煤田地质与勘探，1997（2）：31-33.

[5] 郭建青，周宏飞，王洪胜.分析二维水动力弥散试验数据的改进直线解析法[J].勘察科学技术，2011（1）：14-17.

[6] 丁佳.银川地区潜水含水层弥散参数试验确定方法研究[D].西安：长安大学，2010.

[7] 郑克勋.地下水人工化学连通示踪理论及试验方法研究[D].南京：河海大学，2007.

[8]丁家平，李樟苏.地下水弥散系数的野外试验新方法[J].水利学报，1998（8）：38-42.

[9] 蒋辉.弥散参数计算方法分析[J].勘察科学技术，2013（3）：9-12.

[10] SAUTY J P.An analysis to hydrodispersive transfer in aquifers[J].Water Resources Research，1980，16（1）：145-158.

[11] 赵勇胜.弥散度及其在地下水污染模型中的作用[J].长春地质学院学报，1992（4）：208-211.