秦就國，蔣亞萍

(桂林理工大學，廣西桂林541106)

0 引言

傳統(tǒng)的地下水數(shù)值模擬程序主要局限性在于:(1)傳統(tǒng)地下水數(shù)值模擬程序所應(yīng)用到達西定律，不能刻畫巖溶含水層中普遍存在紊流;(2)傳統(tǒng)的地下水數(shù)值模擬程序處理層流所用的是等價連續(xù)模型(equivalent continuum models(ECM))，將巖溶地區(qū)多重孔隙、各向異性進行等價均一處理。這些直接導(dǎo)致在管道處的水頭，觀測值與模擬值存在明顯的差異［1］。

為了克服傳統(tǒng)地下水數(shù)值模擬程序的這些局限性，考慮到MODFLOW的廣泛應(yīng)用與模塊化設(shè)計結(jié)構(gòu)［2］，美國地質(zhì)勘察局(The U.S.Geological Survey(USGS))基于 MODFLOW －2005(version 1.2.01)開發(fā)管道流數(shù)值模擬程序 CFP(Conduit Flow Process)，用于處理巖溶地區(qū)非達西流、多重孔隙度的復(fù)雜情況。

1 CFP簡介

1.1 CFP 適用范圍

根據(jù)含水層的構(gòu)造情況，CFP設(shè)計成三種運行模式，CFP Model1(CFPM1)、CFPModel2(CFPM2)、CFPModel3(CFPM3)。

CFPM1模式能夠刻畫石灰?guī)r地區(qū)中石灰?guī)r溶解和生物洞穴特征、玄武巖中的熔巖管特征，含水層飽和、半飽和，水流是層流與紊流的條件下都適用。

CFPM2中，加入一個強透水能力的水流層(high－conductivity flow layer)，使得模型可以在層流與紊流之間進行轉(zhuǎn)化。此模式適用3種情形:①多孔介質(zhì)中，在所觀測到的水力坡度下，有可能出現(xiàn)紊流的地方。②一個簡單的次生孔隙度地層特征，如界限明確的橫向地下溶洞。③由大量相互連通的孔隙所組成的水平優(yōu)先流層。

CFPM3是對CFPM1、CFPM2的耦合，在同時包括CFPM1、CFPM2所適用的情形下適用。

1.2 CFP 原理簡介

考慮以下3種含水層構(gòu)造情形:

1.2.1 含水層包含巖溶管道(CFPM1)

1)模型幾何拓撲結(jié)構(gòu)

CFPM1將層流模型(Modflow－2005)與離散的管道網(wǎng)絡(luò)耦合成雙孔隙度模型。管道網(wǎng)絡(luò)由大量管道組成，每個管道由兩個端點組成，端點稱為管道結(jié)點(node)。管道結(jié)點位于MODFLOW有

式中:Δh，hL為沿著管道長度的水頭損失，單位［L］;f為摩擦因子，無量綱;d為管道直徑，單位［L］;V為平均水流速度，單位［L/T］;g為重力加速度，單位［L/T^2］;單位體積流量Q=V·A，其中A是管道水流橫截面積。

所以限差分單元中，一個差分單元最多只有一個管道結(jié)點，一個結(jié)點至少連接1根管道，最多連接6根管道。管道可以在相鄰的差分單元間相連，也可以連接對角的差分單元。管道結(jié)點在水平方向上位于MODFLOW差分單元中心，但在垂向上，管道結(jié)點可以①高于差分單元中心點;②低于差分單元中心點;③剛好位于差分單元中心點。

每個管道結(jié)點都被分配唯一的結(jié)點編號，并由用戶指定高程，每2個管道結(jié)點間的管道也被分配唯一的管道編號。用戶在CFP文件中指出管道結(jié)點在MODFLOW三維模型中的位置i，j，k(行，列，層);并賦予各個管道的屬性:管道直徑d，彎曲度τ、粗糙度kc、管道導(dǎo)水系數(shù)或管道壁滲透系數(shù)等。

2)管道水流量的計算

管道完全充滿水時，Darcy，Weisbach，還有其它水力學專家提出管道水流公式(1)，此公式同時適用層流和紊流:

摩擦因子f是雷諾數(shù)Re的函數(shù)，與管道的粗糙程度相關(guān)。雷諾數(shù)Re是液體密度ρ、動力粘滯系數(shù)μ的函數(shù)，ρ、μ都又做為地下水溫度Tw的函數(shù)，通過地下水平均溫度即可計算出Q。

研究表明，用水力直徑de而不是管道直徑d，在層流與紊流的情形下，計算結(jié)果更準確。根據(jù)水力直徑定義:

其中:A—過水斷面積，P—被水浸潤的管道周長。

假設(shè)管道是水平的，D是液壓深度，即管道內(nèi)水平面到管道底部的深度。根據(jù)幾何學可求出:

其中:

注意:管道中水位高于管道中心點時，用θ2=2π－θ替換式(4)、式(5)中 θ。

圖1 有效直徑計算

3)管道中儲水量變化量計算

在未充滿水的管道內(nèi)，當管道內(nèi)水面上升或下降時，考慮管道內(nèi)儲水量Qs的變化。

管道內(nèi)的水體積 Vt=AΔlτ，故

其中θt可根據(jù)式(5)計算得出。

4)管道與MODFLOW差分單元水量交換的計算

考慮管道與多孔介質(zhì)間存在水量交換，CFP采用線性公式計算MODFLOW差分單元與管道間的交換水量Qex:

其中:

Qex［L3T－1］為單位體積交換量;

αj，i，k在 MODFLOW 單元 j，i，k 處的管道水力傳導(dǎo)系數(shù);

hn位于MODFLOW單元中管道結(jié)點n的水頭;

hj，i，k包圍著 MODFLOW 單元 j，i，k 的水頭。1.2.2 含水層不包含巖溶管道(CFPM2)

一些石灰?guī)r含水層不包含巖溶管道，但具有高孔隙度，可看成是強透水層，地下水流呈紊流狀態(tài)。CFPM2基于以下兩部分實現(xiàn):

1)判斷水平水流是紊流還是層流

CFP用雷諾數(shù)作為界定水流是層流還是紊流的臨界值。由于運行慣性，層流趨向保持層流狀態(tài)，紊流趨向保持紊流狀態(tài)。用戶必須為CFP指定兩個臨界雷諾數(shù)NRe。例如雷諾數(shù)上下限分別為3 000、2 000，若當前水流由層流向紊流過渡，2 000≦Re≦3 000時，水流仍然是層流，＞3 000時，水流轉(zhuǎn)變?yōu)槲闪?同樣的，若當前水流由紊流向?qū)恿鬟^渡，2 000≦Re≦3 000時，水流仍然保持為紊流狀態(tài)，Re＜2 000時，水流變成層流。

2)調(diào)整被MODFLOW－2005所假設(shè)的水力傳導(dǎo)系數(shù)Klam

由于MODFLOW是基于層流假設(shè)，CFPM2模式下，CFP將紊流狀態(tài)下的水力傳導(dǎo)系數(shù)Kturb轉(zhuǎn)換成層流狀態(tài)的水力傳導(dǎo)系數(shù)Klam。引進調(diào)整因子Fadj，無量綱，用來調(diào)整層流水力傳導(dǎo)系數(shù)Klam:

其中:

Δh為水平方向兩個相鄰單元的水頭差;

Δhcrit為臨界水頭差，可根據(jù)臨界雷諾數(shù)求出。

1.2.3 同時包含 A、B(CFPM3)

此情形下稱為CPFM3，即是將情形A、B進行簡單的耦合，不再贅述。

2 管道流建模方法

目前暫時還沒有軟件完全支持對CFP進行圖形化輸入、后處理。對于CFP建模，其參數(shù)輸入主要還依賴于手工完成，用文本編輯器直接對輸入文件以一定的格式進行編輯，指定CFP所需要的各項參數(shù)。

文件采用自由格式，即程序以行為單位進行讀取，程序?qū)?shù)據(jù)的具體位置并不敏感，各個數(shù)據(jù)以空格或TAB隔開，以“#”為首的行是注釋信息，程序不讀取此行。

2.1 管道參數(shù)輸入模塊

不同的模式，CFP輸入文件也會不同:

2.1.1 CFPM1

CFP輸入文件指定當前所處的運行模式(模式1)，描述管道的空間位置結(jié)構(gòu)，CFP文件指出各個管道結(jié)點在MODFLOW差分單元格中的位置及高程，管道的各項參數(shù)(彎曲度、粗糙度、管道滲透系數(shù)、管道直徑)，連續(xù)介質(zhì)與管道間水量交換系數(shù)，平均水溫;用于求解的其它參數(shù)(迭代次數(shù)、松弛因子等);用于層流與紊流轉(zhuǎn)換的臨界雷諾數(shù)。

2.1.2 CFPM2

相對于CFPM1，CFPM2不需要模擬管道流，輸入?yún)?shù)要簡單得多。CFP輸入文件指定當前所處的運行模式(模式2)，描述管道層的總數(shù)，水溫，每一層的平均孔隙直徑、層流與紊流相互轉(zhuǎn)換的兩個臨界雷諾數(shù)。

2.1.3 CFPM3

CFPM3是把CFPM1、CFPM2進行耦合起來，同時包含CFPM1、CFPM2的情況。因此，CFPM3的輸入文件最復(fù)雜，同時包含CFPM1、CFPM2輸入文件所需要的所有參數(shù)。

2.2 管道輸出控制模塊

大部分CFP計算結(jié)果都被寫入到MODFLOW－2005列表文件(List File)中，沒有寫到列表文件中的結(jié)果被保存到CFP的COC文件(Conduit Output Control File)中，以便進行錯誤檢查和后處理。COC文件用于控制管道結(jié)點輸出信息，如輸出頻率、管道個數(shù)等。注意此文件在CFPM1與CFPM3才會被用到。

2.3 管道補給模塊

Conduit Recharge(CRCH)模塊用來處理將面補給進入巖溶管道這一過程，比如降水直接通過落水洞進入巖溶。在CFPM1和CFPM3中，CRCH模塊可以讓部分或所有的補給水進入到某個指定的管道結(jié)點。注意CRCH是RCH模塊的一個擴展，在CFPM1、CFPM3中，Modflow－2005的RCH模塊必須啟用。

3 CFP性能評價

3.1 準確性測試

3.1.1 室內(nèi)基準測試

實驗由一根管道、包圍著管道的連續(xù)介質(zhì)組成［3］，水力坡度在到的范圍內(nèi)調(diào)整以模擬層流，結(jié)果與Hagen－Poiseuille公式計算結(jié)果進行比較，擬合度達到100%;水力坡度由1.00～－1.01的范圍內(nèi)進行調(diào)整以模擬紊流，結(jié)果與Darcy－Weisbach公式計算結(jié)果進行比較，擬合度大于99.86%［3］。室內(nèi)測試表明，在巖溶地區(qū)含水層，MODFLOW－2005 CFP模擬結(jié)果比MODFLOW－2005或 MODFLOW －2000 都要精確［4］。

3.1.2 次區(qū)域研究

Davisetal于2010年，基于 MODFLOW －2000運用 EPM(E-

圖7 不同農(nóng)藝措施處理土壤含鹽量變化

2.2.7 棉花產(chǎn)量統(tǒng)計分析

2012年6－9月，試驗區(qū)所在區(qū)域降水較常年偏多48%。特別是在2012年7月31日至8月1日，試驗區(qū)的降雨量達157 mm，造成大面積內(nèi)澇，對試驗小區(qū)種植的棉花產(chǎn)量造成了一定影響。

棉花產(chǎn)量按照各試驗小區(qū)的實際產(chǎn)量進行統(tǒng)計，棉花產(chǎn)量在不同處理的兩個重復(fù)之間的變化趨勢基本一致，故取其平均值。棉花產(chǎn)量統(tǒng)計結(jié)果顯示，P2處理產(chǎn)量最高，兩個重復(fù)小區(qū)平均收獲籽棉6.3 kg;其次是P3處理，兩個重復(fù)小區(qū)平均收獲籽棉6 kg;對照處理平均收獲棉花5.4 kg。P1、GL1、GL2三個處理，平均收獲棉花基本都在5.7 kg，各處理間棉花產(chǎn)量沒有明顯差異;GL3處理的棉花產(chǎn)量最低為5.3 kg。2號地各農(nóng)藝措施試驗處理的棉花產(chǎn)量統(tǒng)計分析，見圖8。

圖8 不同農(nóng)藝措施處理棉花產(chǎn)量統(tǒng)計

3 結(jié)語

(1)土壤有機改良劑和作物秸稈翻耕等農(nóng)藝措施，可以明顯改善鹽堿土壤的理化性狀，提高鹽堿土壤中有機質(zhì)及其他養(yǎng)分含量，對于黃河三角洲低產(chǎn)鹽堿土壤改良，提高其土地質(zhì)量是成本低廉且行之有效的措施。

(2)作物秸稈屬于黃河三角洲地區(qū)的農(nóng)業(yè)廢棄物，易于獲得。采用作物秸稈粉碎后耕翻入田的措施，可以有效改良鹽堿土壤的物理性狀，增加土壤導(dǎo)水性能，快速降低土壤鹽分、提高作物產(chǎn)量。其田間操作技術(shù)簡單，不需要石膏那樣高昂的成本，也避免了利用粉煤灰、各種礦渣對土壤和地下水的污染。又有效避免了大量秸稈覆蓋帶來的作物病蟲害和減產(chǎn)問題。

(3)根據(jù)對土壤理化性狀的綜合分析，結(jié)合試驗小區(qū)收獲的棉花產(chǎn)量，在黃河三角洲地區(qū)，改良粘質(zhì)氯化物鹽化潮土時，作物秸稈粉碎后施入農(nóng)田的適宜用量為200kg/畝。秸稈施入過多會導(dǎo)致土壤大孔隙數(shù)量偏多，雖然土壤鹽分降低，但是土壤的保水、保肥能力下降，不利于作物根系生長，試驗結(jié)果顯示，對棉花出苗影響較大。

(4)有機肥可以提高土壤有機質(zhì)的含量，改善土壤物理性狀，但是由于含有高濃度的氯離子、鈉、鉀和氨等成分，對土壤鹽分含量的降低沒有明顯的效果，另外還對地下水有一定污染。

(5)土壤中大孔隙的數(shù)量與土壤飽和導(dǎo)水率呈現(xiàn)出明顯的正相關(guān)，通過各種農(nóng)藝措施調(diào)控土壤孔隙狀況，可以明顯改善鹽漬土壤的飽和導(dǎo)水率，有助于灌溉后鹽堿土壤中可溶性鹽分的淋洗，快速降低土壤含鹽量，促進作物生長、提高作物產(chǎn)量。

［1］王詮莊，徐樹貞.麥田秸桿覆蓋的作用及其節(jié)水效應(yīng)的初步研究［J］.干旱地區(qū)農(nóng)業(yè)研究，1989，(2):7 －15.

［2］李新舉，張志國.秸桿覆蓋對鹽漬土水分狀況影響的模擬研究［J］.土壤通報，1999，30(4):176 －177.

［3］王德超，姜軍祥.東營市河口區(qū)鹽漬土的改良治理方法及效果分析［J］.山東國土資源，2005，21(5):36－38.

［4］張建鋒，張旭東，周金星，劉國華，李冬雪.世界鹽堿地資源及改良利用的基本措施［J］.水土保持研究，2005，12(6):28－31.

［5］Doran，John C，Turnbull JW.Austranlian Trees and Shrubs:species for land rehabilitation and farm planting in the tropics［M］.Canberra:Austranlian Centre for International Agricultural Research，1997