青海貴德地區潛水含水層彌散系數研究

皮鍇鴻 康衛東 羅奇斌 蘇賀 唐歡

摘要 以青海貴德地區作為研究區域，以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， Xian， 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.