2.1 概述

MODFLOW 代碼利用以塊為中心的有限差分方法來求解地下水流方程。流域被離散成行、列和層，使得每個(gè)節(jié)點(diǎn)代表多孔材料的矩形塊，稱為單元。生成的有限差分網(wǎng)格中的節(jié)點(diǎn)表示無流量、可變水頭或恒定水頭單元，并且與單元相關(guān)的任何水力屬性都是相對(duì)于相應(yīng)的節(jié)點(diǎn)指定的。

對(duì)于無側(cè)限流模擬，原始 MODFLOW 代碼（McDonald 和 Harbaugh，1988）的 BCF1 包將去飽和可變頭單元轉(zhuǎn)換為無流單元。這會(huì)導(dǎo)致模型模擬中永久排除部分流域。如果電位水平恢復(fù)，MODFLOW 無法將這些單元恢復(fù)為可變頭單元，這可能會(huì)導(dǎo)致代碼產(chǎn)生誤導(dǎo)或錯(cuò)誤的模擬結(jié)果。

McDonald 等人 (1991) 嘗試將無流動(dòng)池重新轉(zhuǎn)換為可變頭池，以便在需要時(shí)允許池重新飽和。他們的再潤(rùn)濕方案在 BCF2 包（塊中心流包，版本 2）中實(shí)施。該方案利用相鄰單元處的頭來確定是否將無流單元重新轉(zhuǎn)換為可變頭單元。不幸的是，BCF2 再潤(rùn)濕選項(xiàng)在再潤(rùn)濕期間或當(dāng)撤回使各個(gè)單元變干時(shí)容易出現(xiàn)收斂和穩(wěn)定性問題（McDonald 等人，1991）。這是因?yàn)橛糜谥匦录せ罡杉?xì)胞的臨時(shí)程序可能嚴(yán)重違反流量和質(zhì)量守恒定律。此外，當(dāng)相鄰單元處的頭部相關(guān)透射率減小到零時(shí)，代碼中使用的諧波透射率平均方案在物理上變得不一致。 Goode 和 Appel (1992) 記錄了與 BCF2 再潤(rùn)濕包相關(guān)的一系列問題以及一系列緩解方案，他們將 BCF3 包添加到 MODFLOW 中，以提供平均透射率的替代程序。鑒于上述再潤(rùn)濕程序遇到的嚴(yán)重困難，需要一種嚴(yán)格的方法來滿足流動(dòng)連續(xù)性要求并允許地下水位在無約束層中自由運(yùn)動(dòng)，而不需要任何形式的強(qiáng)制轉(zhuǎn)換。 HydroGeoLogic, Inc. 開發(fā)了這樣一種創(chuàng)新方法，通過使用 3-D 可變飽和流公式和自動(dòng)生成的偽土壤保水函數(shù)，將非飽和流問題減少到尋求地下水位（即，電池中的壓頭為零或大氣壓）。該公式旨在準(zhǔn)確描繪地下水位，并捕獲無約束系統(tǒng)對(duì)抽水和補(bǔ)給的延遲產(chǎn)量響應(yīng)。

模擬的數(shù)據(jù)要求沒有改變，因?yàn)槭褂昧藗瓮寥狸P(guān)系而不是真實(shí)土壤的本構(gòu)關(guān)系。我們方法的主要優(yōu)點(diǎn)包括：

? 對(duì) 3-D 流進(jìn)行嚴(yán)格處理（即，沒有 Dupuit 假設(shè)），不需要打開和關(guān)閉單元。

? 穩(wěn)健的數(shù)值方案，具有良好的收斂性和穩(wěn)定性特征以及處理細(xì)胞完全去飽和和再飽和的所需能力。

使用有限差分法，將帶有偽土函數(shù)的可變飽和公式編程到 MODFLOW 現(xiàn)有的 BCF3 包中，新版本稱為 BCF4 包。 BCF4 包的新選項(xiàng)將塊間電導(dǎo)計(jì)算為塊水力傳導(dǎo)率、相對(duì)滲透率和平均流通面積的加權(quán)調(diào)和平均值的乘積。以下各節(jié)討論了 BCF4 包中具有擬土壤函數(shù)的可變飽和公式的使用及其應(yīng)用。

BCF4 軟件包還提供使用軸對(duì)稱圓柱坐標(biāo)系執(zhí)行軸對(duì)稱分析的選項(xiàng)。此選項(xiàng)避免使用通過笛卡爾 (x,y,z) 坐標(biāo)的完全 3-D 公式，從而大大節(jié)省計(jì)算量。軸對(duì)稱分析選項(xiàng)適用于可接受軸對(duì)稱幾何形狀的某些情況（即單井流動(dòng)）。

MODFLOW 的模塊化特性和輸入數(shù)據(jù)準(zhǔn)備的結(jié)構(gòu)在 BCF4 包的開發(fā)中得以保留。新功能的數(shù)據(jù)要求與原始 MODFLOW 相同，不需要額外的參數(shù)值。為了保持向后兼容性并允許用戶比較各種計(jì)算方案，BCF4 軟件包還合并了以前的 BCF3 和 BCF2 軟件包的所有選項(xiàng)。在本章中，通過對(duì)公式化、驗(yàn)證和演示示例問題的簡(jiǎn)要討論，向用戶系統(tǒng)地介紹了 BCF4 包的新功能。

2.2A 使用變飽和流動(dòng)公式建立無側(cè)限流動(dòng)方程

水在可變飽和系統(tǒng)中的 3-D 運(yùn)動(dòng)可以表示為（Huyakorn 等人，1986）：

其中：

x、y 和 z 是笛卡爾坐標(biāo) (L)；?

Kxx、Kyy 和 Kzz 分別是沿 x、y 和 z 軸的水力傳導(dǎo)率；

krw 是相對(duì)滲透率，是水飽和度的函數(shù)；

h為水頭（L）；?

W 是每單位體積的體積通量，代表水源和/或匯；?

N 是可排水孔隙度，等于比產(chǎn)量；

Sw是水的飽和度，是壓頭的函數(shù)；?

Ss 為多孔材料的比儲(chǔ)存量；?

t 為時(shí)間（T）。

對(duì)于完全飽和介質(zhì)（即 Sw = 1.0），相對(duì)滲透率為 1，方程 (1) 簡(jiǎn)化為：

方程（2）是 MODFLOW 開發(fā)中使用的基本地下水流方程（McDonald 和 Harbaugh，1988）。因此，可變飽和流量方程簡(jiǎn)化為地下水位以下和受限系統(tǒng)中的傳統(tǒng)地下水流量方程。擬土關(guān)系用于定義函數(shù)關(guān)系：Sw = Sw (R) 和 krw = krw (Sw)，其中 R 是壓頭，定義為 R = h+z，其中 z 在垂直向下方向?yàn)檎怠Ｕ?qǐng)注意，方程 (2) 的以塊為中心的有限差分近似包含塊間水力電導(dǎo)率，這些電導(dǎo)率是作為參與單元的節(jié)點(diǎn)值的加權(quán)調(diào)和平均值獲得的。

因此，對(duì)于 BCF4 包建模選項(xiàng)，用戶指定水平飽和水力傳導(dǎo)率（Kxx 和 Kyy）以及每個(gè)單元的頂部和底部高程。塊間飽和水力傳導(dǎo)率乘以相對(duì)滲透率即可得到塊間有效水力傳導(dǎo)率。垂直塊間泄漏（垂直水力傳導(dǎo)率 Kzz 除以節(jié)點(diǎn)之間的垂直距離 )z）由用戶直接輸入，就像在所有以前的 BCF 包中所做的那樣。這些泄漏可以使用原始 MODFLOW 文檔（McDonald 和 Harbaugh，1988）中討論的各種方法來計(jì)算。使用改進(jìn)的 Picard 方法處理變飽和流量方程中的非線性（Celia 等，1990）。與之前的 BCF 包類似，用戶輸入無量綱主存儲(chǔ)系數(shù) SF1。等式(1)使用的特定存儲(chǔ)(Ss)是通過將存儲(chǔ)系數(shù)SF1除以塊厚度來內(nèi)部計(jì)算的。

總之，BCF4 包中嵌入的具有擬土函數(shù)的可變飽和公式包括比以前的包有兩個(gè)改進(jìn)的方案： (1) 一種解決無側(cè)限流動(dòng)問題的新方法，它將去飽和和再飽和視為控制方程的一個(gè)組成部分但不需要土壤保水?dāng)?shù)據(jù)； (2) 附錄 A 中所示的飽和水力電導(dǎo)率的調(diào)和平均比有效（水頭相關(guān)）透射率的調(diào)和平均更合適。新選項(xiàng)不會(huì)使干電池失活，而是計(jì)算通過非飽和區(qū)域傳輸充電水所需的水頭。請(qǐng)注意，新的 BCF4 軟件包可用于全 3D、準(zhǔn) 3D 和軸對(duì)稱模擬。

對(duì)于軸對(duì)稱流模擬，BCF4 包利用以下形式的飽和地下水流方程：

其中 r 和 z 是軸對(duì)稱圓柱坐標(biāo)系的徑向和垂直坐標(biāo)，Krr 和 Kzz 分別是沿徑向和垂直軸的主要水力傳導(dǎo)率。

當(dāng)調(diào)用軸對(duì)稱選項(xiàng)時(shí)，BCF4包將徑向坐標(biāo)（r）視為水平x坐標(biāo)，而垂直坐標(biāo)（z）不變，并且通過設(shè)置行數(shù)使水平y(tǒng)坐標(biāo)無效網(wǎng)格塊 (NROWS) 為 1。請(qǐng)注意，任何網(wǎng)格塊 ()xi , )zk) 的軸對(duì)稱配置是通過塊繞 z 軸旋轉(zhuǎn)生成的。

2.2B 變飽和水流方程的建立

水在可變飽和地下系統(tǒng)中的 3-D 運(yùn)動(dòng)由方程 (1) 表示。為了解決變飽和流動(dòng)問題，還需要指定相對(duì)滲透率與水相飽和度、壓頭與水相飽和度的關(guān)系。使用兩種替代函數(shù)表達(dá)式來描述相對(duì)滲透率與水飽和度的關(guān)系。這些函數(shù)由（Brooks 和 Corey，1966）給出：

和（van Genuchten，1977）：

當(dāng) n 和 ( 是經(jīng)驗(yàn)參數(shù)，Se 是有效含水飽和度，定義為 Se=(Sw-Swr)/(1-Swr)，Swr 被稱為殘余含水飽和度。BrooksCorey 表達(dá)式產(chǎn)生與 van 類似的曲線Genuchten 函數(shù)，當(dāng)選擇參數(shù) n 使得 n = 1 + 2/( (van Genuchten, 1980)。

按照 Mohanty 等人 (1997) 的建議，通過包含分段連續(xù)水力傳導(dǎo)函數(shù)，進(jìn)一步增強(qiáng)了模擬能力，以包含系統(tǒng)中優(yōu)先流動(dòng)路徑的雙峰描述。因此，如果系統(tǒng)中的飽和度小于臨界值S*，則使用等式(4b)的van Genuchten函數(shù)。然而，當(dāng)飽和度高于 S* 時(shí)，將應(yīng)用指數(shù)電導(dǎo)率曲線，該曲線會(huì)在接近飽和時(shí)快速增加電導(dǎo)率，如下所示：

其中 S* 是臨界或斷點(diǎn)飽和值，其中流量從毛細(xì)管主導(dǎo)變?yōu)榉敲?xì)管主導(dǎo)，反之亦然，kr * 是對(duì)應(yīng)于 S* 的相對(duì)滲透率，* 是代表有效大孔隙度或其他貢獻(xiàn)結(jié)構(gòu)特征的擬合參數(shù)到非毛細(xì)管主導(dǎo)的流動(dòng)。還提供了表格輸入選項(xiàng)，可以容納其他功能形式。

壓頭（R = h-z，其中 z = 垂直向上坐標(biāo)）與水飽和度的關(guān)系由以下函數(shù)描述（van Genuchten，1977，1980）：

其中α和β是經(jīng)驗(yàn)參數(shù)，hc是毛細(xì)管水頭，定義為(hap-Ψ)，其中hap，空氣中的壓頭取大氣壓(=0)，Swr是殘余水相飽和度。參數(shù) β 和 γ 之間的關(guān)系為 γ = 1-1/β。可以在實(shí)驗(yàn)室中測(cè)量給定土壤的保水性和相對(duì)滲透性特性的 Brooks-Corey 和 van Genuchten 函數(shù)。

BCF4 包通過 BCF4 輸入文件中的索引 IREALSL 使用實(shí)際土壤函數(shù) [方程 (4) 或 (5) 和方程 (6)] 調(diào)用可變飽和流模擬。對(duì)于 van Genuchten 土壤函數(shù)，需要額外輸入?yún)?shù) "、$ 和 Swr。如果需要 Brooks Corey 相對(duì)滲透率函數(shù)，則還需要 Brooks Corey 參數(shù) n 作為輸入。其余參數(shù)與討論的參數(shù)相同2.2 節(jié)中的 3-D、橫截面或軸對(duì)稱情況。

2.2C 變飽和空氣流量方程的建立

可變飽和多孔介質(zhì)中的空氣流量可表示為（Huyakorn 等人，1994）：

其中 xi , xj (i, j = 1, 2, 3) 是笛卡爾坐標(biāo)，koij 是固有滲透率張量，kra 是對(duì)空氣的相對(duì)滲透率，μa 是空氣的動(dòng)力粘度，Pa 是空氣相中的壓力，ρa(bǔ) 是空氣密度，Wa 是代表空氣源和/或匯的每單位體積多孔介質(zhì)的空氣體積通量，Sa (=1-Sw) 是空氣相飽和度。

在MODFLOWSURFACT/MODHMS框架下求解式(7)時(shí)，可以方便地將方程轉(zhuǎn)化為MODFLOW類型變量。因此，空氣相需要壓頭和水頭的等效定義：

其中hap是空氣相的壓力水頭，Ma是空氣相的勢(shì)能，Dw是水的密度，g是重力加速度，haM是空氣相的勢(shì)水頭或總水頭。進(jìn)一步地，根據(jù)定義，水力傳導(dǎo)率Kij可表示為： (7)

其中，hap是空氣相的壓力水頭，Φa是空氣相的勢(shì)能，ρw是水的密度，g是重力加速度，haΦ是空氣相的勢(shì)水頭或總水頭。進(jìn)一步地，根據(jù)定義，水力傳導(dǎo)率Kij可表示為：

因此，等式（7）左邊的第一項(xiàng)可以使用等式（8）和（9）來操作并重寫為

方程 (7) 中的最后一項(xiàng)可以以與 Bear (1979, p. 91) 對(duì)水相所做的類似的方式對(duì)空氣相進(jìn)行操作，得到：

其中 Ssa 是空氣的特定體積存儲(chǔ)量，定義類似：

其中 αs 是土壤基質(zhì)的壓縮性，βa 是空氣的壓縮性。對(duì)于大多數(shù)含水層系統(tǒng)來說，Ssa 不是一個(gè)有形參數(shù)，需要進(jìn)一步操作才能根據(jù) Ss（多孔材料的特定存儲(chǔ)量）提出方程，定義為（Bear，1997）

其中 βw 是水的壓縮性。等式（12）和（13）產(chǎn)量。

使用(7)中的方程(10)、(11)和(14)給出了以MODFLOW類型變量和基本參數(shù)表示的空氣相流方程：

僅求解空氣相流方程時(shí)，必須對(duì)系統(tǒng)中的水相做出假設(shè)，因?yàn)樗灰暈楸粍?dòng)相。

MODFLOW-SURFACT/MODHMS 中提供了兩個(gè)選項(xiàng)，用于在求解氣相流時(shí)處理水相。作為第一個(gè)選項(xiàng)，假定水處于靜水平衡狀態(tài)，具有恒定的水頭，等于域中任何 x、y 位置處地下水位的初始高程 ZWT (x,y)。

因此，域內(nèi)任何高度 z 處的水相壓頭 hwp（=之前使用的 Ψ）如下：

其中水壓 hwp 等于地下水位處的環(huán)境氣壓 (=0)。

因此，方程（6）中所需的毛細(xì)管水頭[hc=hap-hwp]可通過使用（19）和（8）獲得：

作為第二個(gè)選擇，可以從同一域上的水流方程的穩(wěn)態(tài)解獲得方程(6)的環(huán)境水壓。因此，通過在環(huán)境空氣條件下提供正確的空氣和水飽和度，可以在空氣相流計(jì)算期間適應(yīng)再充電/放電系統(tǒng)。

系統(tǒng)壓力條件下的空氣密度由理想氣體定律得出：

式中ρstd為大氣壓條件下的標(biāo)準(zhǔn)空氣密度，Patm；Pa為系統(tǒng)內(nèi)空氣的絕對(duì)壓力；并且假設(shè)系統(tǒng)條件下空氣的表觀分子量與標(biāo)準(zhǔn)條件下的相同。重新整理方程（8）并表達(dá)絕對(duì)壓力條件下的氣壓可得出。

其中參考?jí)毫?Patm 出現(xiàn)在參考高程 z=0（即基準(zhǔn)面）處?？諝饷芏葘?duì)系統(tǒng)勢(shì)能的影響通常很小，并且無論基準(zhǔn)面高程如何，都可以針對(duì)平均海平面條件規(guī)定 Patm，而不會(huì)顯著降低精度。

將 (22) 代入 (21) 并重新排列，以 MODFLOW-SURFACT/MODHMS 變量的形式給出空氣相密度：

總結(jié)可變飽和空氣相流的公式，MODFLOWSURFACT/MODHMS 使用方程 (15) 來控制受方程 (6) 的毛細(xì)管飽和關(guān)系和方程 (16) 或方程 (18) 的相對(duì)滲透率行為影響的氣流行為。毛細(xì)管水頭由方程（20）確定，采用水的瞬時(shí)靜水平衡的被動(dòng)水相假設(shè)，或者以環(huán)境 hwp 作為輸入（從理查茲方程的穩(wěn)態(tài)解獲得）作為 hi = hap-hwp 確定。

最后，根據(jù)方程（23）計(jì)算空氣相密度。

BCF4 包使用 BCF4 輸入文件中的索引 IREALSL 調(diào)用空氣相流公式。除了土壤濕度關(guān)系和相對(duì)滲透率函數(shù)所需的參數(shù)外，空氣相流模擬還需要輸入 Dw、μw、μa、$w、$a、Catm、g 和 ZWT (x,y)。其余的仿真輸入?yún)?shù)與 2.2 節(jié)中針對(duì) 3-D、橫截面或軸對(duì)稱情況討論的參數(shù)相同。

BCF4 包使用 BCF4 輸入文件中的索引 IREALSL 調(diào)用空氣相流公式。除了土壤濕度關(guān)系和相對(duì)滲透率函數(shù)所需的參數(shù)外，空氣相流模擬還需要輸入 ρw、μw、μa、βw、βa、Patm、g 和 ZWT (x,y)。其余的仿真輸入?yún)?shù)與 2.2 節(jié)中針對(duì) 3-D、橫截面或軸對(duì)稱情況討論的參數(shù)相同。

2.3A 無側(cè)限流驗(yàn)證示例

BCF4 包通過將數(shù)值結(jié)果與 Neuman (1972) 的解析解進(jìn)行了驗(yàn)證，該解析解針對(duì)無承壓含水層中的井流，該含水層對(duì)抽水表現(xiàn)出顯著延遲的產(chǎn)量響應(yīng)。

所考慮的問題涉及以 3658 ft3/d 的速率從均質(zhì)且各向同性的無承壓含水層中抽水的完全穿透井。含水層的初始飽和厚度為 40 英尺，井在含水層厚度較低的 30 英尺部分進(jìn)行篩選。含水層性質(zhì)如下：

Hydraulic conductivity,?K = 21 ft/d?

Specific yield, Sy = 0.132?

Specific storage, Ss = 1.2 x 10–4 ft–1

建模方法和結(jié)果

上述流動(dòng)問題可使用 BCF4 軟件包中提供的全 3D 和軸對(duì)稱仿真選項(xiàng)來解決。將獲得的數(shù)值結(jié)果與解析解進(jìn)行比較。

對(duì)于完全 3-D 模擬運(yùn)行，模型區(qū)域是一個(gè)尺寸為 2,000 英尺 x 2,000 英尺的正方形，井位于中心。為了利用對(duì)稱性，我們 (22) (23) 僅離散流域的一個(gè)象限并考慮四分之一的井流量 Q。3-D 網(wǎng)格由 20 層組成，每層 2 英尺厚，有 21 行， 21 列。水平網(wǎng)格間距 ()x 和 )y) 從 0.5 英尺到 100 英尺不等。VCONT（垂直導(dǎo)水率除以兩個(gè)相鄰節(jié)點(diǎn)層之間的層間距離）的值計(jì)算為 10.5。每層的恒定主存儲(chǔ)系數(shù)（SF1，定義為特定存儲(chǔ)時(shí)間塊厚度）被指定為 2.4 x 10–4。井流量（象限網(wǎng)格為 Q/4）均勻分布到第 6 層到第 20 層第 1 行、第 1 列的 15 個(gè)井節(jié)點(diǎn)。假設(shè)外側(cè)、底部和頂部邊界為無流動(dòng)邊界條件。強(qiáng)隱式過程 (SIP) 包用于求解矩陣方程。

對(duì)于軸對(duì)稱仿真運(yùn)行，模型域是半徑為 2,000 英尺、厚度為 40 英尺的圓柱體。該網(wǎng)格由 20 層、2 英尺厚、21 行和 1 列組成。水平（徑向）網(wǎng)格間距（△x≡△r ）從0.5英尺到100英尺不等。全井流量Q，均勻分布到第6層至第20層第1行第1列的15個(gè)井節(jié)點(diǎn)。其余信息與完整 3D 模擬運(yùn)行給出的信息相同。

對(duì)距離井 31 英尺、位于底部、中間和頂部節(jié)點(diǎn)層中心的三個(gè)選定觀測(cè)點(diǎn)進(jìn)行了水位下降與時(shí)間關(guān)系的比較。圖 2.1 中繪制的無量綱回撤 sD 和無量綱時(shí)間 ty 定義為：

其中b是初始飽和厚度，其余符號(hào)如前定義。

圖 2.1 將全 3-D 和軸對(duì)稱模擬的數(shù)值結(jié)果與 Neuman 解析解進(jìn)行了比較。對(duì)于底部和中間的觀測(cè)點(diǎn)，數(shù)值解和解析解非常一致。對(duì)于頂部觀測(cè)點(diǎn)，兩個(gè)數(shù)值解非常吻合，但偏離了早期值的解析解。這種偏差可能是由于恒定飽和厚度的簡(jiǎn)化假設(shè)和解析解所使用的近似地下水位邊界條件的影響。請(qǐng)注意，我們的數(shù)值解基于嚴(yán)格的可變飽和流公式。

2.3B 水不飽和流量驗(yàn)證示例

提供了兩個(gè)示例問題來演示 BCF4 包中非飽和區(qū)水流選項(xiàng)的驗(yàn)證和應(yīng)用。第一個(gè)問題考慮的是非常干燥的土柱中水的瞬時(shí)滲透。第二個(gè)問題考慮非飽和矩形土板中的二維流動(dòng)。將這些問題與適用于各自情況的數(shù)值解進(jìn)行比較。

2.3B.1 問題 1 - 非飽和垂直土柱中的瞬時(shí)滲透

Huang 等人 (1996) 詳細(xì)介紹了這個(gè)問題，涉及非飽和垂直土柱中的一維滲透，該土柱長(zhǎng) 160 厘米，離散為 80 層，厚 2 厘米。整個(gè)域的初始水頭為-104 cm，代表極其干燥的條件，并且在 t=0 時(shí)在頂部應(yīng)用 50 cm/d 的恒定滲透通量。

對(duì)典型的粗質(zhì)地土壤進(jìn)行了長(zhǎng)達(dá) 0.8 天的分析，其水力特性如下：

van Genuchten 函數(shù)用于保濕性和相對(duì)滲透性。

建模方法和結(jié)果

選擇 PCG4 求解器來求解收斂公差為 0.05 cm 的矩陣方程組。使用AT04包進(jìn)行自適應(yīng)時(shí)間步進(jìn)，初始△t=0.2×10-7 d，最大△t=0.01 d，時(shí)間步乘數(shù)為1.3，時(shí)間步除因數(shù)為2.7。圖 2.B1 顯示了 0.4 d 和 0.8 d 時(shí)的預(yù)測(cè)水頭剖面。這些結(jié)果與使用 Huang 等人 (1996) 和 HydroGeoLogic (1992) 的有限差分和有限元模型獲得的結(jié)果非常吻合。這個(gè)問題還證明了代碼處理高度非線性情況和極其干燥條件的能力。

2.3B.2 問題 2 - 非飽和矩形土板中的二維流動(dòng)

該問題涉及截面尺寸為水平 15 厘米和垂直 10 厘米的矩形土板中的二維非飽和流動(dòng)。使用 )x=1 cm 和 )z=1 cm 的恒定節(jié)點(diǎn)間距來離散化域。整個(gè)域的初始水頭為-90 厘米。沿著域的左邊緣的頂部4cm被規(guī)定為6cm的邊界水頭，并且域的右邊緣被維持在-90cm的規(guī)定水頭。其余邊界（即頂部邊界、底部邊界和左邊緣下 6 厘米）為無流動(dòng)條件。土壤的水力特性如下：

假設(shè)土壤是均勻且各向同性的。 van Genuchten 關(guān)系用于保濕曲線，而 Brooks Corey 關(guān)系用于定義相對(duì)滲透率函數(shù)。

建模方法和結(jié)果

選擇 PCG4 求解器來求解收斂公差為 0.01 cm 的矩陣方程組。初始時(shí)間步長(zhǎng) 0.01 天擴(kuò)大了 1.2 倍，最大時(shí)間步長(zhǎng)為 0.05 天。圖 2.B2 繪制了時(shí)間值為 0.508 天時(shí)沿流動(dòng)區(qū)域底部和頂部的水頭計(jì)算值。這些結(jié)果與使用 HydroGeoLogic (1992) 有限元模型獲得的結(jié)果非常吻合。

2.3C 氣相流模擬驗(yàn)證示例

通過將數(shù)值結(jié)果與使用綜合多相流模擬器獲得的結(jié)果進(jìn)行比較，驗(yàn)證了 BCF4 包的空氣相流模擬選項(xiàng)（Huyakorn 等人，1994 年）。

Panday 等人 (1994) 詳細(xì)介紹了這里考慮的問題。選擇這個(gè)問題來研究空氣噴射過程的流動(dòng)動(dòng)力學(xué)。本研究使用了 Marley 和 Magee (1992) 提供的相關(guān)現(xiàn)場(chǎng)信息。相關(guān)地點(diǎn)已被包氣區(qū)和飽和區(qū)幾個(gè)泄漏的地下儲(chǔ)罐的汽油污染。地下水位平均深度為地表以下-2.4 m。污染面積-450平方米。通過僅考慮可變飽和氣流（即假設(shè)被動(dòng)液相），可以對(duì)注入井附近的噴射影響半徑和空氣速度進(jìn)行初步評(píng)估。考慮 30 m × 20 m × 3 m 的矩形域，并將其離散為 900 個(gè)節(jié)點(diǎn)，其均勻節(jié)點(diǎn)間距為 )x = 2m、)y = 2m 和 )z = 0.5m。

初始條件對(duì)應(yīng)于包氣區(qū)的主要大氣壓力 (haM=0)。通過在位于 x = 9 m、y = 9 m 和 z = 0.75 m 的節(jié)點(diǎn)處以 10-3m3 s -1 的穩(wěn)定速率注入空氣來創(chuàng)建噴射過程。地下水位位于噴水井下方 0.4 m 處（即 z=0.35 處）。在整個(gè)模擬過程中，表面邊界保持大氣壓條件，底部和側(cè)面邊界被視為無流動(dòng)條件。模擬中使用的土壤和流體屬性采用 Marely 和 Magee (1992) 給出的值，如下所示。

van Genuchten 函數(shù)用于描述空氣-水系統(tǒng)的相對(duì)滲透率和毛細(xì)管壓力。

建模方法和結(jié)果

針對(duì)這個(gè)問題進(jìn)行了瞬態(tài)和穩(wěn)態(tài)仿真。

瞬態(tài)仿真調(diào)用AT04包進(jìn)行自適應(yīng)時(shí)間步進(jìn)，初始時(shí)間步長(zhǎng)為1秒，最小△t=1×10-3s，最大△t=1×105s，時(shí)間步乘數(shù)為1.4，時(shí)間步長(zhǎng)縮減因子為 2。選擇 PCG4 軟件包來求解矩陣方程組，瞬態(tài)情況下最多 20 次非線性迭代，穩(wěn)態(tài)模擬最多 100 次非線性方程。瞬態(tài)和穩(wěn)態(tài)模擬的頭部閉合標(biāo)準(zhǔn)分別為 0.1 和 0.21 m。

圖 2.C1 描述了沿穿過井的水平線的模擬氣壓分布，其中包括兩種模擬以及在 t=7 分鐘和穩(wěn)態(tài)條件下使用 Huyakorn 等人 (1994) 的有限元代碼進(jìn)行的模擬。預(yù)測(cè)和觀察到的壓力分布與這兩個(gè)規(guī)范都相當(dāng)一致。有趣的是，瞬態(tài)模擬僅在 1.16 天內(nèi)就達(dá)到了穩(wěn)態(tài)條件，顯示了空氣相中發(fā)生的快速平衡過程。這些模擬的絕對(duì)氣壓是通過等式（22）從電位頭計(jì)算的，Da是根據(jù)等式（23）計(jì)算的。

2.4 應(yīng)用示例

提供了兩個(gè)示例問題來演示新 BCF4 包的可變飽和流選項(xiàng)的應(yīng)用。在這兩個(gè)示例中，都允許發(fā)生足夠的抽水以使無承壓含水層的部分去飽和。第一個(gè)示例是穩(wěn)態(tài)泵送模擬。第二個(gè)例子是瞬態(tài)模擬，其中在四年的壓力期后關(guān)閉油井以恢復(fù)含水層。該場(chǎng)景也使用 MODFLOW 中現(xiàn)有的潤(rùn)濕選項(xiàng)（McDonald 等，1991）和 STAFF3D 有限元代碼（HydroGeoLogic, Inc., 1992）進(jìn)行建模。下面對(duì)不同方法獲得的結(jié)果進(jìn)行比較和討論。

2.4.1 問題 1——無承壓含水層抽水的穩(wěn)態(tài)模擬

該測(cè)試問題考慮了圖 2.2 所示的 300 英尺厚的無承壓含水層。

建模域是一個(gè)尺寸為 75,000 英尺 x 75,000 英尺的正方形。含水層的頂部和底部海拔分別為 50 英尺和 –250 英尺。該區(qū)域受到 3.0 x 10-9 英尺/秒的均勻、連續(xù)的垂直補(bǔ)給。圖中左（西）邊界為恒水頭邊界，其余邊界為無流邊界。含水層底部 100 英尺處有 15 口井，每口井的抽水速度為 0.95 英尺3/秒。井位如圖2.2所示。含水層參數(shù)包括：

建模方法
該域被均勻離散為 3 層、15 行、15 列的網(wǎng)格塊（即 △z = 100 ft，△x = △y = 5,000 ft）。因此，VCONT（垂直導(dǎo)水率除以兩個(gè)相鄰節(jié)點(diǎn)層之間的層間距離）的值為10-7s-1。頂層和中層西側(cè)規(guī)定恒定水頭為零。代碼所需的初始頭部分布在整個(gè)域中統(tǒng)一設(shè)置為零。 15 口井位于第 3 層。雖然這些井充當(dāng)匯，但地表補(bǔ)給和西部的恒定水頭邊界為含水層系統(tǒng)提供補(bǔ)給，從而維持所需的穩(wěn)定流量條件。

執(zhí)行兩次單獨(dú)的穩(wěn)態(tài)仿真運(yùn)行：(a) 所有三層的輸入 LAYCON 值為 3，(b) 輸入 LAYCON 值為 43。請(qǐng)注意，LAYCON 是一個(gè)控制參數(shù)，指示特定節(jié)點(diǎn)層所需的含水層處理選項(xiàng)。案例 (a) 使用原始 MODFLOW 公式來計(jì)算水平塊間透射率和電導(dǎo)，而案例 (b) 則激活 BCF4 軟件包中包含的新選項(xiàng)。 （參見BCF4包的輸入指令。）對(duì)于這兩種情況，SIP包用于求解矩陣方程。除了各層的輸入 LAYCON 值不同之外，BCF4 包和其他包（Basic、Well、Recharge、SIP 和 Output Packages）的剩余輸入數(shù)據(jù)是相同的。用于模擬的輸入數(shù)據(jù)文件列于附錄 B 中。

為了保持對(duì)流動(dòng)系統(tǒng)的連續(xù)充電，規(guī)定的充電速率應(yīng)用于每個(gè)垂直列中的最高活性單元。如果任何接受再充電的活性電池變干，這對(duì)于情況（a）來說是必要的，以允許連續(xù)滲透。原因是，如果可變水頭電池干燥并且不將補(bǔ)給傳輸?shù)降叵滤?，則可變水頭電池會(huì)轉(zhuǎn)換為無流動(dòng)電池。如果使用新的可變飽和流選項(xiàng)[情況(b)]，則情況并非如此，并且避免了數(shù)據(jù)輸入中的這種復(fù)雜性。

結(jié)果與討論

圖 2.3 至 2.5 分別對(duì)第 1、2 和 3 層由先前 MODFLOW 選項(xiàng) [案例 (a)] 和新的可變飽和流選項(xiàng) [案例 (b)] 產(chǎn)生的穩(wěn)態(tài)水頭分布進(jìn)行了比較。盡管域中的所有單元都作為可變頭輸入，但當(dāng)使用 MODFLOW [案例 (a)] 中的先前 BCF1 選項(xiàng)進(jìn)行模擬時(shí)，第 1 層中大約四分之一的單元被停用（見圖 2.3），因?yàn)槟M水頭值低于該層的底部高程（對(duì)于第 1 層為 –50 英尺）。請(qǐng)注意，可變飽和流選項(xiàng)[情況 (b)] 預(yù)測(cè)去飽和區(qū)域內(nèi)的水頭分布，從而允許將垂直補(bǔ)給傳輸?shù)降叵滤灰约皝碜韵噜弳卧膫?cè)向流分量。

第 2 層和第 3 層的穩(wěn)態(tài)水頭分布（分別如圖 2.4 和 2.5）表明，隨著距去飽和區(qū)距離的增加，兩種模擬選項(xiàng)之間的一致性會(huì)提高。在域的東南角附近，兩次模擬的頭部分布存在顯著差異。對(duì)案例 (a) 和案例 (b) 的結(jié)果進(jìn)行仔細(xì)檢查表明，BCF4 封裝的可變飽和流選項(xiàng)中使用的水平塊間電導(dǎo)表現(xiàn)適當(dāng)。原始方案[情況(a)]中使用的與水頭相關(guān)的透射率的諧波加權(quán)會(huì)導(dǎo)致更大的下降，因?yàn)楫?dāng)下降較大時(shí)，透射率會(huì)向較低的值加權(quán)。從物理角度來看，諧波加權(quán)應(yīng)應(yīng)用于飽和導(dǎo)水率，而不是與水頭相關(guān)的有效透射率。 （演示見附錄 A。）

2.4.2 問題 2——無承壓含水層抽水和采收的瞬態(tài)模擬

本例中的瞬態(tài)分析還考慮了 300 英尺厚的無承壓含水層，如圖 2.2 所示。與問題 1 一樣，建模域?qū)?yīng)于尺寸為 75,000 英尺 x 75,000 英尺的正方形。含水層的頂部和底部分別位于 50 和 –250 英尺的海拔處。整個(gè)域中含水層的初始地下水位為 0 英尺。該域受到 3.0 x 10 –9 英尺/秒的均勻、連續(xù)補(bǔ)給。圖中的左（西）邊界是恒定（零）水頭邊界，其余邊界是無流邊界。如圖 2.2 所示，底層厚度超過 100 英尺，有 15 口井被篩選。

然而，在此問題中，每口井的抽水速率高達(dá) 3.85ft3/s 。四年后，井中的抽水將被關(guān)閉。 （見圖 2.6。）然后允許含水層恢復(fù) 12 年。含水層參數(shù)為：

建模方法

與問題 1 一樣，域被均勻離散為 3 層、15 行和 15 列。 MODFLOW所需的VCONT（垂直導(dǎo)水率除以相鄰節(jié)點(diǎn)層之間的層間距離）的值計(jì)算為10–7s–1。

每層指定恒定的主存儲(chǔ)系數(shù)（SF1 定義為特定存儲(chǔ)時(shí)間塊厚度）10 –4。頂部和中間離散層的西側(cè)規(guī)定恒定水頭為零。最初，假設(shè)域中的頭為零。位于第 3 層的 15 口井在前 4 年允許以 3.85 的速率抽水。在隨后的 12 年期間停止抽水，以便恢復(fù)含水層。

執(zhí)行三個(gè)獨(dú)立的瞬態(tài)仿真運(yùn)行。前兩個(gè)模擬分別利用 (1) 現(xiàn)有的 MODFLOW 再潤(rùn)濕功能和 (2) MODFLOW-SURFACT/MODHMS 的新可變飽和流（BCF4 封裝）選項(xiàng)。這兩個(gè)模擬分別稱為情況（1）和（2）。第三次模擬運(yùn)行 [案例 (3)] 使用有限元 STAFF3D 代碼獨(dú)立執(zhí)行（HydroGeoLogic, Inc., 1992）。對(duì)于基于 MODFLOW 的仿真運(yùn)行，SIP 包用于求解矩陣方程。除了 BCF4 包輸入之外，兩種情況使用的其余輸入數(shù)據(jù)是相同的。請(qǐng)注意，情況 (1) 需要額外的數(shù)據(jù)才能重新潤(rùn)濕干電池。

為了在 MODFLOW 再潤(rùn)濕包 [案例 (1)] 模擬選項(xiàng)中保持連續(xù)充電，表面充電應(yīng)用于每個(gè)垂直列中最高的活性電池（請(qǐng)參閱充電包輸入）。這一規(guī)定是必要的，以允許接收補(bǔ)給的活性單元的連續(xù)滲透（如果有的話）變得干燥，因?yàn)榭勺冾^單元如果干燥則轉(zhuǎn)變?yōu)榉橇鲃?dòng)單元，并且不將補(bǔ)給傳送到地下水位。

結(jié)果與討論

對(duì)于使用原始再潤(rùn)濕能力 [案例 (1)] 和新的可變飽和流選項(xiàng) [案例 (2)] 執(zhí)行的模擬運(yùn)行，第 1 層和第 3 層（代表兩種極端飽和條件）的水頭分布如圖 2.7 所示到 2.10。圖 2.7 比較了抽水周期結(jié)束時(shí)（t=4 年）第一層的水頭分布。對(duì)于情況 (1)，第 1 層中大約四分之一的細(xì)胞因 4 年的連續(xù)泵送而失活，但第 3 層仍然保持完全活躍。

請(qǐng)注意，較長(zhǎng)的泵浦周期會(huì)導(dǎo)致第 3 層去飽和并引發(fā)數(shù)值解的不穩(wěn)定。圖 2.8a 和 2.8b 比較了時(shí)間值為 4 年時(shí)第 3 層的頭部分布。第 3 層中沒有失活的單元。案例 (1) 和 (2) 的水頭值之間存在顯著差異（高達(dá) 20 至 30 英尺），發(fā)生在各個(gè)井附近。遠(yuǎn)離井場(chǎng)的兩個(gè)水頭分布非常吻合。圖 2.9 和 2.10 分別顯示了恢復(fù)期結(jié)束時(shí)（t=16 年）第 1 層和第 3 層的水頭比較。一般來說，MODFLOW 的先前（再潤(rùn)濕）選項(xiàng)會(huì)低估水頭值。兩個(gè)水頭分布之間的差異在流域的東南部更為明顯。

對(duì)三個(gè)選定的觀察點(diǎn)進(jìn)行了水頭與時(shí)間關(guān)系的比較。下面參考圖 2.2 給出了這些點(diǎn)的位置。

三個(gè)觀察點(diǎn)處水頭的時(shí)間變化（圖 2.11）表明，兩種情況 [情況 (1) 和 (2)] 非常一致，直到情況 (1) 的方案將可變水頭單元轉(zhuǎn)換為無流動(dòng)單元。經(jīng)過這些轉(zhuǎn)換后，案例（1）的結(jié)果與案例（2）的結(jié)果有偏差，但遵循相似的趨勢(shì)。如前所述，BCF2 包 [案例 (1)] 的臨時(shí)再潤(rùn)濕程序可能無法產(chǎn)生可靠的數(shù)值解。例如，如果在瞬態(tài)模擬的第一階段，井的抽水速率從 3.85 增加到 3.90ft3/s，MODFLOW 的再潤(rùn)濕選項(xiàng)將無法收斂。然而，新的可變飽和流選項(xiàng)會(huì)產(chǎn)生收斂解。 BCF2 包中使用的再潤(rùn)濕程序失敗是由于迭代之間電池的循環(huán)干燥和潤(rùn)濕（翻轉(zhuǎn)）產(chǎn)生的不穩(wěn)定性。

使用 MODFLOW-SURFACT/MODHMS 可變飽和流選項(xiàng) [案例 (2)] 獲得的結(jié)果也與使用 STAFF3D 代碼獲得的結(jié)果（稱為案例 (3)）進(jìn)行了比較。圖 2.12 顯示了三個(gè)觀察點(diǎn)的頭部時(shí)間變化圖。在整個(gè)模擬期間，情況（2）和（3）之間有很好的一致性。

2.5 以塊為中心的流包的其他擴(kuò)展

對(duì)以塊為中心的流 (BCF4) 包進(jìn)行了一些細(xì)微的增強(qiáng)，以允許模擬具有更大的靈活性。

第一個(gè)修改允許輸入網(wǎng)格塊的垂直水力傳導(dǎo)率，作為輸入相鄰層網(wǎng)格塊之間的泄漏的選項(xiàng)。 BCF4 封裝輸入的第一行中包含了一個(gè)標(biāo)志 IVHYC，以表明這一點(diǎn)。

如果 IVHYC 設(shè)置為 1，則代碼需要輸入所有層的垂直水力電導(dǎo)率，并且根據(jù)每個(gè)網(wǎng)格塊的垂直電導(dǎo)率及其間隔距離的體積加權(quán)調(diào)和平均值在內(nèi)部計(jì)算泄漏。如果 IVHYC 設(shè)置為零，則代碼需要輸入層之間的泄漏，如原始 MODFLOW 中一樣。

第二個(gè)修改允許逐個(gè)單元地輸入水平各向異性，作為每個(gè)模型層具有統(tǒng)一的各向異性值的選項(xiàng)。 BCF4 包輸入的第一行中包含一個(gè)標(biāo)志 IANIXY，以表明這一點(diǎn)。如果 IANIXY 設(shè)置為 1，則代碼需要逐個(gè)單元輸入各向異性。如果 IANIXY 設(shè)置為零，則代碼期望像原始 MODFLOW 中那樣逐層輸入各向異性。

第三種修改允許輸入曲線網(wǎng)格，作為使用矩形網(wǎng)格在區(qū)域平面中離散化域的選項(xiàng)。 BCF4 包輸入的第一行中包含一個(gè)標(biāo)志 ICURVL，以表明這一點(diǎn)。如果 ICURVL 設(shè)置為 1，則使用曲線網(wǎng)格，并且代碼期望為層內(nèi)的所有節(jié)點(diǎn)輸入 DELR 和 DELC。由于層彼此垂直堆疊，因此所有層的 DELR 和 DELC 值都是根據(jù)第一層的值在內(nèi)部計(jì)算的。如果 ICURVL 設(shè)置為零，則像原始 MODFLOW 中一樣使用矩形網(wǎng)格。對(duì)于此 MODHMS 情況，僅輸入 DELR 的列數(shù) (NCOL)，僅輸入 DELC 的網(wǎng)格中的行數(shù) (NROW)。 MODHMS 為非飽和區(qū)水流模擬提供了額外的靈活性，允許通過重力分離垂直平衡 (GSVE) 偽土壤函數(shù)提供的受限/非受限地下水流模擬功能來表示更深的層。此功能對(duì)于非飽和區(qū)域流量模擬非常有用，可以避免地下水位附近的過度垂直離散，這需要準(zhǔn)確表示無承壓水位條件和相關(guān)的含水層流量。因此，當(dāng) IREALSL=1 或 2（不飽和/飽和水流模擬）時(shí)，可以設(shè)置頂部幾層來求解 LAYCON 值為 43 的理查德方程，而下面較厚的層可用于表示無承壓/承壓地下水系統(tǒng)。這是通過在 MODHMS 中為這些較低層設(shè)置 LAYCON=44 來實(shí)現(xiàn)的，這些層將被視為 IREALSL=0。

2.6 BCF4 封裝的輸入指令

2.6.1 概述

BCF4 包所需的輸入數(shù)據(jù)與之前的 BCF 包類似。 BCF4 包采用了新的可變飽和流公式，通過使用偽土壤函數(shù)，通過去飽和和再飽和來嚴(yán)格處理無側(cè)限流問題。出于比較目的，我們還將所有先前記錄的方案包含在新包中。

因此，我們的 BCF4 封裝包括永久停用干電池的 BCF1 無限制封裝（McDonald 和 Harbaugh，1988）、將無流動(dòng)電池轉(zhuǎn)換為可變頭電池的 BCF2 封裝（McDonald 等人，1991）以及 BCF3 封裝，它使用塊間透射率關(guān)系的不同選項(xiàng)來模擬地下水流（Goode 和 Appel，1992）。建議有興趣比較新舊方案的用戶熟悉各自的 BCF 包文檔。

BCF4 的輸入數(shù)據(jù)主要由屬性數(shù)組組成，用于描述水平水力傳導(dǎo)率/透射率、垂直泄漏、單元的頂部和底部標(biāo)高以及瞬態(tài)模擬的存儲(chǔ)系數(shù)和/或比產(chǎn)量。這些屬性數(shù)組取決于每層的輸入 LAYCON 值，如下所述。 BCF2 和 BCF3 無限制存儲(chǔ)方案所需的附加參數(shù)數(shù)組是 WETDRY，當(dāng)該選項(xiàng)的潤(rùn)濕標(biāo)志處于活動(dòng)狀態(tài)時(shí)使用該數(shù)組。 BCF4 內(nèi)存分配包在讀取某些標(biāo)志后為這些數(shù)組保留空間。

然后根據(jù)用戶輸入初始化這些數(shù)組，如下所述。以下頁面提供的輸入說明適用于 BCF4 包。請(qǐng)注意，BCF4 的設(shè)計(jì)旨在保持與其前身的整體兼容性。

要激活新的可變飽和流選項(xiàng)，用戶將層的輸入 LAYCON 指定為 40（如果層受到嚴(yán)格限制）或 43（否則）。嚴(yán)格限制層定義為在整個(gè)模擬期間測(cè)壓頭保持在該層頂部標(biāo)高之上的層。如果沒有先驗(yàn)知識(shí)可用于層中的水頭條件，建議用戶使用輸入 LAYCON 43。BCF4 中實(shí)施的新無側(cè)限流公式消除了輸入準(zhǔn)備中的主要尷尬以及先前單元失活中固有的困難和再潤(rùn)濕計(jì)劃。我們建議用戶改用這種新方法以獲得更可靠的收斂模擬。

BCF4 包的 LAYCON=40 或 43 方案使用相鄰網(wǎng)格塊的水力傳導(dǎo)率的加權(quán)諧波平均來計(jì)算單元之間的電導(dǎo)。通過設(shè)置 LAYCON=50 或 53，可以選擇使用相鄰網(wǎng)格塊的水力傳導(dǎo)率的加權(quán)算術(shù)平均來計(jì)算單元之間的電導(dǎo)率（類似于有限元裝配的執(zhí)行方式），作為額外的靈活性。

2.6.2 BCF4 輸入

BCF4 包的輸入是從 IUNIT (1) 中指定的單元讀取的。下一小節(jié)將提供輸入變量的簡(jiǎn)要說明。數(shù)據(jù)在模擬開始時(shí)讀取一次。由基本包讀取的 IUNIT 數(shù)組確定 MODFLOW 模擬中使用的選項(xiàng)或包。 （有關(guān)詳細(xì)信息，請(qǐng)參閱原始 MODFLOW 文檔的基本包。）請(qǐng)注意，BCF4 所需的額外輸入變量顯示在陰影框中。

以下二維數(shù)組的子集用于描述模型網(wǎng)格每一層的屬性。每層所需的陣列取決于每層的輸入 LAYCON 值（參見第 2 項(xiàng)）、模擬是瞬態(tài)（ISS = 0）還是穩(wěn)態(tài)（ISS 不為 0），以及潤(rùn)濕功能是否處于活動(dòng)狀態(tài)（IWDFLG非 0) 對(duì)于 BCF2 和 BCF3 的潤(rùn)濕方案。請(qǐng)注意，對(duì)于 BCF4 封裝中引入的新的可變飽和流量方案，IWDFLG 應(yīng)設(shè)置為零（即輸入 LAYCON=43）。

如果不需要數(shù)組，則必須將其省略。不存在需要所有數(shù)組的情況。首先讀取第 1 層所需的數(shù)組（第 8-15 項(xiàng)）；然后是第 2 層的數(shù)組

限制為 200 層。如果層數(shù)在40層以下，則使用1條記錄；否則，使用兩個(gè)或多個(gè)記錄。將記錄中未使用的元素留空。 LAYCON 是可以實(shí)現(xiàn) BCF4 封裝新選項(xiàng)的參數(shù)，如下所述。

BFC4 封裝支持 20 個(gè)不同的輸入 LAYCON 值。這些值以 I2 格式輸入。根據(jù)該值，第一個(gè)數(shù)字（十位）（如果有）根據(jù)該數(shù)字的值確定塊間透射率函數(shù)類型的應(yīng)用。輸入 LAYCON 值的第二個(gè)數(shù)字（個(gè)）決定了原始模型中的層特征。第二個(gè)數(shù)字的含義如下（McDonald and Harbaugh，1988）：

上述描述可以總結(jié)如下（Goode and Appel，1992）：

請(qǐng)注意，BCF4 包的新可變飽和流選項(xiàng)是通過模型網(wǎng)格中所有層的輸入 LAYCON = 40、43、50、53 或 54 來實(shí)現(xiàn)的。此外，LAYCON=44 或 54 的值可用于在其他不飽和區(qū)域水流模擬（即 IREALSL …0）中將層表示為無約束/約束層。因此，當(dāng)輸入 LAYCON 值的第二位為 4 時(shí)，其含義與上面討論的 3 相同。此新選項(xiàng)需要 IAPART 的非零值。標(biāo)志 IAPART 從 BAS 包的第 5 項(xiàng)中讀取，指示數(shù)組 BUFF 是否與數(shù)組 RHS 分開。

由輸入 LAYCON 的第一個(gè)數(shù)字表示的因子 10 存儲(chǔ)在 LAYAVG 數(shù)組中。第二個(gè)數(shù)字（個(gè)）存儲(chǔ)在 LAYCON 數(shù)組中。下表給出了輸入 LAYCON 值的 LAYCON 和 LAYAVG 的存儲(chǔ)值：

因此，輸入 LAYCON 值的十位部分存儲(chǔ)在數(shù)組 LAYAVG 中，輸入 LAYCON 值的個(gè)位部分存儲(chǔ)在 LAYCON 數(shù)組中。 LAYAVG 的存儲(chǔ)值決定了計(jì)算塊間透射率的方法（當(dāng) LAYAVG = 10、20、30、40 或 50 時(shí)）。

TRPY – 是一個(gè)一維數(shù)組，包含每層的各向異性因子。它是沿列的透射率或水力傳導(dǎo)率（以所使用的為準(zhǔn)）與沿行的透射率或水力傳導(dǎo)率的比率。對(duì)于各向同性條件，設(shè)置為 1.0。這是一個(gè)單數(shù)組，每層有一個(gè)值。不要為每一層讀取一個(gè)數(shù)組；整個(gè)數(shù)組僅包含一個(gè)數(shù)組控制記錄。如果所有層的值都相同，則可以通過將 LOCAT 設(shè)置為 0 并將 CNSTNT 設(shè)置為適用于所有層的值來指定整個(gè)數(shù)組。 LOCAT 和 CNSTNT 的解釋見下文。

注意：軸對(duì)稱仿真不需要 TRPY 數(shù)組 (IAXIS>0)。

TRPYNU – 是一個(gè)二維數(shù)組，包含每層內(nèi)逐個(gè)單元的各向異性因子。它是沿列的透射率或水力傳導(dǎo)率（以所使用的為準(zhǔn)）與沿行的透射率或水力傳導(dǎo)率的比率。對(duì)于各向同性條件，設(shè)置為 1.0。為每一層讀取一個(gè)數(shù)組。

注意：軸對(duì)稱仿真不需要 TRPYNU 數(shù)組 (IAXIS>0)。

LOCAT – 是填充相關(guān)數(shù)組的值位置的索引。

<0，讀取包含數(shù)組值的未格式化記錄 = 0，將所有數(shù)組值設(shè)置為等于常量 (CNSTNT) >0，讀取包含數(shù)組值的格式化記錄。

CNSTNT – 是一個(gè)常量，如果 LOCAT 等于 0，則所有數(shù)組值都將設(shè)置為該常量；如果 LOCAT 不等于 0，則所有數(shù)組值都會(huì)乘以該常量。

DELR – 是沿行的單元格寬度。讀取每個(gè) NCOL 列的一個(gè)值。這是一個(gè)數(shù)組，每一列都有一個(gè)值。對(duì)于軸對(duì)稱模擬 (IAXIS>0)，DELR 表示圓柱形電池的徑向?qū)挾?()r)。

DELRCV – 是曲線網(wǎng)格沿行的單元寬度。這是一個(gè)二維數(shù)組，對(duì)于一層的每個(gè) NROW 行，每個(gè) NCOL 列包含一個(gè)值。

DELC – 是沿列的單元格寬度。為每一 NROW 行讀取一個(gè)值。這是一個(gè)數(shù)組，每一行都有一個(gè)值。

注意：對(duì)于軸對(duì)稱模擬 (IAXIS>0)，不會(huì)讀取 DELC 數(shù)組。

DELCCV – 是曲線網(wǎng)格沿列的單元寬度。這是一個(gè)二維數(shù)組，其中包含一個(gè)層的每個(gè) NROW 行的每個(gè) NCOL 列的值。

ZWT – 是用作氣流計(jì)算參考的域中的地下水位高程。

PWC – 水相的水頭，用作氣流計(jì)算的參考。

SF1——初級(jí)存儲(chǔ)系數(shù)，可存儲(chǔ)性。對(duì)于瞬態(tài)模擬只讀（穩(wěn)態(tài)標(biāo)志 ISS 為 0）。對(duì)于輸入 LAYCON 等于 1、11、21 或 31，SF1 將始終為具體良率，而對(duì)于輸入 LAYCON 等于 2、12、22、3、13、23、33、40 或 43，SF1 將始終為受限存儲(chǔ)系數(shù)。對(duì)于輸入 LAYCON 等于 0、10 或 20，SF1 通常是受限存儲(chǔ)系數(shù)。主存儲(chǔ)系數(shù)除以塊厚度等于等式(1)和(2)的特定存儲(chǔ)(Ss)。

注 (1)：輸入 LAYCON 值 0、10 或 20 也可用于模擬地下水位條件，其中水位下降預(yù)計(jì)在任何地方都保持為飽和厚度的一小部分，并且上面沒有層或流量從上面看可以忽略不計(jì)。在這種情況下，將為 SF1 輸入特定的屈服值。

TRAN – 是沿網(wǎng)格行的透射率。 TRAN 乘以 TRPY 即可獲得沿列的透射率。對(duì)于由以下輸入 LAYCON 值之一表示的層，TRAN 是只讀的：0、2、10、12、20、22。

HY – 是沿網(wǎng)格行的水力傳導(dǎo)率。 HY 乘以 TRPY 即可獲得沿柱的水力傳導(dǎo)率。對(duì)于由以下輸入 LAYCON 值之一表示的層，HY 是只讀的：1、3、11、13、21、23、31、33、40、43。

VHY – 是每個(gè)單元的垂直水力傳導(dǎo)率。對(duì)于由以下輸入 LAYCON 值之一表示的層，VHY 是只讀的：1、3、11、13、21、23、31、33、40、43。

BOT——含水層底部的高程。對(duì)于由以下輸入 LAYCON 值之一表示的層只讀：1、3、11、13、21、23、31、33、40、43。

VCONT——垂直導(dǎo)水率除以從一層到下一層的厚度。單元的值是水力傳導(dǎo)率除以該單元中的節(jié)點(diǎn)與下面單元中的節(jié)點(diǎn)之間的材料的厚度。因?yàn)榈讓酉旅鏇]有層，所以沒有為底層指定VCONT。

SF2——是二次儲(chǔ)存系數(shù)，即比產(chǎn)量。如果模擬是瞬態(tài)的（穩(wěn)態(tài)標(biāo)志 ISS 為 0）并且該層的輸入 LAYCON 值為以下之一，則 SF2 為只讀：2、3、12、13、22、23、33、43。

請(qǐng)注意，SF2 是求解非飽和區(qū)域中的氣流或理查德方程時(shí)的孔隙率。

TOP——含水層頂部的標(biāo)高。對(duì)于由以下輸入 LAYCON 值之一表示的層，TOP 是只讀的：2、3、12、13、22、23、33、40、43。

WETDRY——是潤(rùn)濕閾值和標(biāo)志的組合，用于指示哪些相鄰單元可以導(dǎo)致單元變濕。如果 WETDRY < 0，則只有干電池下方的電池才會(huì)導(dǎo)致電池變濕。如果 WETDRY > 0，干電池下方的電池和四個(gè)水平相鄰的電池可能會(huì)導(dǎo)致電池變濕。如果 WETDRY 為 0，則電池不能被潤(rùn)濕。 WETDRY 的絕對(duì)值是潤(rùn)濕閾值。當(dāng)干電池的 BOT 和 WETDRY 的絕對(duì)值之和等于或超過相鄰電池的水頭時(shí)，該電池被潤(rùn)濕。如果 IWDFLG 不為 0 并且輸入 LAYCON 值為以下之一，則 WETDRY 為只讀：1、3、11、13、21、23、31、33。

VANAL - 是非飽和土壤的 van Genuchten 參數(shù)。僅當(dāng)圖層的輸入 LAYCON 值為 43 且使用 IREALSL…0 請(qǐng)求真實(shí)土壤濕度函數(shù)時(shí)，VANAL 才為只讀。

VANBT——是非飽和土的van Genuchten參數(shù)$。當(dāng)圖層的輸入 LAYCON 值為 43 且使用 IREALSL…0 請(qǐng)求真實(shí)土壤濕度函數(shù)時(shí)，VANBT 為只讀。

VANSR——是非飽和土的殘余飽和水平。當(dāng)圖層的輸入 LAYCON 值為 43 且使用 IREALSL…0 請(qǐng)求真實(shí)土壤濕度函數(shù)時(shí)，VANSR 為只讀。

BROOK——是計(jì)算非飽和土相對(duì)滲透率的 Brooks-Corey 指數(shù)。當(dāng)圖層的輸入 LAYCON 值為 43 時(shí)，BROOK 為只讀，且實(shí)際土壤濕度函數(shù)與 BrooksCorey 滲透率關(guān)系要求 IREALSL=2 或 4。

S_STAR – 是具有 Mohanty 等人 (1997) 的雙峰電導(dǎo)率函數(shù)的非飽和土壤的臨界飽和參數(shù) S*（參見第 2.2B 節(jié)的方程 5b）。當(dāng)圖層的輸入 LAYCON 值為 43 且 IREALSL=7 時(shí)請(qǐng)求雙峰電導(dǎo)率函數(shù)時(shí)，S_STAR 為只讀。

D_STAR – 是使用 Mohanty 等人 (1997) 的生物模態(tài)電導(dǎo)率函數(shù)時(shí)，非飽和土壤的有效大孔隙度擬合參數(shù)（參見第 2.2B 節(jié)的方程 5b）。當(dāng)圖層的輸入 LAYCON 值為 43 且 IREALSL=7 時(shí)請(qǐng)求雙峰電導(dǎo)率函數(shù)時(shí)，D_STAR 為只讀。

H_VAL – 是保濕和相對(duì)滲透率函數(shù)表格輸入插值點(diǎn)的壓頭值。當(dāng)層的輸入 LAYCON 值為 43 且 IREALSL=5 或 6 時(shí)請(qǐng)求函數(shù)的表格輸入時(shí)，H_VAL 為只讀。

S_VAL – 是與 H_VAL 相關(guān)值的插值點(diǎn)對(duì)應(yīng)的含水飽和度值。當(dāng)層的輸入 LAYCON 值為 43 且 IREALSL=5 或 6 時(shí)請(qǐng)求函數(shù)的表格輸入時(shí)，S_VAL 為只讀。

K_VAL – 是活動(dòng)流動(dòng)相的相應(yīng)相對(duì)滲透率值，對(duì)于具有相關(guān) H_VAL 值的插值點(diǎn)。

當(dāng)層的輸入 LAYCON 值為 43 且 IREALSL=5 或 6 時(shí)請(qǐng)求函數(shù)的表格輸入時(shí)，K_VAL 為只讀。

2.6.3 輸入指令中使用的字段說明

第 3 章：裂縫井 (FWL4) 套件

3.1 概述 ?

HydroGeoLogic, Inc. 在美國(guó)地質(zhì)調(diào)查局模塊化地下水流模型 MODFLOW 中添加了新的井包，可將井筒模擬為高導(dǎo)率裂縫管。這些包是 MODFLOW 中現(xiàn)有井包（WEL1 包）的替代方案，被稱為 Fracture-Well（FWL4 和 FWL5）包。 Fracture-Well 封裝旨在克服與原始 WEL1 封裝相關(guān)的多個(gè)問題。

首先，F(xiàn)WL4 和 FWL5 封裝模擬通過多層孔的流動(dòng)。這是因?yàn)榱芽p管表示允許將與井相關(guān)的含水層單元（節(jié)點(diǎn)）與一維有限直徑圓柱形井單元連接。此外，還嚴(yán)格考慮了該井規(guī)定的總提取率。來自與壓裂井相關(guān)的每個(gè)單獨(dú)節(jié)點(diǎn)的體積通量由代碼自動(dòng)計(jì)算，以求和到從井中提取的總流量。另一方面，原始的 WEL1 包僅在代表井的含水層節(jié)點(diǎn)處實(shí)現(xiàn)源/匯項(xiàng)。

因此，它要求用戶確定每個(gè)節(jié)點(diǎn)對(duì)穿透兩個(gè)或多個(gè)節(jié)點(diǎn)的物理井的貢獻(xiàn)。通常估計(jì)這些節(jié)點(diǎn)通量貢獻(xiàn)，以便根據(jù)貢獻(xiàn)井節(jié)點(diǎn)的透射率來分配規(guī)定的提取率。這種估計(jì)本身可能是一個(gè)嚴(yán)重的錯(cuò)誤，從而歪曲了物理系統(tǒng)。此外，對(duì)于無約束系統(tǒng)，根據(jù) WEL1 包規(guī)定節(jié)點(diǎn)通量可能會(huì)導(dǎo)致計(jì)算困難。然而，F(xiàn)WL4 和 FWL5 軟件包會(huì)自動(dòng)將凈井提取量（由用戶輸入）分配給井節(jié)點(diǎn)。對(duì)于每個(gè)時(shí)間步長(zhǎng)，計(jì)算出的節(jié)點(diǎn)通量可能會(huì)有很大變化，就像各個(gè)層去飽和和再飽和時(shí)的情況一樣。

新FWL4和FWL5套件的另一個(gè)特點(diǎn)是，對(duì)于超抽無側(cè)限系統(tǒng)，當(dāng)井中水位達(dá)到井底時(shí)，總井抽水量會(huì)自動(dòng)調(diào)整。在此階段，底部井節(jié)點(diǎn)的水頭固定在井底條件下，并計(jì)算總抽水量，以對(duì)應(yīng)于該井在這些條件下可以供應(yīng)的可行數(shù)量。如果發(fā)生抽水需求小于井在井底壓力條件下可供應(yīng)的情況（即，在特定時(shí)期減少抽水或增加補(bǔ)給），井節(jié)點(diǎn)處的水頭將反彈為了反映這些新的情況，該井將提供全部的開采需求。當(dāng)水位低于井底時(shí)，以前的 WEL1 套件不會(huì)物理調(diào)整抽水速率。相反，它會(huì)繼續(xù)計(jì)算到非物理領(lǐng)域，井節(jié)點(diǎn)處的水頭值不斷下降到井底標(biāo)高以下。在其他情況下，WEL1 包會(huì)在迭代之間循環(huán)干燥和潤(rùn)濕泵單元，從而導(dǎo)致系統(tǒng)在計(jì)算上變得不穩(wěn)定。

FWL4 和 FWL5 軟件包還實(shí)現(xiàn)了井眼存儲(chǔ)，這在 MODFLOW 的早期版本中被忽略。 FWL4 和 FWL5 軟件包之間的差異在于，F(xiàn)WL4 軟件包假設(shè)疊加，將井流方程添加到多孔基質(zhì)的方程中，而 FWL5 軟件包使用對(duì)數(shù)井函數(shù)來描述從基質(zhì)塊到井節(jié)點(diǎn)的流動(dòng)。因此，如果用戶有興趣使用 FWL4 軟件包準(zhǔn)確預(yù)測(cè)局部水位下降，則仍然需要在每口井附近進(jìn)行精細(xì)的面積離散化。

最后，F(xiàn)WL4 和 FWL5 封裝僅適用于 BCF4 封裝的嚴(yán)格可變飽和模擬選項(xiàng)（即輸入 LAYCON = 40 或 43）。

本章的以下部分進(jìn)一步闡述了 FWL4 和 FWL5 包的功能，并包括對(duì)公式和演示示例問題的簡(jiǎn)要討論。壓裂井公式的數(shù)值實(shí)現(xiàn)也根據(jù)文獻(xiàn)中可用的分析解決方案進(jìn)行了驗(yàn)證，然后是各個(gè)包的輸入說明。

3.2 裂縫-井方程的公式化

通過將代表有限半徑井眼篩選部分的 1-D 線元素疊加到代表多孔介質(zhì)的 3-D 單元上，為 FWL4 包制定壓裂井，如圖 3.1 所示。 Sudicky 等人（1995）詳細(xì)介紹了壓裂井單元的推導(dǎo)和使用，并提供了示范性示例。出于一般性的考慮，我們將他們的技術(shù)應(yīng)用于 BCF4 包中實(shí)現(xiàn)的可變飽和流量公式。該公式進(jìn)一步擴(kuò)展到 FWL5 包，其中使用井函數(shù)在多孔基質(zhì)塊和井節(jié)點(diǎn)之間發(fā)送流量，而不是使用疊加假設(shè)。這將在本節(jié)稍后討論。 FWL4 和 FWL5 套件中通過管狀裂縫的流量均由一維可變飽和流量方程控制，如下所示：

式中，l為沿井篩的距離（L），K′l為井的導(dǎo)水率（LT-1），k′rw為井內(nèi)的相對(duì)滲透率，S′w為井的無量綱飽和厚度。井單元（飽和厚度與單元厚度之比），S's 為井眼特定儲(chǔ)存量（L–1），h' 為井內(nèi)水頭（L）。

偽土壤函數(shù)用于定義壓裂井中的不飽和流動(dòng)行為，如 BCF4 包文檔（第 2 章）中所述，以便井單元內(nèi)的地下水位條件由垂直平衡條件準(zhǔn)確表示。井眼的飽和水力傳導(dǎo)率，根據(jù)層流理論，應(yīng)用 Hagen-Poiseuille 方程估算，可表示為 (Sudicky et al, 1995)：

其中D是水的密度(ML–3)，g是重力加速度，:是水的粘度(ML–1 T–1)，rw是井的半徑(L)。通常，對(duì)于大直徑井，井水力傳導(dǎo)率值比多孔基體的塊間水力傳導(dǎo)率大幾個(gè)數(shù)量級(jí)。因此，仿真結(jié)果對(duì)使用等式(25)可能產(chǎn)生的誤差不敏感。

如前所述，F(xiàn)WL4 包使用井節(jié)點(diǎn)與多孔基體節(jié)點(diǎn)的疊加。另一方面，F(xiàn)WL5 軟件包將它們保持為單獨(dú)的節(jié)點(diǎn)，并使用 Theim 方程井函數(shù)提供矩陣塊與其關(guān)聯(lián)井節(jié)點(diǎn)之間的流動(dòng)。將 Bennett 等人 (1982) 提出的概念應(yīng)用到任意層 k 的裂縫井公式中，基質(zhì)塊與其井節(jié)點(diǎn)之間的通量由下式給出：?

其中 Qmwk 是 k 層井通量矩陣； krk 為相對(duì)滲透率，Kk 和 Thickk 為基質(zhì)塊的水平電導(dǎo)和厚度，hmk 為 k 層基質(zhì)塊中的水頭，hwk 為 k 層井節(jié)點(diǎn)的水頭，Sk 為表皮系數(shù)re 是矩陣塊的有效半徑，定義為（Halford 和 Hanson，2002；Neville 和 Tonkin，2004）。

其中 △x 和 △y 是包含孔的單元格的水平網(wǎng)格塊尺寸。

然后，F(xiàn)WL5 程序包會(huì)同時(shí)求解井中的水頭和多孔基質(zhì)塊節(jié)點(diǎn)，這些節(jié)點(diǎn)彼此不同；因此，多孔基質(zhì)網(wǎng)格塊不需要在井周圍進(jìn)行細(xì)化，以獲得更好的井流量、水頭和提取條件的準(zhǔn)確性。

FWL4 和 FWL5 軟件包要求用戶提供的信息包括：井位、特定井篩選的頂層和最底層的標(biāo)識(shí)、井底高程、井半徑、抽取 (-ve) 或注入 (+ve) ）來自井的速率（L3T–1），以及井筒特定儲(chǔ)存值（L–1）。此外，F(xiàn)WL5 軟件包還需要輸入表面摩擦系數(shù)和在任何應(yīng)力周期激活任何井時(shí)的井內(nèi)初始水頭。 FWL4和FWL5模塊組織幾何信息，使用方程(25)計(jì)算Kl，并將離散裂縫井方程(24)應(yīng)用到各自的井節(jié)點(diǎn)。然后，在每個(gè)應(yīng)力周期，將抽取或注入應(yīng)用于井中最底部的篩選層。該公式將應(yīng)力適當(dāng)?shù)胤植荚诰?jié)點(diǎn)之間。最后，當(dāng)需要時(shí)，井節(jié)點(diǎn)處的邊界條件（即，從規(guī)定的通量到規(guī)定的水頭）的切換自動(dòng)執(zhí)行，如下所述。

每口井按規(guī)定的速率繼續(xù)抽取，直至井內(nèi)水位達(dá)到井底。井的底部標(biāo)高實(shí)際上位于被篩選的最低節(jié)點(diǎn)層中的任何位置。當(dāng)最低井節(jié)點(diǎn)處的測(cè)壓水頭對(duì)應(yīng)于井的底部高程時(shí)，水頭值保持固定，并且允許井產(chǎn)量低于需求并達(dá)到零。如果在后期，含水層系統(tǒng)恢復(fù)并且井產(chǎn)量增加到大于需求，或者如果需求下降到產(chǎn)量以下，則井恢復(fù)以規(guī)定的抽水速率抽水。

輸出文件中報(bào)告了總井產(chǎn)量和每層對(duì)井產(chǎn)量的貢獻(xiàn)。在以下部分中，提供說明性示例來驗(yàn)證和演示 FWL4 和 FWL5 軟件包的應(yīng)用。

3.3 驗(yàn)證示例

通過將數(shù)值結(jié)果與 Theis (1935) 以及 Papadopulos 和 Cooper (1967) 獲得的解析解進(jìn)行比較，驗(yàn)證了壓裂井包 (FWL4) 選項(xiàng)，該解析解針對(duì)具有恒定流量井的均質(zhì)含水層中的水位下降。與 Theis 解不同，Papadopulos 和 Cooper 的解考慮了井筒儲(chǔ)存。

因此，進(jìn)行了兩種模擬，一種有井眼存儲(chǔ)，另一種沒有井眼存儲(chǔ)，以測(cè)試 FWL4 中使用的數(shù)值方案。

Sudicky 等人 (1995) 描述了驗(yàn)證問題。在此問題中，半徑為 0.1 m 的全穿透井位于尺寸為 1,000 m x 1,000 m x 5 m 的均質(zhì)各向同性含水層的中心。該井的抽水速度為 10 –3 m3 /s。

橫向邊界規(guī)定了 20 m 的恒定水頭，并且假定含水層的頂部和底部邊界是不滲透的。含水層的水力傳導(dǎo)率為10–4 m/s，比容為10–4 m–1。使用 FWL4 軟件包執(zhí)行瞬態(tài)仿真。

建模方法和結(jié)果

該域被離散為 5 層，每層 1 英尺厚，具有 35 行和 35 列可變網(wǎng)格間距（從井附近的 111.1 英尺到 0.1 英尺）。含水層的頂部和底部高程分別設(shè)置為 0 和 –5 m。 MODFLOW 要求的 VCONT 值（垂直導(dǎo)水率除以任意兩層之間的層間距離）計(jì)算為 10–4 s–1。含水層的主要蓄水系數(shù)（具體蓄水時(shí)間塊厚度）指定為10 – 4。含水層側(cè)面規(guī)定恒定水頭20 m，并假設(shè)初始水頭均勻?yàn)?0 m。 0.1 m 半徑井位于區(qū)域中心（第 18 行、第 18 列），允許以 10–3 m3 /s 的速率泵送 1000 秒。由于井已完全貫通，井頂和井底分別位于1層和5層，井底標(biāo)高指定為–5 m。為了將井壓降的時(shí)間變化與 Theis 和 Papadopulos 的解析解進(jìn)行比較，在特定井筒儲(chǔ)水量分別為零和 0.1 m-1 的情況下進(jìn)行了兩次瞬態(tài)模擬。對(duì)于這兩種模擬，PCG4 軟件包用于求解矩陣方程。圖 3.2 繪制了井面的壓降與時(shí)間曲線以及 Theis 和 Papadopulos 解析解。如圖所示，使用FWL4 Package獲得的MODFLOWSURFACT/MODHMS結(jié)果與各自的解析解非常吻合，從而驗(yàn)證了FWL4 Package中采用的壓裂井公式。

3.4 應(yīng)用示例

提供了兩個(gè)示例來演示FWL4包的應(yīng)用，第三個(gè)示例演示了FWL5包的驗(yàn)證和應(yīng)用。第一個(gè)示例說明了完全 3D 情況，考慮了從無承壓含水層抽水的井。第二個(gè)示例說明了準(zhǔn) 3-D 情況，考慮了從由限制床分隔的三個(gè)含水層中抽水的井。在這兩種情況下，井都被故意過度抽水，并在第二個(gè)壓力期恢復(fù)含水層。作為對(duì)井邊界條件實(shí)施嚴(yán)格性的全面檢查，系統(tǒng)的響應(yīng)也使用不調(diào)用 FWL4 包的獨(dú)立模擬進(jìn)行評(píng)估。對(duì)于全 3-D 情況，從壓裂井模擬獲得的井處水頭值被強(qiáng)加為驗(yàn)證模擬的時(shí)間相關(guān)水頭條件，并將結(jié)果與從原始模擬獲得的結(jié)果進(jìn)行比較，其中 FWL4 包被使用了。最后一個(gè)示例說明了使用帶有準(zhǔn) 3-D 設(shè)置的 FWL5 軟件包從由約束床分隔的兩個(gè)含水層進(jìn)行井抽水。該問題的受限情況與解析解進(jìn)行比較進(jìn)行驗(yàn)證，該問題的非受限情況則展示了包裝的干燥能力。

3.4.1 問題 1——無承壓含水層抽水的全三維模擬

該問題考慮了圖 3.3 所示的 300 英尺厚的均質(zhì)、各向同性無承壓含水層。建模域是一個(gè) 75,000 英尺 x 75,000 英尺的正方形，含水層的頂部和底部高程分別為 +50 英尺和 –250 英尺。整個(gè)區(qū)域含水層的初始地下水位為零。圖中左（西）邊界為恒水頭邊界，其余邊界為無流邊界。一口直徑為 0.4 英尺的全穿透井和 10 –6 英尺 –1 的井眼特定存儲(chǔ)位于該區(qū)域的中心。 10 ft 3 /s 容量的泵運(yùn)行 2.5 年后，停止從井中抽水。 （見圖 3.4。）然后允許含水層恢復(fù) 4 年。含水層參數(shù)包括：

建模方法

該域被離散為 3 層、15 行和 15 列（△x = △y = 5,000 ft 和 △z = 100 ft）。 VCONT（垂直導(dǎo)水率除以兩個(gè)相鄰層之間的層間距離）的值計(jì)算為 10 –6 s–1。為每層指定 10-4 的恒定主存儲(chǔ)系數(shù)（特定存儲(chǔ)時(shí)間塊厚度）。

頂部和中間離散層的西側(cè)規(guī)定恒定水頭為零。

直徑 0.4 英尺的井位于區(qū)域中心（第 8 行、第 8 列），并在所有 3 層中進(jìn)行屏蔽直至含水層底部。前 2.5 年的提款需求為 10 ft 3 /s。請(qǐng)注意，當(dāng)含水層無法提供指定的抽水速率（即，當(dāng)井中的水位到達(dá)井底時(shí)）時(shí)，F(xiàn)racture-Well (FWL4) 包會(huì)自動(dòng)調(diào)整抽水速率。在第二個(gè)壓力期，抽水會(huì)停止 4 年，以便含水層恢復(fù)。由于井已完全貫通，井的頂部和底部分別位于第 1 層和第 3 層，井底標(biāo)高指定為 –250 英尺。

使用 BCF4（輸入 LAYCON=43 的可變飽和流選項(xiàng)）和 FWL4 封裝執(zhí)行瞬態(tài)仿真。 PCG4 軟件包用于求解矩陣方程。用于模擬的輸入數(shù)據(jù)文件列在附錄 C 中。

結(jié)果與討論

井內(nèi)水位和總井通量隨時(shí)間的變化分別如圖 3.5 和圖 3.6 所示。如圖 3.5 所示，井內(nèi)水位在 2.15 年內(nèi)從初始值（0）逐漸下降至井底（-250 英尺）。水位保持在井底直至停止抽水。在第二個(gè)（恢復(fù)）期間，井中水位從-250 英尺上升到-17.13 英尺?？偩康南鄳?yīng)變化可以在圖3.6 中觀察到。在第一個(gè)應(yīng)力期間，從井中的總抽油量從 10 英尺 3 /秒減少到 7.63 英尺 3 /秒。當(dāng)抽水停止時(shí)，總井通量減少至零。

采用圖 3.5 中描述的瞬態(tài)水頭分布，使用上述相同的水文條件執(zhí)行規(guī)定的井口模擬，但不調(diào)用 FWL4 包。本次仿真的目的是驗(yàn)證之前使用 FWL4 封裝進(jìn)行仿真時(shí)獲得的結(jié)果的準(zhǔn)確性。從之前的模擬中獲得的瞬態(tài)水頭值是針對(duì)無承壓含水層第 3 層、第 8 行和第 8 列（即底部井節(jié)點(diǎn)）的單元指定的（圖 3.3）。

為了適應(yīng)規(guī)定的瞬態(tài)水頭分布，單元被視為具有大界面水力傳導(dǎo)率（109 ft/s）的一般水頭邊界節(jié)點(diǎn)。

圖 3.6 比較了驗(yàn)證井模擬和壓裂井模擬運(yùn)行獲得的模擬結(jié)果。兩組結(jié)果非常吻合，從而驗(yàn)證了 FWL4 套件對(duì)規(guī)定抽水量井進(jìn)行嚴(yán)格處理的準(zhǔn)確性。

3.4.2 問題2——從三含水層系統(tǒng)中抽取的準(zhǔn)三維模擬

該問題涉及由兩個(gè) 10 英尺厚的零滲透性限制層（含水層）分隔的三個(gè)含水層，如圖 3.7 所示。請(qǐng)注意，使用非零泄漏可以輕松處理泄漏的弱透水層。建模域是一個(gè)尺寸為 75,000 英尺的正方形。三個(gè)含水層的頂部和底部高程如下：

整個(gè)域內(nèi)所有含水層的初始測(cè)壓水頭均為零。圖中左（西）邊界為零水頭常水頭邊界，其余邊界為無流邊界。一口直徑為 0.4 英尺的全穿透井和 10 –6 英尺 –1 的井眼特定儲(chǔ)存量位于該區(qū)域的中心。以 5 ft 3 /s 的速度運(yùn)行三年后，將停止從井中抽水（圖 3.8）。然后允許含水層恢復(fù) 9 年。三個(gè)含水層的水力參數(shù)如下：

建模方法

域中的每個(gè)含水層都被視為節(jié)點(diǎn)層，并且限制層（含水層）中的垂直滲漏為零。每個(gè)含水層均勻分為 15 行和 15 列 ()x = )y = 5,000 ft)。第 1 層、第 2 層和第 3 層分別厚 90 英尺、40 英尺和 30 英尺。通過將 VCONT（含水層之間的垂直滲漏）的值指定為零，可以納入限制層對(duì)流量的影響。每個(gè)含水層指定恒定的初級(jí)存儲(chǔ)系數(shù)（特定存儲(chǔ)時(shí)間塊厚度）0.001。在頂部和中間離散層的西側(cè)規(guī)定了恒定的零水頭。一口直徑為 0.4 英尺的井位于域的中心（第 8 行、第 8 列），前 3 年允許以 5 英尺 3 /秒的速率抽水。請(qǐng)注意，裂縫井包會(huì)自動(dòng)將總抽取量分配給各個(gè)含水層。如果總需求不可行（即井內(nèi)水位達(dá)到井底時(shí)），則將井回抽率調(diào)整至含水層系統(tǒng)能夠供給的值。在隨后的 9 年內(nèi)停止抽水以恢復(fù)含水層。由于該井已完全貫通，井的頂部和底部分別位于第 1 層和第 3 層，井底標(biāo)高指定為 -130 ft。井筒具體儲(chǔ)存量為 10–6 ft–1。

使用 BCF4（輸入 LAYCON=43 的可變飽和流選項(xiàng)）和 FWL4 封裝執(zhí)行瞬態(tài)仿真。 PCG4 軟件包用于求解矩陣方程。除了總井通量之外，F(xiàn)WL4 軟件包還打印每個(gè)時(shí)間步長(zhǎng)每個(gè)井節(jié)點(diǎn)的通量和水頭值。

結(jié)果與討論

井節(jié)點(diǎn)水頭隨時(shí)間的變化如圖 3.9 所示。每個(gè)含水層對(duì)井的通量貢獻(xiàn)如圖 3.10 所示。最初，1 號(hào)含水層的貢獻(xiàn)占井取水總量的 90% 以上（圖 3.10）。當(dāng)井節(jié)點(diǎn)的水頭下降到含水層底部標(biāo)高時(shí)，這種貢獻(xiàn)就會(huì)下降。

然后，含水層 2 和含水層 3 開始對(duì)所施加的抽水壓力產(chǎn)生更多影響。

當(dāng)井中的水位到達(dá)井底時(shí)（t = 2.3 年），井的總產(chǎn)量開始降至低于 5 ft 3 /s 的指定需求，并在 3 年時(shí)降至 1.3 ft 3 /s。當(dāng)在第二個(gè)脅迫期（3年后）停止抽水時(shí)，含水層1和含水層2繼續(xù)向井供水，而含水層3從井取水。

圖 3.11 至 3.14 分別顯示了 2.0、3.0、4.5 和 7.5 年模擬的井節(jié)點(diǎn)處測(cè)壓頭的橫截面剖面。三個(gè)含水層中的測(cè)壓頭顯示為沿第 8 行穿過井的東西向垂直剖面（見圖 3.7）。如圖 3.11 所示，在 t = 2 年時(shí)，頂部、中部和底部含水層的井測(cè)壓水位分別為 –39.98 英尺、–89.97 英尺和 –116.00 英尺。圖 3.12 表示第一個(gè)應(yīng)力期結(jié)束時(shí)含水層的狀態(tài)。停止抽水一年半后（t=4.5 年），由于通過井眼從兩個(gè)上覆含水層獲得補(bǔ)給，觀察到含水層 3 顯著恢復(fù)（圖 3.13）。圖 3.14 顯示了第二個(gè)應(yīng)力期（t=7.5 年）中期的測(cè)壓水頭，表明井中水位從 –130 英尺的井底水平上升至 –58.25 英尺。隨著測(cè)壓水頭的增加水平，含水層 3 轉(zhuǎn)換為完全承壓含水層。含水層 2 由來自含水層 1 的水通過井以及從恒水頭邊界向東補(bǔ)給。在第二個(gè)應(yīng)力期結(jié)束時(shí)（t=12 年），井中水位上升至 –23.21 英尺，含水層 2 完全受到限制（另見圖 3.9）。

3.4.3 問題 3——使用 FWL5 軟件包從兩含水層系統(tǒng)中抽取的準(zhǔn)三維模擬

該模擬的問題設(shè)置改編自 Neville 和 Tonkin (2004)，并考慮了跨兩個(gè)含水層的開放井。該區(qū)域空中呈正方形，邊長(zhǎng)為 47,000 英尺，中間有一個(gè)抽水井，每天抽水 62,840 英尺。沿域的橫向邊緣應(yīng)用無流邊界條件。該問題考慮了兩種情況——第一種是受限模擬情況，與 Papadopulos (1966) 的解析解進(jìn)行比較以進(jìn)行驗(yàn)證；第二種是無限制模擬情況，演示了在良好干燥條件下的應(yīng)用。

兩種情況下含水層的頂部和底部高程如下。

頂部含水層的初始水頭統(tǒng)一設(shè)置為 10 英尺。底部含水層的初始水頭統(tǒng)一設(shè)置為 30 英尺。頂部含水層的水力傳導(dǎo)率和比儲(chǔ)存值為 K=100 ft/d，S=0.0001，并且對(duì)于底部含水層，K=400 ft/d，S=0.0001。該井完全穿透兩個(gè)含水層，井眼半徑為 0.5 英尺。

建模方法

使用 FWL4 和 FWL5 軟件包解決了這個(gè)問題。域中的每個(gè)含水層都被視為單獨(dú)的節(jié)點(diǎn)層，通過限制床的垂直泄漏為零。 FWL5 模擬使用 100 列和 100 行的均勻間隔網(wǎng)格，網(wǎng)格間距為 470 英尺。井中的初始水頭與相鄰層中的初始水頭相同。

FWL4 模擬使用可變間距網(wǎng)格，鄰近井的網(wǎng)格間距更細(xì)。第一種情況（受限）模擬持續(xù) 1.2 天，以排除稍后橫向邊界效應(yīng)可能干擾解析解結(jié)果的比較。第二種情況（無側(cè)限）模擬持續(xù)160天，以使井內(nèi)的頭達(dá)到井底條件。

結(jié)果與討論

第一個(gè)（封閉）模擬情況下井內(nèi)水位和通量的時(shí)間變化分別如圖 3.15 和 3.16 所示。 FWL4 和 FWL5 案例的結(jié)果與解析解相匹配，從而驗(yàn)證了 FWL4 和 FWL5 包。

請(qǐng)注意，圖中的比較是使用 MNW 包進(jìn)行的模擬，它與解析解的效果非常好（Neville 和 Tonkin，2004 年）。此外，F(xiàn)WL4 包通過井周圍的精細(xì)網(wǎng)格提供了兩個(gè)含水層之間正確的水頭和通量分配。圖 3.17 和 3.18 顯示了第二個(gè)（無側(cè)限）模擬情況下井內(nèi)水頭和通量的時(shí)間變化。 FWL4 和 FWL5 包的結(jié)果非常匹配，表明具有不同井近似值的兩個(gè)包提供了相似的結(jié)果。此外，兩種包裝都能成功地應(yīng)對(duì)良好的干燥條件。

3.5 FWL4 和 FWL5 包的輸入指令

3.5.1 概述

壓裂井（FWL4 和 FWL5）包是井（WEL1）包的替代品，模擬高導(dǎo)管狀壓裂元件的井。 FWL4 和 FWL5 軟件包通過將與井相關(guān)的所有節(jié)點(diǎn)與具有高電導(dǎo)率的一維有限直徑圓柱形井單元相連來模擬井。從概念上講，壓裂井是通過將代表井眼篩選部分的這些井元素疊加到代表多孔基質(zhì)的 3-D 單元上來在 FWL4 包中制定的。從概念上講，壓裂井在 FWL5 包中被公式化為代表井眼篩選部分的單獨(dú)井元素，這些井元素根據(jù)井函數(shù)（蒂姆方程）從代表多孔基質(zhì)的相關(guān) 3-D 單元接收流量。與 WEL1 封裝相比，F(xiàn)WL4 和 FWL5 封裝的優(yōu)點(diǎn)包括：

? 代碼自動(dòng)計(jì)算與壓裂井相關(guān)的每個(gè)單獨(dú)節(jié)點(diǎn)的抽取量，以求和到用戶指定的井的總抽取量。因此，無需事先估計(jì)每個(gè)節(jié)點(diǎn)對(duì)穿透兩個(gè)或多個(gè)節(jié)點(diǎn)的物理井的貢獻(xiàn)。

? 對(duì)于超抽含水層系統(tǒng)，井中的水位保持在井底水平，井總抽水量自動(dòng)調(diào)整為井在該階段可以供應(yīng)的數(shù)量。對(duì)于多個(gè)含水層，每個(gè)含水層對(duì)油井的貢獻(xiàn)會(huì)自動(dòng)調(diào)整為壓力過大的系統(tǒng)的產(chǎn)量。如果出現(xiàn)抽水需求小于含水層在井底壓力條件下可產(chǎn)出的情況（即，在特定抽水周期內(nèi)減少井抽水量或增加補(bǔ)給量），則井中的水位會(huì)因此而測(cè)壓井節(jié)點(diǎn)處的水頭將反彈以反映新的應(yīng)力。 FWL4 和 FWL5 套件的這一功能甚至對(duì)于單層系統(tǒng)或完全從單層抽水的系統(tǒng)也很有用。

? FWL4 和FWL5 套件用于井眼存儲(chǔ)。應(yīng)該注意的是，在抽水區(qū)域仍然需要對(duì)模擬域進(jìn)行更精細(xì)的面積離散，以準(zhǔn)確地描述 FWL4 程序包井附近的水位下降，但 FWL5 程序包不需要這樣做，因?yàn)檫B接網(wǎng)格的對(duì)數(shù)函數(shù)FWL4 配方中未合并孔細(xì)胞塊，而是使用疊加。

FWL4 封裝只能與 BCF4 封裝的新 LAYCON 選項(xiàng)一起使用（即輸入 LAYCON = 40 或 43）。

3.5.2FWL4輸入

Fracture-Well (FWL4) 包的輸入是從 IUNIT (21) 中指定的單元讀取的。要調(diào)用 FWL4 包，模型網(wǎng)格中的所有層必須使用輸入 LAYCON 值 40 或 43。

FWL4 包的輸入數(shù)據(jù)主要包括井位、每口井被屏蔽的頂層和底層、井的底部高程和半徑、抽取 (-ve) 或注入 (+ve) 速率以及井眼具體信息存儲(chǔ)（L–1）。每口井的信息均在給定壓力期間的一條記錄中指定。應(yīng)力期內(nèi)的記錄數(shù)量取決于域中的井?dāng)?shù)。

MXFWEL——模擬中的最大井?dāng)?shù)。其中包括未堵塞的井，這些井可能不會(huì)抽水，但仍充當(dāng)層之間的垂直連接。

IFWLCB——是一個(gè)標(biāo)志和一個(gè)單元號(hào)。

如果 IFWLCB > 0，則它是每當(dāng)設(shè)置 ICBCFL（請(qǐng)參閱輸出控制）時(shí)將記錄孔的逐個(gè)單元流動(dòng)項(xiàng)的單元編號(hào)。

如果 IFWLCB = 0，則不會(huì)打印或記錄孔的逐個(gè)細(xì)胞流動(dòng)項(xiàng)。

如果 IFWLCB < 0，則每當(dāng)設(shè)置 ICBCFL 時(shí)都會(huì)打印井注入速率。

（如果 IFBLCB=-1，每個(gè) FWL 單元的通量貢獻(xiàn)也會(huì)被打印。）

KR = FCONST *（半徑）2

其中 FCONST =ρg/8μ（參見壓裂井公式的方程 2）

半徑 = 井篩半徑 [L]；?

ρ = 水的密度 [M/L3 ]；?

g = 重力加速度 [L/T 2 ]；和μ=水的動(dòng)力粘度[M/LT]。

FCONST 的值取決于準(zhǔn)備輸入數(shù)據(jù)集時(shí)使用的長(zhǎng)度和時(shí)間單位。下表列出了不同長(zhǎng)度和時(shí)間單位的 FCONST 值：

FCONT 尺寸值[1/LT]

如果仿真中使用的長(zhǎng)度和時(shí)間單位與上表中給出的不同，則必須使用上面給出的關(guān)系計(jì)算 FCONST 的等效值。

IFREWL4——是以自由格式讀取FWL4文件的數(shù)據(jù)項(xiàng)2和3的標(biāo)志。

如果 IFREWL4 = 0，則按照輸入指令中所述的默認(rèn)格式讀取 FWL4 文件。

如果 IFREWL4 = 1，則以自由格式讀取 FWL4 文件的數(shù)據(jù)項(xiàng) 2 和 3。

MUTWL4——是一個(gè)禁止在輸出文件中打印 FWL4 輸入的標(biāo)志。

如果 MUTWL4 = 0，則將每個(gè)應(yīng)力周期的 FWL4 輸入讀取打印到輸出文件。

如果 MUTWL4 = 1，則抑制在每個(gè)應(yīng)力周期將 FWL4 輸入打印到輸出文件，但仍會(huì)打印警告。

ITMP——是給定壓力期內(nèi)的井?dāng)?shù)（包括所有未堵塞的井）。

注：除非鉆探新井或堵塞舊井，否則 ITMP 保持不變且等于 MXFWEL。換句話說，未放電/補(bǔ)給的現(xiàn)有井必須包含在每個(gè)應(yīng)力周期的數(shù)據(jù)中，并將其 Q 值設(shè)置為零。

LSTART——是包含孔篩頂部的模型單元的層號(hào)。

LEND——是包含井底的模型單元的層號(hào)。

行——是包含孔的模型單元的行號(hào)。

列——是包含孔的模型單元的列號(hào)。

高程——是井底的高程。從物理上講，該標(biāo)高應(yīng)位于井的底層 (LEND) 內(nèi)。

半徑——是井的半徑，其單位與模擬中使用的任何其他長(zhǎng)度單位相同。將井的半徑設(shè)置為零可以有效地?cái)嚅_井裂縫單元，因?yàn)榱芽p井電導(dǎo)率（方程（1）和（2）中的K R ）變?yōu)榱恪?（由于井體積為零，井的儲(chǔ)存量也變?yōu)榱悖?/strong>

Q——是進(jìn)入井中的體積通量（L3 T–1）。正值表示井中注入水，負(fù)值表示井中排出水。

儲(chǔ)存量——井眼特定儲(chǔ)存量（S′s），由以下公式給出：S′s=（1/Ls）/（rc2/rw2），假設(shè)含水層飽和模擬，其中Ls為篩網(wǎng)總長(zhǎng)度，rc為直立水管（或套管）半徑，rw為篩網(wǎng)半徑。請(qǐng)注意，井眼儲(chǔ)存的影響對(duì)于大直徑井來說相當(dāng)大，特別是在模擬的早期階段。通過將存儲(chǔ)期限設(shè)置為零，可以忽略井眼存儲(chǔ)效應(yīng)。對(duì)于無側(cè)限模擬，因井去飽和而產(chǎn)生的井筒存儲(chǔ)自動(dòng)納入方程（4）中的井。

3.5.3 輸入指令中使用的字段說明

?3.5.4 斷裂井 (FWL5) 封裝輸入

FWL5 包的輸入是從 IUNIT (37) 中指定的單位讀取的，按照流量模擬的要求（在文件中默認(rèn)擴(kuò)展名 WL5）。首先提供所有井的輸入井屬性，包括井位置、每口井被屏蔽的頂層和底層、底部高程、井半徑和井眼特定存儲(chǔ)。

在模擬開始時(shí)，所有孔均設(shè)置為非活動(dòng)狀態(tài)。然后，對(duì)于每個(gè)應(yīng)激期，可以激活或隨后停用任意數(shù)量的井。對(duì)于在任何壓力時(shí)期激活的每口井，都會(huì)尋求井的抽水速率的附加信息（抽水為負(fù)，注入為正）、表示集膚效應(yīng)電導(dǎo)的包含/讀數(shù)的標(biāo)志以及指示輸入的標(biāo)志井中的初始水頭條件。如果這些標(biāo)志處于活動(dòng)狀態(tài)，則需要額外的相應(yīng)輸入。請(qǐng)注意，每當(dāng)將非活動(dòng)狀態(tài)激活時(shí)，都需要輸入井中的初始水頭，以防止非物理值 (HDRY) 表達(dá)井的壓力狀態(tài)。 FWL5 軟件包執(zhí)行這些檢查，如果檢測(cè)到在沒有關(guān)聯(lián)初始水頭輸入的情況下激活了井，則將中止模擬，并在輸出列表文件中顯示適當(dāng)?shù)南ⅰ４送?，為每個(gè)物種的傳輸模擬指定了注入流體的濃度。請(qǐng)注意，傳輸所需的額外輸入變量顯示在陰影框中。

請(qǐng)注意，LSTART 和 LEND 是井 LWELL 的起始層和結(jié)束層，并通過數(shù)據(jù)項(xiàng) 2 為所有井定義。

** 僅當(dāng) INHD = 1 時(shí)才輸入第 6 項(xiàng)（該孔的 ITMP 和 IBOUNDF5 也必須大于零）。

6. 數(shù)據(jù)： 格式：HDFWL5 (K), K=LSTART, LEND 8F10.3 ** 僅當(dāng) INHD = 1 且 ITRAN 0 時(shí)才輸入第 7 項(xiàng)。（對(duì)于 ≠ 此孔，ITMP 和 IBOUNDF5 也必須大于零）。

7. 數(shù)據(jù)： 格式：CONFWL5 (K, LSPEC), K=LSTART, LEND 8F10.3 注：對(duì)于模擬的 LSPEC = 1、NSPECI 污染物種類，應(yīng)重復(fù)第 7 項(xiàng)（對(duì)于每個(gè)種類）。

3.5.5 FWL5 輸入指令中使用的字段說明

MXFWEL5——模擬中的最大井?dāng)?shù)。其中包括可能不活躍（匯或源）但仍充當(dāng)層之間垂直連接的井。

IFWL5CB——是流量模擬的標(biāo)志和單元號(hào)。

如果 IFWL5CB > 0，則它是每當(dāng)設(shè)置 ICBCFL（請(qǐng)參閱輸出控制）時(shí)將記錄孔的逐個(gè)單元流動(dòng)項(xiàng)的單元編號(hào)。文件的默認(rèn)擴(kuò)展名為 CW5。

如果 IFWL5CB = 0，則不會(huì)打印或記錄孔的逐個(gè)細(xì)胞流動(dòng)項(xiàng)。

如果 IFWL5CB < 0，則每當(dāng)設(shè)置 ICBCFL 時(shí)都會(huì)打印井注入率。

（如果 IFWL5CB=-1，每個(gè) FWL 單元的通量貢獻(xiàn)也會(huì)被打印。）

FCONST5——用于計(jì)算井的水力傳導(dǎo)率的乘數(shù)，計(jì)算公式為：?

?KR = FCONST5 *（半徑）2

其中 FCONST = ρg/8μ

半徑 = 井篩半徑 [L]；

?ρ = 水的密度 [M/L3 ]；

g = 重力加速度 [L/T 2 ]；和

μ= 水的動(dòng)力粘度[M/LT] 。

FCONST5 的值取決于準(zhǔn)備輸入數(shù)據(jù)集時(shí)使用的長(zhǎng)度和時(shí)間單位。下表列出了不同長(zhǎng)度和時(shí)間單位的 FCONST5 值：

FCONT5 Value in dimensions[1/LT]

FCONT5 尺寸值[1/LT]

第 4 章：補(bǔ)給滲流面邊界 (RSF4) 條件包

4.1 概述

MODFLOW 代碼使用 Recharge (RCH1) 包適應(yīng)地下水補(bǔ)給邊界條件。然而，RCH1 包在表示無約束系統(tǒng)方面有一定的局限性。如果無承壓含水層飽和到地表，含水層吸收補(bǔ)給的能力就會(huì)降低，剩余的水會(huì)作為地表徑流流失。 RCH1 包無法處理這種情況。相反，它繼續(xù)為含水層提供補(bǔ)給（如在封閉系統(tǒng)中），水頭不斷上升至遠(yuǎn)高于地面。 MODFLOW-SURFACT 的新 RSF4 軟件包旨在消除這種物理上不切實(shí)際的情況。如果地下水位低于用戶規(guī)定的水池（積水）高度，RSF4 套件允許將所提供的補(bǔ)給到地下水系統(tǒng)中。

水池標(biāo)高對(duì)應(yīng)于無積水時(shí)的地面標(biāo)高。

如果在任何時(shí)候地下水位達(dá)到水池高度，模擬只允許進(jìn)行盡可能多的補(bǔ)給以維持規(guī)定的水池條件。剩余充值不被接受到模擬域中。模擬的輸出反映了由于系統(tǒng)飽和至積水高度而導(dǎo)致系統(tǒng)體積補(bǔ)給通量 (L3 /T) 的減少。請(qǐng)注意，RSF4 包與 Streamflow-Routing 包（STR1 或 SFR1）分離，因此，該剩余補(bǔ)給不會(huì)以任何方式路由到規(guī)定的流中。

此外，對(duì)于 MODHMS 模擬，無側(cè)限滲流選項(xiàng)與該域處于活動(dòng)狀態(tài)的地表水流域相沖突，因此應(yīng)在這些位置設(shè)置非常高的積水標(biāo)高，以防止該包與地表水流的干擾，如果無側(cè)限滲流滲流選項(xiàng)用于 RSF4 包中。

新 RSF4 軟件包的一項(xiàng)附加功能可以方便地輸入瞬態(tài)補(bǔ)給/降水，而與 MODFLOW 的“應(yīng)力周期”概念無關(guān)。使用此選項(xiàng)時(shí)，可以通過單獨(dú)的文件輸入充值時(shí)間序列，該文件確定域上的充值如何隨時(shí)間變化，而不管壓力周期設(shè)置如何，以提供輸入便利，并且通常在充值時(shí)提供較小的輸入文件大小變化很快，并且經(jīng)常在模擬周期內(nèi)變化。

除了地表補(bǔ)給-排水條件外，RSF4 軟件包還可用于通過規(guī)定滲流面邊界節(jié)點(diǎn)的高程來模擬滲流面邊界條件。因此，邊界處出現(xiàn)零通量，直到地下水位達(dá)到滲流面高程。如果地下水位達(dá)到滲流面高程，則含水層排水以維持滲流面高程（即大氣壓力）條件。如果在模擬過程中地下水位下降到滲流面高程條件以下，則滲流面邊界處將保持零通量。

RSF4 軟件包的無側(cè)限滲流選項(xiàng)僅適用于 BCF4 軟件包的嚴(yán)格可變飽和模擬選項(xiàng)（即輸入 LAYCON = 40 或 43）。下面介紹 RSF4 包的功能，包括對(duì)公式和演示示例問題的簡(jiǎn)要討論。

4.2 RSF4 邊界條件包的制定

4.3 快速、瞬態(tài)充電條件驗(yàn)證和應(yīng)用示例

?4.4 驗(yàn)證和應(yīng)用示例

?4.4.1 問題 1——無承壓含水層中平行排水溝的流量

4.4.2 問題2——方形路堤滲漏

4.5 Rsf4 封裝輸入指令

4.5.1 RSF4 輸入

4.5.2 瞬態(tài)充電時(shí)間序列文件

第 5 章：自適應(yīng)時(shí)間步進(jìn)和輸出控制 (ATO4) 包

5.1 概述

在 MODFLOW 代碼中，時(shí)域（用于瞬態(tài)流模擬）使用后向差分公式進(jìn)行離散化，時(shí)間步值通過指定以下參數(shù)預(yù)先確定：時(shí)間步數(shù) (NSTP)、應(yīng)力周期長(zhǎng)度（PERLEN），以及使用幾何級(jí)數(shù)增加時(shí)間步長(zhǎng)的乘法因子（TSMULT）。下面列出了使用這種方案的限制。

? 如果矩陣求解器或非線性迭代方案（對(duì)于無約束情況）未能在給定時(shí)間步內(nèi)收斂，則計(jì)算將中止。

? 即使獲得了收斂的數(shù)值結(jié)果，預(yù)選的時(shí)間步長(zhǎng)對(duì)于整體解決方案來說可能效率低下。

? 必須保持時(shí)間步長(zhǎng)的幾何級(jí)數(shù)。結(jié)果，用戶在分配特定時(shí)間值來檢查流動(dòng)系統(tǒng)方面缺乏靈活性。

因此，只能在時(shí)間級(jí)數(shù)上的時(shí)間值請(qǐng)求輸出，這需要由用戶提前確定。

為了克服上述困難并提高解決方案的效率，開發(fā)了自適應(yīng)時(shí)間步進(jìn)和輸出控制包（ATO4）。自適應(yīng)時(shí)間步長(zhǎng)方案根據(jù)給定計(jì)算的系統(tǒng)的預(yù)期非線性來選擇時(shí)間步長(zhǎng)大小。如果預(yù)期的非線性不顯著，則選擇較大的時(shí)間步長(zhǎng)以積極地推進(jìn)模擬。如果預(yù)期的非線性很嚴(yán)重，則選擇較小的時(shí)間步長(zhǎng)以確保該時(shí)間步長(zhǎng)的收斂。如果解在給定的時(shí)間步長(zhǎng)內(nèi)未能收斂，則進(jìn)一步減小時(shí)間步長(zhǎng)，并重復(fù)求解。

對(duì)于自動(dòng)時(shí)間步進(jìn)方案，模擬結(jié)果的時(shí)間值是未知的；因此，ATO4 軟件包還包括用于輸出控制的新模塊，該模塊確定所需輸出的時(shí)間值。調(diào)用 ATO4 包后，不再需要使用 MODFLOW 之前的輸出控制 (OC) 包。請(qǐng)注意，原始 MODFLOW 的輸出控制實(shí)現(xiàn)起來很繁瑣且麻煩。另一方面，新 (ATO4) 方案的輸出控制則簡(jiǎn)單明了。用戶只需輸入任何給定壓力期間所需輸出的時(shí)間值（以先前使用的時(shí)間單位），以及所請(qǐng)求輸出的詳細(xì)信息（即，頭、回撤、質(zhì)量預(yù)算和預(yù)算中的一項(xiàng)或全部）。細(xì)胞間的流動(dòng)項(xiàng)是否被輸出）。 ATO4 包自動(dòng)識(shí)別輸出的選定時(shí)間值，并調(diào)整時(shí)間步長(zhǎng)以在這些打印時(shí)間執(zhí)行計(jì)算。因此，例如，用戶可以在 15 個(gè)月的壓力期模擬的 10、11 和 15 個(gè)月后檢查結(jié)果。 ATO4 包不僅在整個(gè)模擬過程中提供高效的時(shí)間步進(jìn)，而且還確保在 10、11 和 15 個(gè)月（請(qǐng)求輸出的時(shí)間值）執(zhí)行計(jì)算。

最后，ATO 包提供了從非零起始時(shí)間值重新啟動(dòng)仿真的選項(xiàng)。 MODFLOW 的 OC 包不允許重新啟動(dòng)，并且開發(fā)新的輸入文件來執(zhí)行重新啟動(dòng)可能會(huì)變得極其繁瑣，特別是在涉及多個(gè)壓力周期時(shí)。

5.2 ATO4 包中使用的自適應(yīng)時(shí)間步進(jìn)的制定

5.3 驗(yàn)證示例

5.4 ATO4 封裝的輸入指令

5.4.1 概述

5.4.2 ATO4 輸入

5.4.3 輸入指令中使用的字段說明

第 6 章：預(yù)條件共軛梯度 (PCG4) 包

6.1 概述

MODFLOW 代碼生成一個(gè)以有限差分形式描述地下水流系統(tǒng)的方程組。這組代數(shù)方程通常通過使用切片連續(xù)過松弛 (SSOR) 或強(qiáng)隱式過程 (SIP) 包來求解。對(duì)于復(fù)雜的現(xiàn)場(chǎng)模擬，這些求解例程可能缺乏效率和魯棒性。 Hill (1994) 提供了預(yù)條件共軛梯度 (PCG2) 包，為 MODFLOW 提供了替代求解器。

PCG2 中的預(yù)處理方案是 Saad (1985) 提出的最小二乘多項(xiàng)式預(yù)處理器或 Meyer 等人 (1989) 的最優(yōu)切比雪夫多項(xiàng)式預(yù)處理器。選擇這些方案主要是出于計(jì)算機(jī)存儲(chǔ)方面的考慮，并且在大規(guī)模現(xiàn)場(chǎng)研究中通常表現(xiàn)不佳。認(rèn)識(shí)到這一點(diǎn)，我們合并了 PCG4 包，其中包含基于部分 LU 分解的 PCG 求解器作為預(yù)處理器。計(jì)算機(jī)存儲(chǔ)要求比以前使用 MODFLOW 的求解器要求的要多。然而，PCG4 求解器簡(jiǎn)單、穩(wěn)健且高效。該求解器已用于許多不同規(guī)模和復(fù)雜程度的地下水建模項(xiàng)目。本章的其余部分致力于 PCG4 包的制定和輸入其使用說明。

6.2 標(biāo)準(zhǔn)PCG方案的制定

6.3 PCG4 封裝的輸入指令

6.3.1 PCG4 輸入

6.3.2 輸入指令中使用的字段說明

第 7 章：帶回溯 NRB1 包的牛頓-拉夫森線性化

7.1 概述 .?

以水或空氣為活性相的嚴(yán)格非飽和區(qū)流動(dòng)模擬通常是高度非線性的，需要小時(shí)間步長(zhǎng)和多次迭代才能解決實(shí)際情況。對(duì)于某些情況，一步穩(wěn)態(tài)模擬可能無法與無約束情況下的 MODFLOW 皮卡德迭代方案收斂，或者對(duì)于非飽和流（空氣或水）情況下的 MODFLOW-SURFACT/MODHMS（第 2 章）的修改皮卡德方案可能無法收斂，或用偽土壤函數(shù)確定無承壓含水層中的地下水位高程，其中不需要嚴(yán)格的非飽和流模擬。提供二次收斂的Newton-Raphson方案可以極大地緩解收斂困難。這里討論的 NRB1 包將 Newton-Raphson 方案與回溯算法集成到 MODFLOW 中以控制步長(zhǎng)，以用于非線性情況，尤其是當(dāng) Picard 迭代過多或完全失敗時(shí)，這會(huì)受益。

NRB1 包包括 Newton-Raphson 線性化和回溯方案，以穩(wěn)定牛頓迭代。該公式被設(shè)計(jì)為與當(dāng)前在每次迭代中求解所有節(jié)點(diǎn)處的頭的 MODFLOW 方案兼容。

此外，該公式不需要對(duì) MODFLOW 當(dāng)前模擬的任何邊界條件進(jìn)行特殊處理。最后，回溯方案限制了任何迭代中殘差的增加，而欠松弛技術(shù)有助于迭代之間解決方案的振蕩行為。這些問題將在本章的其余部分討論

7.2 牛頓拉夫森方案的制定

7.3 NRB1 封裝的輸入指令

第 8 章：觀察節(jié)點(diǎn) (OBS1) 包

8.1 概述

MODFLOWSURFACT/MODHMS 中包含觀察節(jié)點(diǎn)包，以簡(jiǎn)化域內(nèi)任何節(jié)點(diǎn)處瞬態(tài)模擬突破曲線的后處理。流量模擬輸出顯示模擬所請(qǐng)求節(jié)點(diǎn)處的壓力-時(shí)間關(guān)系，傳輸輸出包括模擬的每種污染物的濃度-時(shí)間關(guān)系。創(chuàng)建一個(gè)包含突破輸出的單獨(dú)輸出文件，可以輕松對(duì)其進(jìn)行后處理以提供突破曲線

8.2 OBS1 包輸入指令

8.2.1 OBS1 輸入

8.2.2 輸入指令中使用的字段說明

8.3 OBS1 包的輸出

第 9 章：參考資料

標(biāo)簽：