不同網(wǎng)格加密方法在河流回補(bǔ)地下水模擬中的對(duì)比

崔偉哲，郝奇琛，陳 康，陳 飛，唐世南，朱玉晨，曹勝偉

(1. 河北地質(zhì)大學(xué) 水資源與環(huán)境學(xué)院，河北 石家莊 050031； 2. 中國(guó)地質(zhì)科學(xué)院水文地質(zhì)環(huán)境地質(zhì)研究所，河北 石家莊 050061； 3. 水利部水利水電規(guī)劃設(shè)計(jì)總院，北京 100120)

0 引 言

為緩解地下水超采導(dǎo)致的地下水位下降和水資源枯竭等問題，很多國(guó)家實(shí)施了人工回補(bǔ)工程[1-3]。2019年5月，水利部印發(fā)了《2019年度河湖生態(tài)補(bǔ)水方案及試點(diǎn)河段后續(xù)補(bǔ)水計(jì)劃》，明確根據(jù)南水北調(diào)、引黃、引灤等水源條件，向京津冀地區(qū)河湖生態(tài)補(bǔ)水約22.1×108m3，人工回補(bǔ)工程得到高度重視[4]。

地下水模擬是研究地下水運(yùn)動(dòng)的有效方法，可用來進(jìn)行地下水回補(bǔ)模擬研究。國(guó)內(nèi)外許多學(xué)者采用數(shù)值模擬的方法對(duì)人工回補(bǔ)過程、效果及優(yōu)化問題進(jìn)行了研究[5-10]，然而對(duì)于河流這種以線狀回補(bǔ)模擬的精度問題很少提及。與面狀回補(bǔ)模擬不同，河流附近網(wǎng)格剖分方式的不同將會(huì)對(duì)模擬結(jié)果產(chǎn)生影響，進(jìn)而影響回補(bǔ)效果評(píng)估。

GMS(Groundwater Modeling System)軟件及其內(nèi)置的MODFLOW程序是地下水模擬的主流程序[11-13]，從10.0版本開始支持MODFLOW-USG(Unstructured Grid)程序，它是基于體積有限差分法(Control Volume Finite Difference,CVFD)對(duì)網(wǎng)格進(jìn)行剖分[14]。相對(duì)于傳統(tǒng)的MODFLOW程序，MODFLOW-USG程序可以對(duì)局部重點(diǎn)區(qū)域(如河流、開采井)進(jìn)行非結(jié)構(gòu)化網(wǎng)格加密，并且更加靈活。Krcmar等使用 MODFLOW-USG程序?qū)λ孤宸タ薌bely礦區(qū)進(jìn)行模擬，精細(xì)刻畫出與露天礦坑、構(gòu)造斷裂有關(guān)的地質(zhì)非均質(zhì)性[15]。Shokri等通過非結(jié)構(gòu)化網(wǎng)格對(duì)地下水水力梯度較大的非承壓含水層進(jìn)行數(shù)值模擬[16]。苑希民等利用非結(jié)構(gòu)化網(wǎng)格對(duì)黃河南岸灌區(qū)潰口等處進(jìn)行局部加密模擬[17]。周念清等結(jié)合算例模型使用MODFLOW-USG程序?qū)ρ芯繀^(qū)域的源匯項(xiàng)、地下水位和流場(chǎng)分布進(jìn)行數(shù)值模擬[18]。

非結(jié)構(gòu)化網(wǎng)格加密方法包括四叉樹網(wǎng)格(Q-tree Grid)、嵌套網(wǎng)格(Nested Grid)以及泰森多邊形(Voronoi)，目前還沒有專門針對(duì)回補(bǔ)模擬過程中不同網(wǎng)格加密方法精度與效果的對(duì)比研究。其中,泰森多邊形法將網(wǎng)格剖分為不規(guī)則多邊形，與其他兩種規(guī)則矩形的剖分方法很難在網(wǎng)格剖分粗細(xì)上進(jìn)行精確對(duì)比。為此，本文重點(diǎn)對(duì)比四叉樹網(wǎng)格和嵌套網(wǎng)格加密方法，通過與兩種結(jié)構(gòu)化網(wǎng)格(粗、細(xì)結(jié)構(gòu)網(wǎng)格)加密方法的模擬結(jié)果對(duì)比，分析線狀回補(bǔ)模擬中不同網(wǎng)格加密方法的優(yōu)劣，為地下水回補(bǔ)模擬加密方法選擇提供參考。

1 非結(jié)構(gòu)化網(wǎng)格加密方法簡(jiǎn)介

1.1 四叉樹網(wǎng)格加密方法

四叉樹網(wǎng)格加密方法最早應(yīng)用在計(jì)算機(jī)領(lǐng)域，主要用于研究圖像分割[19]。其原理是先將圖像看成是一個(gè)正方形單元，如果該單元內(nèi)有不同性質(zhì)的多邊形，則將單元分成4個(gè)大小相同的二級(jí)單元(圖1)，這種逐級(jí)一分為四的方法一直繼續(xù)到預(yù)定的最高分辨率[20]。Yerry等基于這一原理將四叉樹網(wǎng)格加密方法應(yīng)用于網(wǎng)格生成[21]。丁勝勇等利用四叉樹網(wǎng)格加密方法對(duì)混凝土的界面過渡區(qū)進(jìn)行局部網(wǎng)格加密，給出了更為合理的反映界面過渡區(qū)組分的混凝土細(xì)觀模型[22]。

圖1 四叉樹網(wǎng)格和嵌套網(wǎng)格加密方法對(duì)比Fig.1 Comparison of Q-tree Grid and Nested Grid Refinement Methods

隨著四叉樹網(wǎng)格加密方法的不斷改進(jìn)，其應(yīng)用領(lǐng)域也越來越廣。2001年，英國(guó)牛津大學(xué)Rogers等首次將四叉樹網(wǎng)格應(yīng)用到二維潛水模型[23]。2002年，劉曉東等建立了基于四叉樹網(wǎng)格的Godunov型二維水流數(shù)學(xué)模型，通過實(shí)驗(yàn)證明四叉樹網(wǎng)格比傳統(tǒng)結(jié)構(gòu)網(wǎng)格在復(fù)雜流動(dòng)區(qū)域的模型中具有更高的效率，計(jì)算值較精確[20]。2003年，華祖林等基于四叉樹網(wǎng)格分層對(duì)重點(diǎn)污水排放區(qū)域進(jìn)行模擬[24]。2020年，Bozovic等發(fā)現(xiàn)利用四叉樹網(wǎng)格加密方法的區(qū)域地下水水力梯度最為明顯，能較真實(shí)地刻畫出地下水的三維流動(dòng)[25]。

1.2 嵌套網(wǎng)格加密方法

嵌套網(wǎng)格是Mehl等開發(fā)的一種局部網(wǎng)格加密方法[26]。其原理是將模型分為子模型(加密區(qū))和母模型(非加密區(qū))兩個(gè)獨(dú)立的MODFLOW程序，通過一系列的迭代耦合，并在交界面?zhèn)鬟f水頭和流量信息，實(shí)現(xiàn)模型嵌套[27-28]。

1994年，林振寶等基于FAC(Fast Adaptive Composite)方法，將全局粗結(jié)構(gòu)網(wǎng)格解與局部細(xì)結(jié)構(gòu)網(wǎng)格解之間交替迭代得到的局部加密嵌套網(wǎng)格應(yīng)用到氣藏?cái)?shù)值模擬中[29]。1999年，孫雪梅提出了有限元單元網(wǎng)格加密方法，該方法只需要對(duì)模型所需部位的初始網(wǎng)格進(jìn)行局部加密[30]。2013年，周巍等利用三層嵌套網(wǎng)格較精確地模擬了珠江口的主要水動(dòng)力過程和鹽度分布[31]。2015年，楊楊等對(duì)模型嵌套源代碼進(jìn)行二次開發(fā)，使得該程序可以對(duì)邊界不規(guī)則的區(qū)域進(jìn)行局部嵌套網(wǎng)格加密，在保證計(jì)算精度的情況下縮短了運(yùn)行時(shí)間[28]。2017年，齊歡等通過GMS軟件建立了趵突泉泉域的地下水流數(shù)值模型，利用嵌套網(wǎng)格加密方法對(duì)重點(diǎn)區(qū)域進(jìn)行局部加密[27]。

1.3 非結(jié)構(gòu)化網(wǎng)格加密的優(yōu)勢(shì)

在區(qū)域性地下水模型中，網(wǎng)格通常會(huì)剖分比較大[32]，在有開采井或河流這些水位變化明顯的地方，模型不能很好地模擬出局部水位變化，這就需要對(duì)重點(diǎn)區(qū)域進(jìn)行加密剖分。非結(jié)構(gòu)化網(wǎng)格加密要比結(jié)構(gòu)化網(wǎng)格加密更加靈活，在確保計(jì)算精度的情況下能有效減小計(jì)算量。韋未等通過兩個(gè)算例證明局部網(wǎng)格加密可在確保計(jì)算精度的情況下,很大程度地減少計(jì)算量[33]。童少青等分別用結(jié)構(gòu)化網(wǎng)格全區(qū)域加密和非結(jié)構(gòu)化網(wǎng)格局部加密建立黑河中游盆地地下水流動(dòng)模型，發(fā)現(xiàn)非結(jié)構(gòu)化網(wǎng)格加密方法能夠提高模型仿真度與計(jì)算結(jié)果精度[34]。

2 地下水回補(bǔ)算例模型

算例模型是根據(jù)水利部地下水管理與保護(hù)項(xiàng)目《河湖地下水回補(bǔ)試點(diǎn)補(bǔ)水效果》數(shù)值模擬中的拒馬河模型概化而來。除模型范圍外，算例模型中的水文地質(zhì)結(jié)構(gòu)、參數(shù)、邊界條件等在概化過程中盡量保持與原始模型具有一定相似性，以求較為客觀地反映實(shí)例情況。

2.1 算例模型建立

假設(shè)河流所在的區(qū)域面積為10 km×10 km的含水層，地下水總體流向?yàn)樽晕飨驏|。根據(jù)研究區(qū)的水文地質(zhì)條件，將含水層概化為非均質(zhì)各向異性潛水、中層承壓水、深層承壓水等3層含水系統(tǒng)。其厚度分別為50、100、100 m，水平滲透系數(shù)分別為50、20、7 m·d-1，垂向滲透系數(shù)分別為0.08、0.04、0.02 m·d-1。含水層組間的水力聯(lián)系強(qiáng)弱由模擬層間的垂向滲透系數(shù)控制。第1層給水度為0.08，第2層、第3層儲(chǔ)水率均為0.000 1。東側(cè)和西側(cè)與外界存在明顯的水量交換，將西側(cè)設(shè)定為流量邊界，東側(cè)設(shè)定為給定水頭邊界，其他設(shè)定為通用水頭邊界。主要補(bǔ)給來源為大氣降水，降雨強(qiáng)度為0.003 m·d-1。主要排泄項(xiàng)為人工開采，共布有9口開采井(圖2)，其中6口井位于第1層，2口井位于第2層，1口井位于第3層，每口井的開采量設(shè)定為4 000 m3·d-1。模擬期為1年，半個(gè)月設(shè)置為一個(gè)應(yīng)力期，1 d為一個(gè)步長(zhǎng)。河道位于模型的第1層，河道的平均回補(bǔ)量設(shè)定為1.5×105m3·d-1，回補(bǔ)時(shí)間為1年。算例模型的水文地質(zhì)結(jié)構(gòu)、初始條件、邊界條件等如圖2所示。

圖2 地下水回補(bǔ)算例模型結(jié)構(gòu)及設(shè)置Fig.2 Structure and Setting for the Groundwater Recharge Model

2.2 4種網(wǎng)格剖分方案

圖3 4種網(wǎng)格剖分方案Fig.3 Four Grid Generation Methods

為對(duì)比不同加密方法的誤差，共設(shè)置了4種網(wǎng)格剖分方案。方案一(粗結(jié)構(gòu)網(wǎng)格)：網(wǎng)格剖分為500 m×500 m，網(wǎng)格數(shù)量為1 200[圖3(a)]。方案二(細(xì)結(jié)構(gòu)網(wǎng)格)：為盡量精細(xì)刻畫地下水運(yùn)動(dòng)過程，網(wǎng)格剖分為15.625 m×15.625 m，網(wǎng)格數(shù)量為1 228 800[圖3(b)]。方案三(四叉樹網(wǎng)格)：局部河道處最小網(wǎng)格為125 m×125 m，由河道向兩側(cè)漸變?yōu)?00 m×500 m，網(wǎng)格數(shù)量為2 586[圖3(c)]。方案四(嵌套網(wǎng)格)：局部河道處網(wǎng)格大小為125 m×125 m，其他區(qū)域網(wǎng)格大小為500 m×500 m，網(wǎng)格數(shù)量為3 135[圖3(d)]。

2.3 模擬情景

為了便于對(duì)比河流回補(bǔ)條件下的差別，設(shè)置了兩種模擬情景，即無河道回補(bǔ)與有河道回補(bǔ)。兩種情景中除有無河流回補(bǔ)的差別之外，其他條件均保持一致。以細(xì)結(jié)構(gòu)網(wǎng)格模型作為標(biāo)準(zhǔn)模型，通過兩種情景的模擬，對(duì)比研究不同網(wǎng)格加密方法對(duì)模擬結(jié)果的影響，包括末流場(chǎng)、水位及水量等方面，以此表征不同模型計(jì)算結(jié)果的精確程度[35]。

3 模擬結(jié)果分析及討論

3.1 末流場(chǎng)對(duì)比

圖(a)中數(shù)據(jù)為細(xì)結(jié)構(gòu)網(wǎng)格模型的水位；圖(b)、(c)、(d)中數(shù)據(jù)為4種網(wǎng)格模型的水位圖4 模擬末流場(chǎng)對(duì)比Fig.4 Comparisons of Simulated Final Flow Fields

對(duì)比無河流回補(bǔ)與有河流回補(bǔ)條件下1年后細(xì)結(jié)構(gòu)網(wǎng)格模型的淺層末流場(chǎng)[圖4(a)]，可直觀地看出在河流回補(bǔ)條件下，河道附近的水頭值上升明顯，回補(bǔ)效果較為顯著。通過對(duì)比4種網(wǎng)格模型回補(bǔ)1個(gè)月、6個(gè)月、12個(gè)月后的末流場(chǎng)[圖4(b)～(d)]可以看出，整體上4種網(wǎng)格剖分方案模擬的末流場(chǎng)大體一致，大部分地區(qū)等水位線基本重合，但在河道附近，模擬的末流場(chǎng)形態(tài)有一定差異，特別是回補(bǔ)初期，粗結(jié)構(gòu)網(wǎng)格模型的誤差較大。隨著回補(bǔ)時(shí)間的延長(zhǎng)，河道附近的誤差變小，但遠(yuǎn)離河道的區(qū)域開始出現(xiàn)一定的誤差，特別是在邊界處。

為能更精確地看出網(wǎng)格模型之間的末流場(chǎng)差異，采用定量化評(píng)價(jià)指標(biāo)，以細(xì)結(jié)構(gòu)網(wǎng)格模型的末流場(chǎng)為標(biāo)準(zhǔn)，計(jì)算出位于河道500 m范圍內(nèi)粗結(jié)構(gòu)網(wǎng)格、嵌套網(wǎng)格、四叉樹網(wǎng)格模型與標(biāo)準(zhǔn)網(wǎng)格模型末流場(chǎng)的標(biāo)準(zhǔn)差(表1)。河流回補(bǔ)1個(gè)月后，粗結(jié)構(gòu)網(wǎng)格模型的標(biāo)準(zhǔn)差要比其他兩種非結(jié)構(gòu)化網(wǎng)格模型大0.2左右。回補(bǔ)12個(gè)月后，粗結(jié)構(gòu)網(wǎng)格模型的標(biāo)準(zhǔn)差與四叉樹網(wǎng)格模型的差距依然保持在0.2以上，而與嵌套網(wǎng)格模型的標(biāo)準(zhǔn)差差距逐漸降為0.14。在整個(gè)回補(bǔ)過程中，與粗結(jié)構(gòu)網(wǎng)格模型相比，兩種非結(jié)構(gòu)化網(wǎng)格模型的末流場(chǎng)都要更接近標(biāo)準(zhǔn)末流場(chǎng)。而隨回補(bǔ)時(shí)間的增加，四叉樹網(wǎng)格模型與嵌套網(wǎng)格模型的標(biāo)準(zhǔn)差也在逐漸拉大，回補(bǔ)12個(gè)月后四叉樹網(wǎng)格模型的標(biāo)準(zhǔn)差要比嵌套網(wǎng)格模型小0.07。從兩種非結(jié)構(gòu)化網(wǎng)格模型的標(biāo)準(zhǔn)差對(duì)比來看，在回補(bǔ)期末，四叉樹網(wǎng)格模型的末流場(chǎng)要比嵌套網(wǎng)格模型更接近標(biāo)準(zhǔn)末流場(chǎng)。

表1 不同網(wǎng)格模型標(biāo)準(zhǔn)差對(duì)比

根據(jù)上述結(jié)論，統(tǒng)計(jì)出位于河道500 m范圍內(nèi)兩種非結(jié)構(gòu)化網(wǎng)格模型與標(biāo)準(zhǔn)網(wǎng)格模型的3期(回補(bǔ)1個(gè)月、6個(gè)月、12個(gè)月)誤差，并統(tǒng)計(jì)出誤差所占的比例分布(圖5)。河流回補(bǔ)1個(gè)月后，兩種非結(jié)構(gòu)化網(wǎng)格模型的誤差基本都保持在0.1 m內(nèi)(分別為98.83%、99.77%)。隨著回補(bǔ)時(shí)間的增加，兩種非結(jié)構(gòu)化網(wǎng)格模型的誤差在0.1 m以上的比例也都在逐漸增加。回補(bǔ)6個(gè)月后，嵌套網(wǎng)格模型在誤差為0.5 m以上的比例是四叉樹網(wǎng)格模型的2.7倍，并且此時(shí)嵌套網(wǎng)格模型已經(jīng)開始出現(xiàn)誤差大于1 m的區(qū)域，四叉樹網(wǎng)格模型從回補(bǔ)初到回補(bǔ)末，誤差一直保持在1 m內(nèi)。這也能夠解釋隨著回補(bǔ)時(shí)間的增加，四叉樹網(wǎng)格模型與嵌套網(wǎng)格模型的誤差差距逐漸變大的原因，同時(shí)也進(jìn)一步說明四叉樹網(wǎng)格模型在長(zhǎng)期模擬過程中的優(yōu)勢(shì)。

3.2 水位相對(duì)上升差異

3.2.1 4種網(wǎng)格剖分方案對(duì)比

圖6為4種網(wǎng)格模型回補(bǔ)條件下的水位相對(duì)上升情況(上升值為同一種剖分方式下，有回補(bǔ)條件下模擬的末流場(chǎng)水位減去無回補(bǔ)條件下模擬的末流場(chǎng)水位)。整體上看，4種網(wǎng)格剖分方案的模擬結(jié)果大體一致，水位上升范圍沿河道呈條帶狀分布，最遠(yuǎn)影響距離約為5.9 km(以水位變差1 m計(jì)算)，回補(bǔ)河段上游影響距離較大。

從形態(tài)上看，粗結(jié)構(gòu)網(wǎng)格模型的模擬結(jié)果與其他3種網(wǎng)格模型的模擬結(jié)果有一定差異，主要表現(xiàn)為等水位線不夠圓滑，下游一端的梭形形態(tài)未體現(xiàn)出來。細(xì)結(jié)構(gòu)網(wǎng)格模型與其他兩種非結(jié)構(gòu)化網(wǎng)格模型的模擬結(jié)果較相似，形態(tài)較圓滑。

從最高變幅看，粗結(jié)構(gòu)網(wǎng)格模型的模擬結(jié)果與其他3種網(wǎng)格模型也有較大差別。細(xì)結(jié)構(gòu)網(wǎng)格模型的最大水位變幅約為9 m，且呈條帶狀分布在河道兩側(cè)，9 m等水位線圈閉的面積為0.74 km2，其他兩種局部加密的非結(jié)構(gòu)化網(wǎng)格模型模擬結(jié)果也類似，9 m等水位線圈閉的面積分別為0.68、0.66 km2；但粗結(jié)構(gòu)網(wǎng)格模型的9 m等水位線圈閉的面積只有0.33 km2，相比細(xì)結(jié)構(gòu)網(wǎng)格模型小了55%，且圈閉范圍不連續(xù)，甚至8 m等水位線也不圈閉。

綜上所述，相比粗結(jié)構(gòu)網(wǎng)格模型，其他3種網(wǎng)格模型均能有效提高模擬精度，但三者也略有差別，以下將從不同距離的水位變差、單個(gè)觀測(cè)孔的水位過程線以及均衡情況進(jìn)一步論述。

為了更直觀體現(xiàn)不同方案的差別，統(tǒng)計(jì)出4種網(wǎng)格剖分方案中距離河道30～800 m的水位變差(圖7)。從整體變化趨勢(shì)上看，4種網(wǎng)格模型的水位變差變化趨勢(shì)基本相同。隨著與河道距離的增加，水位變差不斷減小，受河流回補(bǔ)影響越來越小，同時(shí)水位變差速率也逐漸減小。但在距河道300 m范圍內(nèi)，粗結(jié)構(gòu)網(wǎng)格模型的模擬結(jié)果與其他3種網(wǎng)格模型有較大差別。在距離河道30 m范圍內(nèi)，細(xì)結(jié)構(gòu)網(wǎng)格模型的水位變差約為8.41 m，兩種局部加密的非結(jié)構(gòu)化網(wǎng)格模型模擬結(jié)果與其較接近，分別為8.21、8.17 m，而粗結(jié)構(gòu)網(wǎng)格模型的水位變差只有7.89 m，比細(xì)結(jié)構(gòu)網(wǎng)格模型小6%。在回補(bǔ)效果評(píng)估過程中，這種差異一定程度上會(huì)導(dǎo)致評(píng)估效果的不準(zhǔn)確。

圖6 回補(bǔ)條件下的水位相對(duì)上升情況Fig.6 Relative Increases of Water Level Under Recharge Conditions

圖7 距離河道不同范圍的水位變差對(duì)比Fig.7 Comparison of Water Level Variations Under Different Distances with River

在距離河道上游100 m、中游200 m處設(shè)置兩個(gè)觀測(cè)孔(圖2)，對(duì)比4種網(wǎng)格模型的地下水位動(dòng)態(tài)(圖8)。兩種局部加密的非結(jié)構(gòu)化網(wǎng)格模型與細(xì)結(jié)構(gòu)網(wǎng)格模型的水位動(dòng)態(tài)更接近。在回補(bǔ)期末，距離河道100 m觀測(cè)孔的細(xì)結(jié)構(gòu)網(wǎng)格模型的水位變差要比粗結(jié)構(gòu)網(wǎng)格模型的水位變差高1.21 m。從時(shí)間變化來看，隨著時(shí)間的推移，受河道回補(bǔ)的影響，水位逐漸升高，升高速率逐漸變緩，但4種網(wǎng)格模型的水位變差的差距卻在擴(kuò)大。從距離河道不同范圍的觀測(cè)孔來看，距離河道200 m觀測(cè)孔的4種網(wǎng)格模型的水位動(dòng)態(tài)差距要小于距離河道100 m處的觀測(cè)孔。

距離河道不同范圍的水位變差對(duì)比(圖7)綜合分析表明：在有河道回補(bǔ)的地下水模型中，河道所在網(wǎng)格剖分較粗時(shí)，由回補(bǔ)引起的附近水位變化會(huì)被平緩化，在模型中表現(xiàn)不突出；當(dāng)對(duì)河道所在網(wǎng)格進(jìn)行加密后，河道所在網(wǎng)格體積減小，其附近網(wǎng)格數(shù)量相對(duì)增加，地下水位的變化受河道回補(bǔ)影響也就更加明顯；局部加密的非結(jié)構(gòu)化網(wǎng)格模型能夠更加精細(xì)刻畫出局部的水位變化。

圖8 距離河道100、200 m觀測(cè)孔水位動(dòng)態(tài)對(duì)比Fig.8 Comparisons of Water Level Change for Observation Wells Under the Distances of 100 and 200 m with River

3.2.2 不同邊界條件對(duì)誤差的影響

邊界條件是數(shù)值模擬中一項(xiàng)重要設(shè)置，以上研究結(jié)果都是基于西側(cè)補(bǔ)給區(qū)邊界為二類邊界(流量邊界)得出的，因此，需要通過改變不同的邊界條件來分析其對(duì)誤差的影響。在補(bǔ)給區(qū)改用一類邊界(給定水頭邊界)、三類邊界(通用水頭邊界)，模型的其他條件或參數(shù)均保持不變。以細(xì)結(jié)構(gòu)模型在距離河道不同范圍的水位變差為標(biāo)準(zhǔn)，分別計(jì)算出不同邊界條件下，其他3種網(wǎng)格模型模擬的水位變差與細(xì)結(jié)構(gòu)網(wǎng)格模型之間的誤差(圖9)。

圖9 不同邊界條件下水位變差的誤差對(duì)比Fig.9 Comparisons of Water Level Variation Errors Under Different Boundary Conditions

3類邊界條件下，3種網(wǎng)格模型之間的誤差變化趨勢(shì)基本相同。距離河道30 m處，粗結(jié)構(gòu)網(wǎng)格模型的誤差均在0.5 m以上，兩種非結(jié)構(gòu)化網(wǎng)格模型的誤差相對(duì)較小，保持在0.2～0.3 m；隨著與河道距離的增加，粗結(jié)構(gòu)網(wǎng)格模型與非結(jié)構(gòu)化網(wǎng)格模型的誤差逐漸減小，距離河道400 m后，3種網(wǎng)格模型的誤差基本趨于一致。這說明不管在什么邊界條件下，非結(jié)構(gòu)化網(wǎng)格模型的模擬精度都要高于粗結(jié)構(gòu)網(wǎng)格模型。同時(shí)，從整體上看，3類邊界條件下，四叉樹網(wǎng)格模型的誤差要比嵌套網(wǎng)格模型偏小0.01～0.04 m，模擬精度要稍高于嵌套網(wǎng)格模型。

3.2.3 不同加密級(jí)次對(duì)誤差的影響

四叉樹網(wǎng)格模型河道附近最小網(wǎng)格為125 m×125 m，最大網(wǎng)格為500 m×500 m，屬于三級(jí)加密。將其改為六級(jí)加密，局部最小網(wǎng)格剖分為15.625 m×15.625 m(網(wǎng)格數(shù)量為15 771)，同時(shí)將嵌套網(wǎng)格模型的最小網(wǎng)格也剖分為15.625 m×15.625 m(網(wǎng)格數(shù)量為133 167)。進(jìn)一步加密后的模型參數(shù)與之前的參數(shù)均保持一致，對(duì)比非結(jié)構(gòu)化網(wǎng)格模型與標(biāo)準(zhǔn)網(wǎng)格模型在距離河道不同范圍的水位變差(圖10)。

從模擬精度上看，進(jìn)一步加密后的四叉樹網(wǎng)格模型和嵌套網(wǎng)格模型在距離河道30 m處的水位變差相對(duì)于之前的模型分別僅增加了0.04 m和0.03 m，模擬精度提升不大。在距離河道300 m后，非結(jié)構(gòu)化網(wǎng)格模型的水位變差基本保持一致。

考慮到模型的網(wǎng)格數(shù)量，進(jìn)一步加密后的嵌套網(wǎng)格模型網(wǎng)格數(shù)量要比之前高42倍，四叉樹網(wǎng)格模型網(wǎng)格數(shù)量也要比之前高6倍。顯然，為提高模擬精度再次加密局部網(wǎng)格的非結(jié)構(gòu)化網(wǎng)格模型意義不大。另外，對(duì)局部網(wǎng)格再次加密后，嵌套網(wǎng)格模型要比四叉樹網(wǎng)格模型的網(wǎng)格數(shù)量多7倍，再次說明四叉樹網(wǎng)格模型具有多級(jí)加密且級(jí)次越多，優(yōu)勢(shì)越明顯的特點(diǎn)。

圖10 網(wǎng)格加密前后在距離河道不同范圍的水位變差對(duì)比Fig.10 Comparison of Water Level Variations Before and After Grid Refinement Under Different Distances with River

3.3 地下水側(cè)向徑流量變化

在河道兩側(cè)200 m范圍建立均衡區(qū)，對(duì)比第2.2節(jié)中4種網(wǎng)格剖分方案回補(bǔ)期末均衡區(qū)內(nèi)的側(cè)向徑流量(表2)。從表2可以看出，細(xì)結(jié)構(gòu)網(wǎng)格模型與兩種非結(jié)構(gòu)化網(wǎng)格模型的模擬結(jié)果較為接近，但這3種網(wǎng)格模型與粗結(jié)構(gòu)網(wǎng)格模型的模擬結(jié)果差異較大，粗結(jié)構(gòu)網(wǎng)格模型比細(xì)結(jié)構(gòu)網(wǎng)格模型的側(cè)向徑流量少約15%。上述結(jié)果進(jìn)一步驗(yàn)證了前述論證，即回補(bǔ)模擬過程中，粗結(jié)構(gòu)網(wǎng)格模型的水位升高值要比其他3種網(wǎng)格模型低，故水力梯度也小，最終導(dǎo)致側(cè)向徑流量明顯偏小。

表2 4種網(wǎng)格剖分方案?jìng)?cè)向徑流量對(duì)比

3.4 模型運(yùn)行效率

在實(shí)際應(yīng)用中，地下水模型運(yùn)行效率尤為重要。由于軟件本身和計(jì)算機(jī)能力的限制，在建模過程中，當(dāng)網(wǎng)格數(shù)量過多時(shí)，很容易出現(xiàn)模型運(yùn)行效率低下甚至無法運(yùn)行的情況，給局部網(wǎng)格加密帶來很大困難。本次研究模型為前述建立的河流回補(bǔ)模型，對(duì)4種網(wǎng)格剖分方案的運(yùn)行效率進(jìn)行了對(duì)比分析，采用的計(jì)算機(jī)CPU為4核1.8 GHz，RAM為8 GB。計(jì)算結(jié)果如表3所示。

細(xì)結(jié)構(gòu)網(wǎng)格模型相對(duì)于粗結(jié)構(gòu)網(wǎng)格模型進(jìn)行了全區(qū)域整體加密，在模擬結(jié)果上更精確，但由于剖分的網(wǎng)格數(shù)量過多，運(yùn)行時(shí)間高達(dá)5 h，并且運(yùn)行后會(huì)產(chǎn)生40 GB的結(jié)果文件，占據(jù)了大量計(jì)算機(jī)存儲(chǔ)空間。兩種非結(jié)構(gòu)化網(wǎng)格模型只是在河道附近區(qū)域進(jìn)行局部加密，相對(duì)粗結(jié)構(gòu)網(wǎng)格模型，網(wǎng)格數(shù)量增加不大，由1 200分別增加到2 586和3 135，運(yùn)行時(shí)間有一定程度的增加。另外，兩種非結(jié)構(gòu)化網(wǎng)格模型之間相比較，四叉樹網(wǎng)格模型(三級(jí)加密)要比嵌套網(wǎng)格模型(局部為125 m×125 m)節(jié)省18%的運(yùn)行時(shí)間。這是由于嵌套網(wǎng)格模型是對(duì)河道所在區(qū)域均勻加密，而四叉樹網(wǎng)格模型采用的是漸變網(wǎng)格剖分方法進(jìn)行局部加密，網(wǎng)格數(shù)量相對(duì)較少。四叉樹網(wǎng)格模型(六級(jí)加密)的運(yùn)行時(shí)間僅是嵌套網(wǎng)格模型(局部為15.625 m×15.625 m)的15%，節(jié)省了85%的運(yùn)行時(shí)間，具有多級(jí)加密的優(yōu)勢(shì)。

表3 不同網(wǎng)格剖分方案運(yùn)行時(shí)間對(duì)比

4 結(jié) 語(yǔ)

(1)通過對(duì)比分析不同網(wǎng)格剖分方案下模擬的末流場(chǎng)、距離河道不同范圍回補(bǔ)前后水位變差、單個(gè)觀測(cè)孔的水位過程線、均衡情況、模型運(yùn)行效率等，表明兩種局部加密的非結(jié)構(gòu)化網(wǎng)格模型以及細(xì)結(jié)構(gòu)網(wǎng)格模型比粗結(jié)構(gòu)網(wǎng)格模型的模擬精度高，更能準(zhǔn)確刻畫因河道回補(bǔ)導(dǎo)致的地下水位變化。對(duì)于長(zhǎng)期模擬的模型而言，四叉樹網(wǎng)格模型相對(duì)于嵌套網(wǎng)格模型更具優(yōu)勢(shì)，且適用于不同類型的邊界條件。

(2)兩種局部加密的非結(jié)構(gòu)化網(wǎng)格模型比細(xì)結(jié)構(gòu)網(wǎng)格模型網(wǎng)格數(shù)量少，運(yùn)行時(shí)間短，在沒有大幅降低模擬精度的前提下，運(yùn)行效率更高，且四叉樹網(wǎng)格模型具有多級(jí)加密的優(yōu)勢(shì)，加密級(jí)次越多，優(yōu)勢(shì)越明顯，運(yùn)行時(shí)間相對(duì)更短。

(3)在模擬河流生態(tài)補(bǔ)水對(duì)地下水影響時(shí)，四叉樹網(wǎng)格加密方法的模擬精度和運(yùn)行效率相對(duì)較高，具有一定優(yōu)勢(shì)，是值得推薦的一種局部網(wǎng)格加密方法。