馬 倩，李海龍

(1.生物地質與環(huán)境地質國家重點實驗室/中國地質大學(北京)水資源與環(huán)境學院，北京 100083;2.地下水循環(huán)與環(huán)境演化教育部重點實驗室/中國地質大學(北京)，北京 100083)

有限元法作為求解地下水流和溶質運移對流-彌散方程的常用數值方法，在水流模擬中有很多優(yōu)點。它有固定的網格，通常滿足質量守恒定律;可以精確、高效地處理以彌散為主的問題;便于編寫程序并得以實現等［1］。但是，對于許多野外常見的以對流為主的問題，求解易引起顯著的數值振蕩。為此，許多學者進行了大量的計算實踐和理論分析指出［2］:對于以對流占優(yōu)的溶質運移問題，有限差分法和有限元法計算時皆不可避免地出現數值振蕩。因此，掌握消除地下水流和溶質運移問題中濃度振蕩的方法和理論依據是很有意義的。

Henry問題研究的是承壓含水層中咸淡水靜力平衡問題，在地下水流和溶質運移領域具有極高的代表性［3］。其常用的概念模型為:承壓含水層均質各向同性，長2m高1m;上下頂板、底板為隔水邊界;內陸邊界為恒定的淡水流量邊界;海岸邊界為定水頭定濃度的海水邊界，海水水位1m［4］。它是一個理想的概念模型，為提高其實際適用性及科研價值(如提高模型參數敏感性)，很多學者在Henry問題的基礎上，研究水文地質條件改變后的海水入侵問題，即變異Henry問題。Frind［5］研究的海水入侵問題，水動力彌散系數與原Henry問題不同;Huyakorn等［6］將右邊界分為定濃度和零濃度通量邊界條件兩部分;Simpon和Clement［3］發(fā)現減小Henry問題中左邊界的內陸淡水補給可以增加密度效應;Abarca等［7］綜合考慮含水層滲透系數的各向異性及右邊界水流方向對濃度的影響，發(fā)現變異后的Henry問題對密度變化更敏感。本研究對內陸和海岸邊界條件改變后的Henry問題進行模擬求解，具體為:左邊界為定水頭邊界、右邊界水位為2m，發(fā)現該變異Henry問題對網格剖分要求更高，對網格Peclet數更為敏感(詳見2.3.3)，為消除濃度振蕩提供了實例，同時，驗證和說明了網格Peclet數在數值計算中的重要性。

MARUN是專門用于地下水流和溶質運移模擬的數值程序，采用Galerkin二維有限元法求解對流-彌散方程，其可靠性已經得到各方面的驗證［8～10］。Li等［11］利用該程序對均質海灘在海潮作用下的地下水鹽動態(tài)進行了有效的定量模擬，Li等［12］、Guo 等［13］、Xia等［14］成功定量擬合了阿拉斯加不同海灣地下水位、鹽分和示蹤劑濃度的野外實測值，同時對海灘地下水動態(tài)水化學特征進行了定量刻畫。

本研究以有限元法數值計算為例，采用MARUN模擬程序對變異Henry問題進行模擬求解，得到了不同網格剖分及水動力彌散系數特選節(jié)點處的濃度穿透曲線，發(fā)現了濃度振蕩現象的存在。同時，分析找到了濃度振蕩的原因及合理的消除方法。模擬結果及分析表明網格Peclet數能夠有效地判定給定的網格剖分是否會引起濃度振蕩，對有限元法數值計算的網格剖分具有指導意義。

1 模型介紹及網格剖分

1.1 模型介紹

變異Henry問題的概念模型具體為:承壓含水層均質各向同性，長2m高1m;上下頂板、底板為隔水邊界;內陸邊界為定水頭邊界，淡水水頭值為2.0632m;海岸邊界為定水頭定濃度的海水邊界，海水水位為2m。原Henry問題的左邊界流量為6.6×10-5m/s，本文提出的變異Henry問題所對應的左邊界流量為1.2×10-4m/s，流速大約是原Henry問題的兩倍。模擬區(qū)域及相關模型參數取值見圖1。

圖1 變異Henry問題的模擬區(qū)域及參數取值Fig．1 Model domain and parameters used in the modified Henry problem

Henry問題使用的水動力彌散系數為D1=1.8857×10-5m2/s，本研究考慮不同水動力彌散系數對模擬結果的影響，故加入水動力彌散系數D2=4.0×10-5m2/s。

1.2 網格剖分

本研究選用MARUN模擬程序采用Galerkin有限單元法計算模擬區(qū)域內各節(jié)點的濃度。剖分區(qū)域為長2m高1m的二維規(guī)則矩形，故采用直角三角形單元剖分，長度方向上的單元數是高度方向上的2倍，保證單元在整個區(qū)域上均勻分布;且所有的三角形單元為“等腰直角三角形”，不存在鈍角三角形單元，是理想的三角形網格。研究共選取兩種網格剖分:［Ⅰ］20行、40列;［Ⅱ］40行、80列。網格剖分Ⅰ共861個節(jié)點，1 600個單元，網格間距0.05m;網格剖分Ⅱ共3 321個節(jié)點，6 400個單元，網格間距0.025m。

2 結果與討論

2.1 模擬求解

使用MARUN程序對以下三種情況的變異Henry問題進行模擬求解，得到各情況模擬區(qū)域濃度的時空變化情況。各節(jié)點的初始濃度為0mg/L，初始水頭由左右邊界水頭線性插值得到。設置初始時間為0s，時間步長以1.1倍的因子增大或縮小，初始時間步長為10s，最大為 1 200s。

Case1:采用網格剖分Ⅰ，水動力彌散系數為D1=1.8857×10-5m2/s;

Case2:采用網格剖分Ⅱ，水動力彌散系數為D1=1.8857×10-5m2/s;

Case3:采用網格剖分Ⅰ，水動力彌散系數為D2=4.0 ×10-5m2/s。

Case1與Case2僅網格剖分不同，模型參數、初始條件及邊界條件均相同。Case1與Case3僅水動力彌散系數不同，其他模型參數、網格剖分、初始條件及邊界條件均相同。

2.2 結果

淡水從內陸邊界流入，海水從海岸邊界驅入含水層，與排泄的淡水混合。內陸和海岸邊界條件始終不變，所以入侵海水與淡水流場最終可以達到動態(tài)平衡狀態(tài)，即得到承壓含水層中鹽水楔形體與淡水流場間的穩(wěn)態(tài)解。本研究模擬區(qū)域內節(jié)點初始濃度均為0mg/L，故各節(jié)點的濃度應先由0mg/L逐漸增加，達到穩(wěn)定狀態(tài)后不再隨時間變化［7］。

有限元法在處理以對流為主的地下水流和溶質運移問題時，網格剖分往往引起數值計算的數值振蕩，當出現很陡的濃度鋒面時，也就是以對流為主時，振蕩尤為明顯［1］。為更清楚地研究和說明濃度振蕩的普遍存在，選取流速較小的節(jié)點A(x=1.5m，z=0.05m)處的濃度穿透曲線進行研究。

2.2.1 不同網格剖分

如圖2所示，Case1、Case2為不同網格剖分節(jié)點A濃度穿透曲線的對比，兩種情況的節(jié)點A濃度由初始時刻的0mg/L近直線上升，至2h左右不再呈現明顯上升趨勢(圖2a);Case1節(jié)點濃度呈現隨時間在某值附近振蕩，并未穩(wěn)定于某一濃度值;Case2節(jié)點濃度最終穩(wěn)定于0.4mg/L(圖2b)。

如前所述，入侵海水與淡水流場達到動態(tài)平衡后，濃度不再隨時間變化。而Case1節(jié)點濃度并未穩(wěn)定，說明使用網格剖分Ⅰ有限元法數值計算出現了數值振蕩，表現為節(jié)點處的濃度振蕩;Case2節(jié)點濃度達到穩(wěn)定，說明使用網格剖分Ⅱ有限元法數值計算未出現數值振蕩，表現為節(jié)點處濃度最終穩(wěn)定于某一濃度值，而不隨時間振蕩。Case1節(jié)點A濃度數值解明顯偏離Case2濃度數值解，表明使用網格剖分Ⅰ有限元數值計算存在數值誤差。即使用網格剖分Ⅰ有限元數值計算時既存在數值誤差，又出現了濃度振蕩。

圖2 三種情況的節(jié)點A處濃度穿透曲線Fig．2 Breakthrough curve of point A at(a)initial time，and(b)whole time of three cases

2.2.2 不同水動力彌散系數

Case1、Case3為不同水動力彌散系數節(jié)點A濃度穿透曲線的對比，兩種情況的節(jié)點A濃度由初始時刻的0mg/L近直線上升，至2h左右不再呈現明顯上升趨勢(圖2a);Case1濃度呈現隨時間在某值附近振蕩，并未穩(wěn)定于某一濃度值;Case3濃度最終穩(wěn)定于0.38mg/L(圖2b)。說明水動力彌散系數為 D1=1.8857×10-5m2/s時，有限元法數值計算出現了數值振蕩，表現為節(jié)點處的濃度振蕩;水動力彌散系數增大為D2=4.0×10-5m2/s時，有限元法數值計算得到了鹽水楔形體與淡水流場的穩(wěn)態(tài)解，數值解未出現數值振蕩，表現為節(jié)點處濃度最終穩(wěn)定于某一濃度值。另外，水動力彌散系數為D2時，若對網格加密后(用網格Ⅱ)進行模擬計算，則粗細網格剖分兩種情況下的節(jié)點處濃度穿透曲線幾乎完全重合，網格加密不再影響計算的精度。說明此時的網格剖分Ⅰ已滿足模擬計算精度要求，得到了精確的穩(wěn)態(tài)數值解。

2.3 討論

2.3.1 濃度振蕩消除方法

Case1、Case2僅網格剖分不同，數值計算得出的特選節(jié)點處濃度穿透曲線特征明顯不同。Case1節(jié)點A濃度振蕩，并未穩(wěn)定于某一濃度值;Case2節(jié)點A濃度最終達到穩(wěn)定。即使用網格剖分Ⅰ時數值計算出現濃度振蕩，網格剖分Ⅱ在網格剖分Ⅰ的基礎上將網格單元整體加密一倍，消除了濃度振蕩。

Case1、Case3僅水動力彌散系數不同，特選節(jié)點處濃度穿透曲線特征明顯不同。Case1節(jié)點濃度振蕩，Case3節(jié)點濃度最終達到穩(wěn)定。即水動力彌散系數為D1時數值計算出現濃度振蕩，而水動力彌散系數D2在D1的基礎上增加約一倍，消除了濃度振蕩。

故有限元法求解地下水流和溶質運移問題時，若出現濃度數值解在某值附近振蕩的現象，可以通過加密網格將其消除，同時須考慮水動力彌散系數設置是否合理，通過增加水動力彌散系數也可消除濃度振蕩。另外，Case1與Case2采用不同網格處理同一問題，說明同一問題對網格精度的要求不同;Case1與Case3均采用網格剖分Ⅰ，僅模型參數不同，說明同一網格對不同模型參數的有效性不同。已有研究結果表明［3～7］，Henry問題在采用網格剖分Ⅰ時能得到精確的數值解，并不出現濃度振蕩，說明即使是研究區(qū)域相同，不同的邊界條件對網格精度的要求不同。網格剖分要具體問題具體分析，不能單純認為邊界條件相同，網格就適用于不同模型參數取值。

2.3.2 理論依據

網格Peclet數(Pe)是衡量對流項在溶質運移問題中主導程度的無量綱參數:

式中:v——達西速度;

D——水動力彌散系數;

Δx——網格間距。

對于純對流問題，D→0，Pe→∞。隨著物理彌散變得明顯，Peclet數則變小。同時，Peclet數還依賴于網格間距Δx;Δx減小時，Peclet數也減小。即網格剖分越細，Peclet數也越小，這意味著采用較小的網格間距可以減小數值振蕩。若Δx對應的Pe＜2，振蕩就消除了［1］。

三種情況下對應的最大網格Peclet數如下:

Case1:采用網格剖分Ⅰ，Δx=0.05m;D=1.8857×10-5m2/s。速度v最大值為0.0015m/s，故模擬區(qū)域內最大網格Peclet數為Pe1=3.97。

Case2:采用網格剖分Ⅱ，Δx=0.025m;D=1.8857×10-5m2/s。速度v最大值為0.0015m/s，故模擬區(qū)域內最大網格Peclet數為Pe2=1.98。

Case3:采用網格剖分Ⅰ，Δx=0.05m;D=4.0×10-5m2/s。速度v最大值為0.0015m/s，故模擬區(qū)域內最大網格Peclet數為Pe3=1.88。

變異Henry問題的內陸和海岸邊界條件不變，模擬區(qū)域內各點的濃度最終應達到動態(tài)平衡狀態(tài)。而Case1濃度數值解出現振蕩現象，與真實的物理意義不符，此時Pe1＞2;Case2濃度數值解最終穩(wěn)定于某一濃度值，達到了動態(tài)平衡狀態(tài)，符合Henry問題真實的物理意義，此時Pe2＜2;Case3濃度數值解最終穩(wěn)定于某一濃度值，符合真實的物理意義，此時Pe3＜2。即不滿足Pe＜2時，使用給定網格剖分進行數值計算時，數值解會出現濃度振蕩(Case1);滿足Pe＜2時，使用給定的網格剖分進行數值計算時，數值解不出現濃度振蕩(Case2、Case3)，說明Peclet數能夠有效地判定給定的網格剖分是否會引起濃度振蕩。使用有限元法進行數值計算時，必須考慮Peclet數，使其盡量小于2。

2.3.3 變異Henry問題與Henry問題的對比

Henry問題與變異Henry問題的區(qū)別是，Henry問題的內陸邊界為定流量邊界，流量為6.6×10-5m/s，變異Henry問題為定水頭邊界，水頭值為2.0632m;Henry問題的海水水位為1m，變異Henry問題的海水水位為2m;變異Henry問題左邊界平均流速為1.2×10-4m/s，大約是原Henry問題左邊界流速的兩倍，故變異Henry問題的對流作用更強，為滿足Peclet數小于2，要求網格剖分更密，即對網格剖分精度要求更高，更適于檢驗網格剖分是否會引起濃度振蕩。海水入侵范圍較Henry問題小，50%等氯線與底板交于1.55m處(Henry問題50%等氯線與底板交于1.40m處［3～7］)。變異Henry問題的求解為消除地下水流和溶質運移問題數值求解中的數值振蕩提供了實例。

地下水數值模擬軟件常用的求解方法是有限元法和有限差分法，二者在數值計算時均可能出現數值振蕩。若出現數值解在某值附近振蕩的現象，可以通過加密網格將其消除，同時需要考慮模型參數設置是否合理，增加水動力彌散系數也可消除。進行網格剖分時，即使研究區(qū)域相同，不同的邊界條件、不同的水動力彌散系數對網格精度的要求不同，同一網格對不同模型參數的有效性也不同，要具體問題具體分析，不可一概而論。

3 結論

(1)濃度振蕩可以通過加密網格、增加水動力彌散系數消除。

(2)即使研究區(qū)域相同，不同的邊界條件、不同的水動力彌散系數對網格精度的要求不同，同一網格對不同模型參數的有效性也不同。

(3)網格Peclet數能夠有效地判定給定的網格剖分是否會引起濃度振蕩，對網格剖分具有指導意義。

(4)在地下水流和溶質運移有限元數值計算中，要綜合考慮模型參數的合理性，分析地下水流場，對流作用較強的區(qū)域要加密剖分，滿足Peclet數小于2，這樣既可以消除濃度振蕩，又保證計算速度。

［1］ 鄭春苗，Gordon D Bennett.地下水污染物遷移模擬［M］.2版.北京:高等教育出版社，2009:109-118.［ZHENG C M，Gordon D Bennett.Applied Contaminant TransportModeling［M］.2nd ed.Beijing:Higher Education Press，2009:109-118.(in Chinese)］

［2］ 馮紹元.解地下水中溶質運移問題的特征有限元法［J］.武漢水利電力學院學報，1991，24(4):410-417.［Feng S Y.Method of Characteristic Finite Element for Solving Solute Transport in an Aquifer［J］.Journal of Wuhan institute of water conservancy and electric power，1991，24(4):410-417.(in Chinese)］

［3］ Simpson M J，T P Clement.Improving the worthiness of the Henry problem as a benchmark for densitydependent groundwaterflow models［J］. Water Resources Research，2004，40:W01504.

［4］ Henry H R.Effects of dispersion on salt encroachment in coastal aquifers［J］.US Geological Survey，1964，Water-Supply Paper 1613-C.

［5］ Frind E O.Simulation of long-term transient densitydependent transport in groundwater［J］.Advances in Water Resources，1982，5:73-88.

［6］ Huyakorn P，Andersen P F，Mercer J，et al.Saltwater intrusion in aquifers:development and testing of a three-dimensional finite element model［J］.Water Resources Research，1987，23(2):293-312.

［7］ Abarca E，Carrera J，Vila X S，et al.Anisotropic dispersvie Henry problem［J］.Advances in Water Resources，2006，30:913-926.

［8］ Boufadel M C，M T Suidan，A D Venosa.A numerical model for density-and-viscosity-dependent flows in two-dimensionalvariably-saturated porous media［J］.Journal of Contaminant Hydrology，1999a，36:1-20.

［9］ Boufadel M C，M T Suidan，A D Venosa.Numerical modeling of water flow below dry salt lakes:effect of capillarity and viscosity［J］.Journal of Hydrology，1999b，221:55-74.

［10］ Boufadel M C.A mechanistic study of nonlinear solute transport in a groundwater-surface water system under steady state and transient hydraulic conditions［J］.Water Resources Research，2000，36:2549-2565.

［11］ Li Hailong，Boufadel M C，Weaver J W.Tideinduced seawater-groundwater circulation in shallow beach aquifers［J］.Journal of Hydrology，2008，352:211-224.

［12］ Li Hailong，Boufadel M C.Long-term persistence of oil from the Exxon Valdez spill in two-layer beaches［J］.Nature Geoscience，2010，3(2):96-99.

［13］ Guo Qiana， Li Hailong， Boufadel M C， et al.Hydrodynamics in a gravel beach and its impact on the Exxon Valdez oil［J］. JournalofGeophysical Research，2010，115:C12077.

［14］ Xia Yiqiang，Li Hailong，Boufadel M C，et al.Hydrodynamic factors affecting the persistence of the Exxon Valdez oil in a shallow bedrock beach［J］.Water Resources Research，2010，46:W10528.