專利名稱：一種地下水環(huán)境的建模及數(shù)值模擬方法
技術領域：
本發(fā)明涉及一種數(shù)值模擬方法，特別是關于一種地下水環(huán)境的建模及數(shù)值模擬方法。
背景技術：
地下水環(huán)境突發(fā)污染事故需要快速掌握污染物在地下水環(huán)境中的遷移轉化規(guī)律， 確定污染物擴散范圍，進而迅速實施有針對性的治理與防范。由此要求地下水環(huán)境數(shù)值模 型具有快速高效的建模與計算能力。目前，地下水環(huán)境的數(shù)值模擬方法主要有歐拉法、拉格朗日法、混合歐拉_拉格朗 日法等。基于這些方法構建了 MT3D、HST3D、WORM、SWIFT、USGS2D-M0C等地下水環(huán)境數(shù)值模 擬軟件。這些軟件各有其優(yōu)缺點，但它們或因局限于飽和帶或非飽和帶單一含水層污染物 遷移模擬；或因構建模型程序復雜，難以適應快速建模的需要；或因計算單元較多，計算效 率低，難以應用于地下水環(huán)境突發(fā)污染事故應急數(shù)值模擬等，因此，不能對地下水環(huán)境突發(fā) 污染事故進行全面、高效、準確的模擬計算。
發(fā)明內(nèi)容
針對上述問題，本發(fā)明的目的是提供一種地下水環(huán)境的建模及數(shù)值模擬方法，該 方法可對地下水環(huán)境的突發(fā)污染事故進行全面、高效、準確的模擬計算。為實現(xiàn)上述目的，本發(fā)明采取以下技術方案一種地下水環(huán)境的建模及數(shù)值模擬 方法，其包括以下步驟1)前處理設置兩個屬性文件，一是模擬區(qū)土壤屬性文件，包括研 究區(qū)邊界、X與Y方向網(wǎng)格間距、區(qū)域土壤滲透系數(shù)、土壤有效孔隙度、污染物在此種土壤中 的縱向和橫向的彌散度、土壤非飽和帶擴散系數(shù)和溶質(zhì)吸附系數(shù)；二是地表高程與地下水 水位屬性文件；應用IDL軟件的空間插值函數(shù)KRIG2D0或TRI_SURF()對地表高程、地下水 位進行空間網(wǎng)格插值，同時劃分空間有限差分數(shù)值模擬網(wǎng)格，得到兩個相同的差分網(wǎng)格矩 陣地表高程矩陣matrixl和地下水位矩陣matrix2 ；應用SIZEO函數(shù)確定研究區(qū)網(wǎng)格屬 性，包括單元網(wǎng)格數(shù)、單元網(wǎng)格行列數(shù)；應用隨機數(shù)值生成函數(shù)FINDGEN0生成隨機數(shù)，然 后將網(wǎng)格計算單元的節(jié)點坐標生成序列數(shù)組；2) 二維水流模型計算二維水流模型計算包 括溶質(zhì)非飽和帶遷移模擬的地下水埋深矩陣計算以及飽和帶二維地下水流場計算部分； 利用前處理得到的地表高程矩陣與地下水位矩陣相減得非飽和帶厚度Iv即地下水位埋深 矩陣 matrixhO h0:matrixhO = matrixl_matrix2該矩陣直接用于非飽和帶溶質(zhì)遷移模擬計算；飽和帶二維溶質(zhì)遷移模型的水動力 條件是地下水流場的分布，該流場由二維水流模型計算得到，應用二維水流模型計算地下 水流場①計算單元節(jié)點流速利用達西定律
KIV =—
其中，ν為流速；K為滲透系數(shù)；I為水力坡度；ne為有效孔隙度；②計算單元邊界 面流速求得節(jié)點流速，根據(jù)網(wǎng)格中相鄰兩個節(jié)點的流速，取兩個節(jié)點流速的平均值，計算 得到網(wǎng)格單元邊界面上流速分量的IDL矩陣；3)建立非飽和帶一維溶質(zhì)運移模型；
AHμ、h =- ⑴
vZ 其中，ti為溶質(zhì)在非飽和帶中的運移總時間；ΔΗ為非飽和帶厚度，AH = h0;Vz = D/Δ ζ表示假設非飽和帶厚度為Δ ζ = Im時，污染物垂向擴散速度，D為垂向擴散系數(shù)，量 綱為L2/T ；4)建立飽和帶溶質(zhì)遷移二維數(shù)值模型
DdC d ^ 5C、 θ .n 3C、 θ /ΤΛ dC、 δ ΓΛ 5C,丟(VxC)-要(VyC)+ ^cs-A1C-X1^-C ox oy θθ其中，R為等溫吸附系數(shù)，C為溶解濃度，單位為ML—3，、為飽和帶中溶質(zhì)遷移數(shù)值 模擬的時間，Dij為彌散系數(shù)張量，單位為L2T-1, χ、y為計算距離L，vx和Vy為地下水流速， 單位為Lr1Mqs為源/匯處單位體積含水層流量(，單位為Μ3!1—1，θ為含水層孔隙度，Cs為 源/匯項濃度，單位為ML_3，X1為溶解相的反應速率常數(shù)，單位為Γ1，λ 2為吸附相的反應 速率常數(shù)，單位為Γ1，Pb為空隙介質(zhì)體積密度，單位為ML_3，f為吸附濃度，單位為MT1 ；5) 將公式(2)中的各項式進行近似表示在規(guī)則間距的節(jié)點中心網(wǎng)格中，公式(2)右邊第一項 在單元(i，j)近似為3 (Dxx dC)1 ’川/2)— C'->)“、cu 一cU-I)(3) QxxxQx Ax Ax公式(3)表示χ方向上，由濃度梯度引起從χ方向上進入單元(i，j)的凈彌散通 量；公式(2)右邊第二項近似為
(《)
1 Ay(iJ+V2) (^i+lj+l/2 — ^j-lj+l/l) — DxyiJJ-鴉(CMJ_y2 ~ C,_lJ_y2 ) Δ^IAy(4) = 4AxAy [ A沖 J+l/2) (CMJ+l + CMJ _ CMJ+1 — Ci-XJ )-Dxyii J_V2)(CMj + CM j_x -C…-C1^1)]公式(4)表示由y方向的濃度梯度引起的從χ方向進入單元(i，j)的凈彌散通 量；同理，得到由y方向的濃度梯度引起的從y方向上進入單元(i，j)的凈彌散通量，即公 式⑵右邊的第三項為
「00221 —(D 工^ 丄 '+徹⑷+1" — C'，》-Z^d)(C'" —U(5)由χ方向的濃度梯度引起的從y方向進入單元(i，j)的凈彌散通量，即公式(2) 右邊的第四項，寫為
a dc、— (^―) oy οχ
6
__1 ^yx(i+l/2J) (Q+l/2,y+l _ QH-1/2J-1) ~ ^( -1/2,；) (^-1/2,7+1 ~ ^/-1/2,7-1)
~ 4y2Ax(6)= J^c [^(-1/2,7) (q+U+i + C>,M ~ — CU-i)—Dyx(l_l/2 J) (Cl J+l + CM J+1 - Ci H - Ct_l H)]公式⑵右邊第五項為從χ方向進入單元(i，j)的凈對流通量，近似為 ^ 1— (vxC) —{vx(iJ+V2) [(1 -a)Cu + aCtJ+x\-( )v^^l-^C^+aC^]}第六項為從y方向進入單元(i，j)的凈對流通量，近似為
51— (VyC) — {vy(M,2J) [(ι _ cc)Cu + aCMJ ]-(8)v>(,_1/2j)[(l- )CMJ+aC;J]}其中，α為空間加權因子；中心加權，α等于0. 5 ；采用上風加權方案
ο 若 v>0
a — {公式(2)右邊第七項是由源/匯引起的流入或流出單元(i，j)凈質(zhì)量通量；公式 (2)右邊第八項和第九項是單元(i，j)內(nèi)一級化學反應引起的質(zhì)量損失或增加量；吸附濃 度6表示為溶解濃度的函數(shù)，具體取決于所用的吸附等溫線；公式⑵的左邊是儲存于單元(i，j)的質(zhì)量變化率，近似為
A 廣 Cn+l — C"R^L^R^d——(9)
dtAt其中延遲因子R由吸附等溫線確定；6)對以上各分量差分計算過程詳細實現(xiàn)過程 如下①彌散系數(shù)D 由邊界流速分量及公式(10)，利用IDL矩陣運算計算彌散系數(shù)分量；
22
ρ, — aLVx(iJ±l/2)aTVy(,J±l/2) n*
2j^(Uiyi) — —+ ~~ 十 u
"J±l/2V' J±l/2
Π_ aLVy(i±\l2j) , aTVxU±l/2,j) ,
uMM^J) ~ ~“十十 U
V'±l/2,yVi±i/2J
Γ,-(rv ry \ V<U±V2)Vy(iJ±V2)
(10)
Dyx(,±\!2j) =CaL ~ατ)
VU±V2 Vx(i±/2J)Vy(i±l/2J)
V'±V2j 其中為有效分子擴散系數(shù)v,=.fv2 +V2(11)
·* V (·,*) .κ·.·)②彌散項設置初始溶質(zhì)濃度場，即確定污染事件發(fā)生位置及該物質(zhì)的飽和溶解
7度作為模型的初始濃度分布矩陣C ；由彌散項計算公式(3)、公式(4)、公式(5)及公式(6)， 結合彌散系數(shù)計算結果矩陣，利用IDL矩陣運算計算各彌散項；③對流項由對流項計算 公式(7)、公式(8)，結合單元格邊界流速分量、土壤彌散度屬性，利用IDL矩陣運算計算各 對流項；④公式(2)的第七項，源匯項由于快速處理模型、源匯項設置為0 ；⑤公式(2)的 第八、第九項，化學反應項根據(jù)各化學物質(zhì)的等溫吸附曲線獲取；⑥單元節(jié)點濃度由公 式(9)及各個分量矩陣計算結果，利用IDL矩陣運算計算單元節(jié)點后一時刻濃度場；7)返 回步驟5)中的計算過程，即實現(xiàn)一個時間段后，區(qū)域流場及溶質(zhì)濃度場的計算；在步驟5) 前加一個總時段的循環(huán)控制語句，即實現(xiàn)各時段區(qū)域流場及溶質(zhì)濃度場的計算；8)結果輸 出經(jīng)過上述步驟實時輸出地表高程矩陣matrixl、地下水位矩陣matrix2、地下水埋深矩 陣matrixhO、地下水流場、溶質(zhì)濃度場和時間項。所述步驟4)中，根據(jù)中心加權方案，從節(jié)點濃度獲得界面濃度 +υ+1/「 +1"(12)本發(fā)明由于采取以上技術方案，其具有以下優(yōu)點1、本發(fā)明建立的一維模型和二 維模型的結構簡單，模型核心計算全部采用矩陣計算，矩陣公式與算法公式完全一致，易于 理解與建模。2、本發(fā)明由于充分運用了 IDL平臺提供的函數(shù)庫，對于空間插值、矩陣計算等 原本在Fortran等平臺下十分繁瑣的工作，本模型多數(shù)僅需一行代碼即可完成，全部模型 僅有約180行代碼，模型代碼精煉。3、本發(fā)明由于核心處理采用了基于IDL矩陣及函數(shù)計 算，模型計算效率大為提高。4、本發(fā)明針對地下水環(huán)境突發(fā)污染事件，選定模擬區(qū)域后，可 快速插值生成指定大小的計算網(wǎng)格；并根據(jù)污染區(qū)域周邊少數(shù)觀測井的地表高程、潛水水 位、地表水水位，快速插值生成近似于實際條件的地形、水位及水位埋深等；通過快速建立 的模型可指定污染物位置、濃度等。5、本發(fā)明采用有限差分方法計算單元溶質(zhì)的遷移擴散， 完全擺脫了以往的Fortran語言或C語言的數(shù)組循環(huán)計算，取而代之的是IDL平臺下的矩 陣運算方法，計算簡單高效。本發(fā)明基于IDL平臺，可高效、快速、簡便地模擬地下水環(huán)境， 因此，可廣泛用于地下水環(huán)境突發(fā)污染事件的數(shù)值模擬、預測、應急預測等過程中。
圖1是本發(fā)明流程示意2是本發(fā)明計算由χ以及y方向的濃度梯度引起的χ方向上的彌散通量示意3是本發(fā)明快速處理模塊地表高程插值結果示意4是本發(fā)明快速處理模塊地下水位插值結果示意5是本發(fā)明快速處理模塊潛水埋深計算結果示意6是本發(fā)明實施例中地下水環(huán)境應急快速預警模型模擬結果示意7是本發(fā)明實施例中當t = 1000天時區(qū)域流場與污染場疊加分布示意圖
具體實施例方式下面結合附圖和實施例對本發(fā)明進行詳細的描述。本發(fā)明是一種耦合土壤非飽和帶、飽和帶溶質(zhì)遷移過程的地下水環(huán)境數(shù)值模型的 快速構建與準確模擬方法，其基于以下思想由于污染事件通常發(fā)生于地表，污染物首先經(jīng)
8過含水層上部非飽和帶垂向運動，然后進入含水層飽和帶；本發(fā)明對在非飽和帶含水層中 的溶質(zhì)遷移采用垂向一維模型，而在飽和含水層中，由于溶質(zhì)垂向遷移擴散速率往往較橫 向遷移擴散速率小一個數(shù)量級，因此，基于快速模擬計算忽略飽和含水層的垂向運動，采用 了水平二維數(shù)值模擬方法。本發(fā)明包括以下步驟，如圖1所示1)前處理設置兩個屬性文件，一是模擬區(qū)土壤屬性文件(Soilpro文件)，包括研究區(qū)邊界 (模擬區(qū)范圍)、X與Y方向網(wǎng)格間距、區(qū)域土壤滲透系數(shù)、土壤有效孔隙度、污染物在此種 土壤中的縱向和橫向的彌散度、土壤非飽和帶擴散系數(shù)和溶質(zhì)吸附系數(shù)，如表1所示；二是 地表高程與地下水水位屬性文件(Wellpro文件)，如表2所示。應用IDL軟件的空間插值函數(shù)KRIG2D0或TRI_SURF()對地表高程、地下水位進 行空間網(wǎng)格插值，同時劃分空間有限差分數(shù)值模擬網(wǎng)格，得到兩個相同的差分網(wǎng)格矩陣地 表高程矩陣(matrixl)和地下水位矩陣(matrix2)。應用SIZEO函數(shù)確定研究區(qū)網(wǎng)格屬 性，包括單元網(wǎng)格數(shù)、單元網(wǎng)格行列數(shù)等。應用隨機數(shù)值生成函數(shù)FINDGEN()生成隨機數(shù)， 然后將網(wǎng)格計算單元的節(jié)點坐標生成序列數(shù)組。表ISoilpro的文件格式 表2Wellpro的文件格式 2) 二維水流模型計算二維水流模型計算包括溶質(zhì)非飽和帶遷移模擬的地下水埋深矩陣計算以及飽和 帶二維地下水流場計算部分。利用前處理得到的地表高程矩陣與地下水位矩陣相減可得非飽和帶厚度Iv即地 下水位埋深矩陣matrixhO h0:matrixhO = matrixl_matrix2該矩陣可直接用于非飽和帶溶質(zhì)遷移模擬計算。飽和帶二維溶質(zhì)遷移模型的水動力條件是地下水流場的分布，該流場可由二維水 流模型計算得到，應用二維水流模型計算地下水流場①計算單元節(jié)點流速利用達西定律
9KIV =—其中，ν為流速；K為滲透系數(shù)；I為水力坡度；ne為有效孔隙度。根據(jù)達西定律和表1，利用IDL矩陣計算節(jié)點流速代碼如下 ②計算單元邊界面流速求得節(jié)點流速，根據(jù)網(wǎng)格中相鄰兩個節(jié)點的流速，取兩個節(jié)點流速的平均值，計算 得到網(wǎng)格單元邊界面上流速分量的IDL矩陣，計算代碼如下
VelYX=DBLARR(nDG
, nDG[l]-l) ；Y面上(縱向面)X方向網(wǎng)格面中心點實際流速
；滲透系數(shù)公式中的Vx(i, j+-l/2) VelYX
= (VelX[0:nDG
-1, 0:nDG[l]-2]+VelX[0:nDG
-l, l:nDG[l]_l])*0. 5 3)非飽和帶一維溶質(zhì)運移模型計算非飽和帶計算采用簡化一維溶質(zhì)運移模型，僅計算溶質(zhì)由地表運移至含水層的時
間，計算公式如下
AH…tx =--U)
vZ其中，、為溶質(zhì)在非飽和帶中的運移總時間；ΔΗ = Iltl為非飽和帶厚度；vz = D/ Δ ζ表示假設非飽和帶厚度為Δ ζ = Im時，污染物垂向擴散速度，D為垂向擴散系數(shù)，量綱 為(L2/T)。4)飽和帶溶質(zhì)遷移二維數(shù)值模型計算采用有限差分法的歐拉方法，計算地下水溶質(zhì)遷移擴散過程。模型飽和含水層溶 質(zhì)對流_彌散遷移采用二維對流_彌散方程如下
Dec d ,n ac, d ,n sc. d ,n dc、 d ,n dc. _8] ⑵I (vxC)-備(vyC)
dxoyθθ其中，R為等溫吸附系數(shù)，C為溶解濃度(ML—3)，、為飽和帶中溶質(zhì)遷移數(shù)值模擬 的時間，Dij為彌散系數(shù)張量(L2T_0，χ、y為計算距離(L)，vx和Vy為地下水流速(L Γ1)，，
10qs為源/匯處單位體積含水層流量(Μ3!1—1)，θ為含水層孔隙度，Cs為源/匯項濃度(ML—3)， λ工為溶解相的反應速率常數(shù)(Γ1)，λ 2為吸附相的反應速率常數(shù)(Γ1)，P b為空隙介質(zhì)體 積密度(ML—3)，5為吸附濃度(MM—1)。該模型充分考慮了彌散、遷移、源匯項、化學吸附等影 響。
5)將公式(2)中的各項式進行近似表示，以便采用有限差分法編程實現(xiàn)。 在規(guī)則間距的節(jié)點中心網(wǎng)格中，公式⑵右邊第一項在單元(i，j)近似為
0091]
0092]
0093]—(D SC) .-■ 1 d^-J+V2')(C',·/+1 ~C'J)“J-i/2)(C。. ~Cu-i)(3) QxxxQx Ax Ax
0094]如圖2(a)所示，公式(3)表示x方向上，由濃度梯度引起從x方向上進入單元(i， j)的凈彌散通量。
0095]公式(2)右邊第二項可近似為
d …dc、
0096]—(Dxy—) ox oy
Q097-,1 Dxy(iJ+l/2) (CMJ+l/2 一 Ci-\J+/2) — Dxy(iJ-l/2) (CMJ—l/2 ~ Ci-\J-l/2 )
IAy(4)
0098]=’川 + Cmj - C,_u+1 - Cm 7)-
0099]Dxy(,j-y2) (CM,J + CMJ-\ — Ci-l,j 一 Ci—l’j-l)]
0100]如圖2(b)所示，公式⑷表示由y方向的濃度梯度引起的從χ方向進入單元(i， j)的凈彌散通量。其中用到中心加權方案，從節(jié)點濃度獲得界面濃度，圖中以 標記出
0101]C,.+u+1/2 = C'+U+12+C'+U(5)
0102]同理，可以得到由y方向的濃度梯度引起的從y方向上進入單元(i，j)的凈彌散 通量，即公式(2)右邊的第三項為
Q103JdC) -■ 1 DMM/2J) (CMJ — CU) ~ Dyy{i-l/2j) (CiJ — C,-IJ)
dy yySy AyAy
0104] 由χ方向的濃度梯度引起的從y方向進入單元(i，j)的凈彌散通量，即公式(2) 右邊的第四項，可寫為
dy y dx
、1 ^yx(i+l/2J) (Q+1/2J+1 一 ^/+IAv-I) ‘ ^γχ( -1/2,β (^/-1/2,7+1 一 Q-l/2,y-l )
0105]
0106]
Ay2Δχ(7)
0107] = ^^ [^(1+i/2j) (Ci+1 J+1 + Ci j+l - CMJ_, - Ct ]_x)-
01 08] DyX[,-V2,j) (C, J+1 + C,-xJ+\ - Cij-x - Cm^1 )]
0109] 公式(3)、公式(4)、公式(6)及公式(7)中的彌散系數(shù)分量的計算，若給定縱向彌 散度與橫向彌散度α τ，各彌散系數(shù)計算公式如下凡(,>;±1/2)=^ 1 + ^^ + 廣
V'J±V2ViJ±V2
11CN 101908100 A
說明 書8/10頁
2 2 Γ π π- aLVy(i±V2j) , aTVx(i±l/2J) ,Dyym/2j)-—-+ —-+ L)
Vi±V2JVi±l/2JDxy(i,j±\ll) = (aL ~ατ) X('J±i/2) y{'J±V2)
VU±V2
Γ Π-fry ry Λ V ^>±V2J)V yQ±i/2J)νγχ(,±1/2,β ~ KaL ~ aT )-
V'±l/2,7其中為有效分子擴散系數(shù)。vtt=.fv2 W~
/"N 8 /Is
(9)公式(2)右邊第五項為從χ方向進入單元(i，j)的凈對流通量，近似為
Q1
&[(1 +( ο)VX(IJ_1/2) [(1 -+ aCu ]}第六項為從y方向進入單元(i，j)的凈對流通量，可近似為
Q1— (VyC) — {vr(i+1/2J) [(1 — a)CtJ + aCMJ ] 一(n)
其中α為空間加權因子。中心加權，α等于0.5;采用上風加權方案
Γθ 若 V>0 a = <公式⑵右邊第七項是由源/匯引起的流入或流出單元(i，j)凈質(zhì)量通量，第八 項和第九項是單元(i，j)內(nèi)一級化學反應引起的質(zhì)量損失或增加量。吸附濃度G可以表示 為溶解濃度的函數(shù)，具體取決于所用的吸附等溫線。
公式(2)的左邊是儲存于單元(i，j)的質(zhì)量變化率，可以近似為
Λ廣廣"+1 _ f^n
R-^Rlj(12)
dtAt
其中延遲因子R由吸附等溫線確定。
6)對以上各分量差分計算過程詳細實現(xiàn)過程如下
①彌散系數(shù)D 由邊界流速分量及公式(8)，利用IDL矩陣運算可以計算彌散系數(shù)
分量；②彌散項設置初始溶質(zhì)濃度場，即確定污染事件發(fā)生位置及該物質(zhì)的飽和溶解 度作為模型的初始濃度分布矩陣C。由彌散項計算公式(3)、公式(4)、公式(6)及公式(7)， 結合彌散系數(shù)計算結果矩陣，利用IDL矩陣運算可以計算各彌散項；③對流項由對流項計算公式(10)、公式(11)，結合單元格邊界流速分量、土壤彌 散度屬性，利用IDL矩陣運算可以計算各對流項；④公式(2)的第七項，源匯項由于快速處理模型，源匯項可以暫時設置為O ；⑤公式(2)的第八、第九項，化學反應項根據(jù)各化學物質(zhì)的等溫吸附曲線獲?。?br>
12
⑥單元節(jié)點濃度由公式(12)及各個分量矩陣計算結果，利用IDL矩陣運算可以 計算單元節(jié)點后一時刻濃度場；計算單元節(jié)點后一時刻濃度場的IDL矩陣計算代碼如下cinit[l:nDG
-l,l:nDG[l]-l] = cinit [1 :nDG
_1，1 :nDG[l] _1]+$(deltT/RR)*(Dcxx
+Dcxy
+Dcyy
+$Dcyx
+Vqcx
+Vqcy
+Lambda)其中，cinitG為濃度場矩陣，deltT為時間步長，RR為延遲因子，DcijG為彌散 項，Vqcx[]為對流項，Lambda為等溫吸附項。7)然后返回步驟5)中的計算過程，即可實現(xiàn)一個時間段后，區(qū)域流場及溶質(zhì)濃度 場的計算；在步驟5)前加一個總時段的循環(huán)控制語句，即可實現(xiàn)各時段區(qū)域流場及溶質(zhì)濃 度場的計算。8)結果輸出經(jīng)過上述步驟可以實時輸出地表高程矩陣matrixl、地下水位矩陣 matrix2、地下水埋深矩陣matrixhO、地下水流場、溶質(zhì)濃度場和時間項(時間段)。上述步驟8)輸出的結果以及二維水流模型結果可與地表水環(huán)境模型接口連接， 用于地表水環(huán)境，數(shù)值模擬河道側向入流地下水的溶質(zhì)邊界條件，且輸出結果可根據(jù)地表 水環(huán)境模型輸入文件的格式需要設計輸出相應的文件格式。下面列舉一具體實施例1、水流模擬假設任意一個地方發(fā)生污染事件，利用前處理模塊可以快速得到模型研究區(qū)域邊 界，設定模型剖分網(wǎng)格精度、土壤屬性等，如表3所示，為土壤屬性即Soilpro的文件。根據(jù) 研究區(qū)提供的5眼地下水位觀測井的資料，如表4所示，為地下水位和地表高程即Wellpr0 的文件，可以插值得到該區(qū)域地表高程，如圖3所示。地下水水位插值結果，如圖4所示；地 下水埋深插值結果，如圖5所示。圖3 5中的(a)圖為采用Kringing方法的差值結果圖， (b)圖所示為采用五次樣條光滑差值方法的差值結果圖。表3Soilpro 的文件 2、溶質(zhì)遷移模擬假設污染發(fā)生位置坐標為(20238486,3599112)，假設該污染物在含水層中飽和溶 解度為10000mg/L。通過模擬計算，可得到以下結果(1)污染物由地表運移通過非飽和帶時間為32. 11天。該處地下水埋深為6. 422m。(2)污染物進入含水層5天后向東北方向最遠擴散距離為140m(擴散標準為區(qū)域 污染物濃度達到0. 01mg/L)；污染物進入含水層10天后向東北方向最遠擴散距離為141m ； 污染物進入含水層30天后向東北方向最遠擴散距離為214m ；污染物進入含水層100天后 向東北方向最遠擴散距離為278m ；污染物進入含水層500天后向東北方向最遠擴散距離為 390m ；污染物進入含水層1000天后向東北方向最遠擴散距離為450m。各時間節(jié)點污染物 濃度分布圖，如圖6所示。(3)通過地下水流場與第1000天污染物濃度場疊加圖，如圖7所示，可以看出，模 型模擬污染物遷移擴散趨勢完全符合流場變化規(guī)律。上述實例的關鍵技術在于實例研究中，研究區(qū)共計剖分51506個計算單元，計算時間步長1天，時間總長度 為1000天，模型輸出為1000天的各節(jié)點地表高程、水位、水位埋深、溶質(zhì)濃度。計算機硬件 配置為 Intel (R) Core (TM) 2 Quad Q6600 2. 40GHz CPU，內(nèi)存 4GB (PC2-5300 DDR2 666MHz)， 硬盤為希捷 ST3500820AS (498GB, 7200 轉 / 分)。操作系統(tǒng)為 Windows XP (32bit/SP2)。模 型計算耗時50. 18秒。上述各實施例僅用于說明本發(fā)明，其中各部件的結構、連接方式等都是可以有所 變化的，凡是在本發(fā)明技術方案的基礎上進行的等同變換和改進，均不應排除在本發(fā)明的 保護范圍之外。
1權利要求
一種地下水環(huán)境的建模及數(shù)值模擬方法，其包括以下步驟1)前處理設置兩個屬性文件，一是模擬區(qū)土壤屬性文件，包括研究區(qū)邊界、X與Y方向網(wǎng)格間距、區(qū)域土壤滲透系數(shù)、土壤有效孔隙度、污染物在此種土壤中的縱向和橫向的彌散度、土壤非飽和帶擴散系數(shù)和溶質(zhì)吸附系數(shù)；二是地表高程與地下水水位屬性文件；應用IDL軟件的空間插值函數(shù)KRIG2D()或TRI_SURF()對地表高程、地下水位進行空間網(wǎng)格插值，同時劃分空間有限差分數(shù)值模擬網(wǎng)格，得到兩個相同的差分網(wǎng)格矩陣地表高程矩陣matrix1和地下水位矩陣matrix2；應用SIZE()函數(shù)確定研究區(qū)網(wǎng)格屬性，包括單元網(wǎng)格數(shù)、單元網(wǎng)格行列數(shù)；應用隨機數(shù)值生成函數(shù)FINDGEN()生成隨機數(shù)，然后將網(wǎng)格計算單元的節(jié)點坐標生成序列數(shù)組；2)二維水流模型計算二維水流模型計算包括溶質(zhì)非飽和帶遷移模擬的地下水埋深矩陣計算以及飽和帶二維地下水流場計算部分；利用前處理得到的地表高程矩陣與地下水位矩陣相減得非飽和帶厚度h0，即地下水位埋深矩陣matrixh0h0matrixh0＝matrix1 matrix2該矩陣直接用于非飽和帶溶質(zhì)遷移模擬計算；飽和帶二維溶質(zhì)遷移模型的水動力條件是地下水流場的分布，該流場由二維水流模型計算得到，應用二維水流模型計算地下水流場①計算單元節(jié)點流速利用達西定律 <mrow><mi>v</mi><mo>=</mo><mfrac> <mi>KI</mi> <msub><mi>n</mi><mi>e</mi> </msub></mfrac> </mrow>其中，v為流速；K為滲透系數(shù)；I為水力坡度；ne為有效孔隙度；②計算單元邊界面流速求得節(jié)點流速，根據(jù)網(wǎng)格中相鄰兩個節(jié)點的流速，取兩個節(jié)點流速的平均值，計算得到網(wǎng)格單元邊界面上流速分量的IDL矩陣；3)建立非飽和帶一維溶質(zhì)運移模型； <mrow><msub> <mi>t</mi> <mn>1</mn></msub><mo>=</mo><mfrac> <mi>ΔH</mi> <msub><mi>v</mi><mi>z</mi> </msub></mfrac><mo>-</mo><mo>-</mo><mo>-</mo><mrow> <mo>(</mo> <mn>1</mn> <mo>)</mo></mrow> </mrow>其中，t1為溶質(zhì)在非飽和帶中的運移總時間；ΔH為非飽和帶厚度，ΔH＝h0；vz＝D/Δz表示假設非飽和帶厚度為Δz＝1m時，污染物垂向擴散速度，D為垂向擴散系數(shù)，量綱為L2/T；4)建立飽和帶溶質(zhì)遷移二維數(shù)值模型 <mrow><mi>R</mi><mfrac> <mrow><mo>∂</mo><mi>C</mi> </mrow> <mrow><mo>∂</mo><msub> <mi>t</mi> <mn>1</mn></msub> </mrow></mfrac><mo>=</mo><mfrac> <mo>∂</mo> <mrow><mo>∂</mo><mi>x</mi> </mrow></mfrac><mrow> <mo>(</mo> <msub><mi>D</mi><mi>xx</mi> </msub> <mfrac><mrow> <mo>∂</mo> <mi>C</mi></mrow><mrow> <mo>∂</mo> <mi>x</mi></mrow> </mfrac> <mo>)</mo></mrow><mo>+</mo><mfrac> <mo>∂</mo> <mrow><mo>∂</mo><mi>x</mi> </mrow></mfrac><mrow> <mo>(</mo> <msub><mi>D</mi><mi>xy</mi> </msub> <mfrac><mrow> <mo>∂</mo> <mi>C</mi></mrow><mrow> <mo>∂</mo> <mi>y</mi></mrow> </mfrac> <mo>)</mo></mrow><mo>+</mo><mfrac> <mo>∂</mo> <mrow><mo>∂</mo><mi>y</mi> </mrow></mfrac><mrow> <mo>(</mo> <msub><mi>D</mi><mi>yy</mi> </msub> <mfrac><mrow> <mo>∂</mo> <mi>C</mi></mrow><mrow> <mo>∂</mo> <mi>y</mi></mrow> </mfrac> <mo>)</mo></mrow><mo>+</mo><mfrac> <mo>∂</mo> <mrow><mo>∂</mo><mi>y</mi> </mrow></mfrac><mrow> <mo>(</mo> <msub><mi>D</mi><mi>yx</mi> </msub> <mfrac><mrow> <mo>∂</mo> <mi>C</mi></mrow><mrow> <mo>∂</mo> <mi>x</mi></mrow> </mfrac> <mo>)</mo></mrow><mo>-</mo><mo>-</mo><mo>-</mo><mo>-</mo><mrow> <mo>(</mo> <mn>2</mn> <mo>)</mo></mrow> </mrow> <mrow><mfrac> <mo>∂</mo> <mrow><mo>∂</mo><mi>x</mi> </mrow></mfrac><mrow> <mo>(</mo> <msub><mi>v</mi><mi>x</mi> </msub> <mi>C</mi> <mo>)</mo></mrow><mo>-</mo><mfrac> <mo>∂</mo> <mrow><mo>∂</mo><mi>y</mi> </mrow></mfrac><mrow> <mo>(</mo> <msub><mi>v</mi><mi>y</mi> </msub> <mi>C</mi> <mo>)</mo></mrow><mo>+</mo><mfrac> <msub><mi>q</mi><mi>s</mi> </msub> <mi>θ</mi></mfrac><msub> <mi>C</mi> <mi>s</mi></msub><mo>-</mo><msub> <mi>λ</mi> <mn>1</mn></msub><mi>C</mi><mo>-</mo><msub> <mi>λ</mi> <mn>2</mn></msub><mfrac> <msub><mi>ρ</mi><mi>b</mi> </msub> <mi>θ</mi></mfrac><mover> <mi>C</mi> <mo>‾</mo></mover> </mrow>其中，R為等溫吸附系數(shù)，C為溶解濃度，單位為ML 3，t1為飽和帶中溶質(zhì)遷移數(shù)值模擬的時間，Dij為彌散系數(shù)張量，單位為L2T 1，x、y為計算距離L，vx和vy為地下水流速，單位為LT 1，，qs為源/匯處單位體積含水層流量(，單位為M3T 1，θ為含水層孔隙度，Cs為源/匯項濃度，單位為ML 1，λ1為溶解相的反應速率常數(shù)，單位為T 1，λ2為吸附相的反應速率常數(shù)，單位為T 1，ρb為空隙介質(zhì)體積密度，單位為ML 3，為吸附濃度，單位為MM 1；5)將公式(2)中的各項式進行近似表示在規(guī)則間距的節(jié)點中心網(wǎng)格中，公式(2)右邊第一項在單元(i，j)近似為 <mrow><mfrac> <mo>∂</mo> <mrow><mo>∂</mo><mi>x</mi> </mrow></mfrac><mrow> <mo>(</mo> <msub><mi>D</mi><mi>xx</mi> </msub> <mfrac><mrow> <mo>∂</mo> <mi>C</mi></mrow><mrow> <mo>∂</mo> <mi>x</mi></mrow> </mfrac> <mo>)</mo></mrow><mo>≈</mo><mfrac> <mn>1</mn> <mi>Δx</mi></mfrac><mfrac> <mrow><msub> <mi>D</mi> <mrow><mi>xx</mi><mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>+</mo> <mn>1</mn> <mo>/</mo> <mn>2</mn> <mo>)</mo></mrow> </mrow></msub><mrow> <mo>(</mo> <msub><mi>C</mi><mrow> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>+</mo> <mn>1</mn></mrow> </msub> <mo>-</mo> <msub><mi>C</mi><mrow> <mi>i</mi> <mo>,</mo> <mi>j</mi></mrow> </msub> <mo>)</mo></mrow><mo>-</mo><msub> <mi>D</mi> <mrow><mi>xx</mi><mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>-</mo> <mn>1</mn> <mo>/</mo> <mn>2</mn> <mo>)</mo></mrow> </mrow></msub><mrow> <mo>(</mo> <msub><mi>C</mi><mrow> <mi>i</mi> <mo>,</mo> <mi>j</mi></mrow> </msub> <mo>-</mo> <msub><mi>C</mi><mrow> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>-</mo> <mn>1</mn></mrow> </msub> <mo>)</mo></mrow> </mrow> <mi>Δx</mi></mfrac><mo>-</mo><mo>-</mo><mo>-</mo><mrow> <mo>(</mo> <mn>3</mn> <mo>)</mo></mrow> </mrow>公式(3)表示x方向上，由濃度梯度引起從x方向上進入單元(i，j)的凈彌散通量；公式(2)右邊第二項近似為 <mrow><mfrac> <mo>∂</mo> <mrow><mo>∂</mo><mi>x</mi> </mrow></mfrac><mrow> <mo>(</mo> <msub><mi>D</mi><mi>xy</mi> </msub> <mfrac><mrow> <mo>∂</mo> <mi>C</mi></mrow><mrow> <mo>∂</mo> <mi>y</mi></mrow> </mfrac> <mo>)</mo></mrow> </mrow> <mrow><mo>≈</mo><mfrac> <mn>1</mn> <mi>Δx</mi></mfrac><mfrac> <mrow><msub> <mi>D</mi> <mrow><mi>xy</mi><mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>+</mo> <mn>1</mn> <mo>/</mo> <mn>2</mn> <mo>)</mo></mrow> </mrow></msub><mrow> <mo>(</mo> <msub><mi>C</mi><mrow> <mi>i</mi> <mo>+</mo> <mn>1</mn> <mo>,</mo> <mi>j</mi> <mo>+</mo> <mn>1</mn> <mo>/</mo> <mn>2</mn></mrow> </msub> <mo>-</mo> <msub><mi>C</mi><mrow> <mi>i</mi> <mo>-</mo> <mn>1</mn> <mo>,</mo> <mi>j</mi> <mo>+</mo> <mn>1</mn> <mo>/</mo> <mn>2</mn></mrow> </msub> <mo>)</mo></mrow><mo>-</mo><msub> <mi>D</mi> <mrow><mi>xy</mi><mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>-</mo> <mn>1</mn> <mo>/</mo> <mn>2</mn> <mo>)</mo></mrow> </mrow></msub><mrow> <mo>(</mo> <msub><mi>C</mi><mrow> <mi>i</mi> <mo>+</mo> <mn>1</mn> <mo>,</mo> <mi>j</mi> <mo>-</mo> <mn>1</mn> <mo>/</mo> <mn>2</mn></mrow> </msub> <mo>-</mo> <msub><mi>C</mi><mrow> <mi>i</mi> <mo>-</mo> <mn>1</mn> <mo>,</mo> <mi>j</mi> <mo>-</mo> <mn>1</mn> <mo>/</mo> <mn>2</mn></mrow> </msub> <mo>)</mo></mrow> </mrow> <mrow><mn>2</mn><mi>Δy</mi> </mrow></mfrac><mo>-</mo><mo>-</mo><mo>-</mo><mrow> <mo>(</mo> <mn>4</mn> <mo>)</mo></mrow> </mrow> <mrow><mo>=</mo><mfrac> <mn>1</mn> <mrow><mn>4</mn><mi>ΔxΔy</mi> </mrow></mfrac><mo>[</mo><msub> <mi>D</mi> <mrow><mi>xy</mi><mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>+</mo> <mn>1</mn> <mo>/</mo> <mn>2</mn> <mo>)</mo></mrow> </mrow></msub><mrow> <mo>(</mo> <msub><mi>C</mi><mrow> <mi>i</mi> <mo>+</mo> <mn>1</mn> <mo>,</mo> <mi>j</mi> <mo>+</mo> <mn>1</mn></mrow> </msub> <mo>+</mo> <msub><mi>C</mi><mrow> <mi>i</mi> <mo>+</mo> <mn>1</mn> <mo>,</mo> <mi>j</mi></mrow> </msub> <mo>-</mo> <msub><mi>C</mi><mrow> <mi>i</mi> <mo>-</mo> <mn>1</mn> <mo>,</mo> <mi>j</mi> <mo>+</mo> <mn>1</mn></mrow> </msub> <mo>-</mo> <msub><mi>C</mi><mrow> <mi>i</mi> <mo>-</mo> <mn>1</mn> <mo>,</mo> <mi>j</mi></mrow> </msub> <mo>)</mo></mrow><mo>-</mo> </mrow> <mrow><msub> <mi>D</mi> <mrow><mi>xy</mi><mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>-</mo> <mn>1</mn> <mo>/</mo> <mn>2</mn> <mo>)</mo></mrow> </mrow></msub><mrow> <mo>(</mo> <msub><mi>C</mi><mrow> <mi>i</mi> <mo>+</mo> <mn>1</mn> <mo>,</mo> <mi>j</mi></mrow> </msub> <mo>+</mo> <msub><mi>C</mi><mrow> <mi>i</mi> <mo>+</mo> <mn>1</mn> <mo>,</mo> <mi>j</mi> <mo>-</mo> <mn>1</mn></mrow> </msub> <mo>-</mo> <msub><mi>C</mi><mrow> <mi>i</mi> <mo>-</mo> <mn>1</mn> <mo>,</mo> <mi>j</mi></mrow> </msub> <mo>-</mo> <msub><mi>C</mi><mrow> <mi>i</mi> <mo>-</mo> <mn>1</mn> <mo>,</mo> <mi>j</mi> <mo>-</mo> <mn>1</mn></mrow> </msub> <mo>)</mo></mrow><mo>]</mo> </mrow>公式(4)表示由y方向的濃度梯度引起的從x方向進入單元(i，j)的凈彌散通量；同理，得到由y方向的濃度梯度引起的從y方向上進入單元(i，j)的凈彌散通量，即公式(2)右邊的第三項為 <mrow><mfrac> <mo>∂</mo> <mrow><mo>∂</mo><mi>y</mi> </mrow></mfrac><mrow> <mo>(</mo> <msub><mi>D</mi><mi>yy</mi> </msub> <mfrac><mrow> <mo>∂</mo> <mi>C</mi></mrow><mrow> <mo>∂</mo> <mi>y</mi></mrow> </mfrac> <mo>)</mo></mrow><mo>≈</mo><mfrac> <mn>1</mn> <mi>Δy</mi></mfrac><mfrac> <mrow><msub> <mi>D</mi> <mrow><mi>yy</mi><mrow> <mo>(</mo> <mi>i</mi> <mo>+</mo> <mn>1</mn> <mo>/</mo> <mn>2</mn> <mo>,</mo> <mi>j</mi> <mo>)</mo></mrow> </mrow></msub><mrow> <mo>(</mo> <msub><mi>C</mi><mrow> <mi>i</mi> <mo>+</mo> <mn>1</mn> <mo>,</mo> <mi>j</mi></mrow> </msub> <mo>-</mo> <msub><mi>C</mi><mrow> <mi>i</mi> <mo>,</mo> <mi>j</mi></mrow> </msub> <mo>)</mo></mrow><mo>-</mo><msub> <mi>D</mi> <mrow><mi>yy</mi><mrow> <mo>(</mo> <mi>i</mi> <mo>-</mo> <mn>1</mn> <mo>/</mo> <mn>2</mn> <mo>,</mo> <mi>j</mi> <mo>)</mo></mrow> </mrow></msub><mrow> <mo>(</mo> <msub><mi>C</mi><mrow> <mi>i</mi> <mo>,</mo> <mi>j</mi></mrow> </msub> <mo>-</mo> <msub><mi>C</mi><mrow> <mi>i</mi> <mo>-</mo> <mn>1</mn> <mo>,</mo> <mi>j</mi></mrow> </msub> <mo>)</mo></mrow> </mrow> <mi>Δy</mi></mfrac><mo>-</mo><mo>-</mo><mo>-</mo><mrow> <mo>(</mo> <mn>5</mn> <mo>)</mo></mrow> </mrow>由x方向的濃度梯度引起的從y方向進入單元(i，j)的凈彌散通量，即公式(2)右邊的第四項，寫為 <mrow><mfrac> <mo>∂</mo> <mrow><mo>∂</mo><mi>y</mi> </mrow></mfrac><mrow> <mo>(</mo> <msub><mi>D</mi><mi>yx</mi> </msub> <mfrac><mrow> <mo>∂</mo> <mi>C</mi></mrow><mrow> <mo>∂</mo> <mi>x</mi></mrow> </mfrac> <mo>)</mo></mrow> </mrow> <mrow><mo>≈</mo><mfrac> <mn>1</mn> <mi>Δy</mi></mfrac><mfrac> <mrow><msub> <mi>D</mi> <mrow><mi>yx</mi><mrow> <mo>(</mo> <mi>i</mi> <mo>+</mo> <mn>1</mn> <mo>/</mo> <mn>2</mn> <mo>,</mo> <mi>j</mi> <mo>)</mo></mrow> </mrow></msub><mrow> <mo>(</mo> <msub><mi>C</mi><mrow> <mi>i</mi> <mo>+</mo> <mn>1</mn> <mo>/</mo> <mn>2</mn> <mo>,</mo> <mi>j</mi> <mo>+</mo> <mn>1</mn></mrow> </msub> <mo>-</mo> <msub><mi>C</mi><mrow> <mi>i</mi> <mo>+</mo> <mn>1</mn> <mo>/</mo> <mn>2</mn> <mo>,</mo> <mi>j</mi> <mo>-</mo> <mn>1</mn></mrow> </msub> <mo>)</mo></mrow><mo>-</mo><msub> <mi>D</mi> <mrow><mi>yx</mi><mrow> <mo>(</mo> <mi>i</mi> <mo>-</mo> <mn>1</mn> <mo>/</mo> <mn>2</mn> <mo>,</mo> <mi>j</mi> <mo>)</mo></mrow> </mrow></msub><mrow> <mo>(</mo> <msub><mi>C</mi><mrow> <mi>i</mi> <mo>-</mo> <mn>1</mn> <mo>/</mo> <mn>2</mn> <mo>,</mo> <mi>j</mi> <mo>+</mo> <mn>1</mn></mrow> </msub> <mo>-</mo> <msub><mi>C</mi><mrow> <mi>i</mi> <mo>-</mo> <mn>1</mn> <mo>/</mo> <mn>2</mn> <mo>,</mo> <mi>j</mi> <mo>-</mo> <mn>1</mn></mrow> </msub> <mo>)</mo></mrow> </mrow> <mrow><mn>2</mn><mi>Δx</mi> </mrow></mfrac><mo>-</mo><mo>-</mo><mo>-</mo><mrow> <mo>(</mo> <mn>6</mn> <mo>)</mo></mrow> </mrow> <mrow><mo>=</mo><mfrac> <mn>1</mn> <mrow><mn>4</mn><mi>ΔyΔx</mi> </mrow></mfrac><mo>[</mo><msub> <mi>D</mi> <mrow><mi>yx</mi><mrow> <mo>(</mo> <mi>i</mi> <mo>+</mo> <mn>1</mn> <mo>/</mo> <mn>2</mn> <mo>,</mo> <mi>j</mi> <mo>)</mo></mrow> </mrow></msub><mrow> <mo>(</mo> <msub><mi>C</mi><mrow> <mi>i</mi> <mo>+</mo> <mn>1</mn> <mo>,</mo> <mi>j</mi> <mo>+</mo> <mn>1</mn></mrow> </msub> <mo>+</mo> <msub><mi>C</mi><mrow> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>+</mo> <mn>1</mn></mrow> </msub> <mo>-</mo> <msub><mi>C</mi><mrow> <mi>i</mi> <mo>+</mo> <mn>1</mn> <mo>,</mo> <mi>j</mi> <mo>-</mo> <mn>1</mn></mrow> </msub> <mo>-</mo> <msub><mi>C</mi><mrow> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>-</mo> <mn>1</mn></mrow> </msub> <mo>)</mo></mrow><mo>-</mo> </mrow> <mrow><msub> <mi>D</mi> <mrow><mi>yx</mi><mrow> <mo>(</mo> <mi>i</mi> <mo>-</mo> <mn>1</mn> <mo>/</mo> <mn>2</mn> <mo>,</mo> <mi>j</mi> <mo>)</mo></mrow> </mrow></msub><mrow> <mo>(</mo> <msub><mi>C</mi><mrow> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>+</mo> <mn>1</mn></mrow> </msub> <mo>+</mo> <msub><mi>C</mi><mrow> <mi>i</mi> <mo>-</mo> <mn>1</mn> <mo>,</mo> <mi>j</mi> <mo>+</mo> <mn>1</mn></mrow> </msub> <mo>-</mo> <msub><mi>C</mi><mrow> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>-</mo> <mn>1</mn></mrow> </msub> <mo>-</mo> <msub><mi>C</mi><mrow> <mi>i</mi> <mo>-</mo> <mn>1</mn> <mo>,</mo> <mi>j</mi> <mo>-</mo> <mn>1</mn></mrow> </msub> <mo>)</mo></mrow><mo>]</mo> </mrow>公式(2)右邊第五項為從x方向進入單元(i，j)的凈對流通量，近似為 <mrow><mfrac> <mo>∂</mo> <mrow><mo>∂</mo><mi>x</mi> </mrow></mfrac><mrow> <mo>(</mo> <msub><mi>v</mi><mi>x</mi> </msub> <mi>C</mi> <mo>)</mo></mrow><mo>≈</mo><mfrac> <mn>1</mn> <mi>Δx</mi></mfrac><mo>{</mo><msub> <mi>v</mi> <mrow><mi>x</mi><mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>+</mo> <mn>1</mn> <mo>/</mo> <mn>2</mn> <mo>)</mo></mrow> </mrow></msub><mo>[</mo><mrow> <mo>(</mo> <mn>1</mn> <mo>-</mo> <mi>α</mi> <mo>)</mo></mrow><msub> <mi>C</mi> <mrow><mi>i</mi><mo>,</mo><mi>j</mi> </mrow></msub><mo>+</mo><mi>α</mi><msub> <mi>C</mi> <mrow><mi>i</mi><mo>,</mo><mi>j</mi><mo>+</mo><mn>1</mn> </mrow></msub><mo>]</mo><mo>-</mo><mo>-</mo><mo>-</mo><mo>-</mo><mrow> <mo>(</mo> <mn>7</mn> <mo>)</mo></mrow> </mrow> <mrow><msub> <mi>v</mi> <mrow><mi>x</mi><mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>-</mo> <mn>1</mn> <mo>/</mo> <mn>2</mn> <mo>)</mo></mrow> </mrow></msub><mo>[</mo><mrow> <mo>(</mo> <mn>1</mn> <mo>-</mo> <mi>α</mi> <mo>)</mo></mrow><msub> <mi>C</mi> <mrow><mi>i</mi><mo>,</mo><mi>j</mi><mo>-</mo><mn>1</mn> </mrow></msub><mo>+</mo><mi>α</mi><msub> <mi>C</mi> <mrow><mi>i</mi><mo>,</mo><mi>j</mi> </mrow></msub><mo>]</mo><mo>}</mo> </mrow>第六項為從y方向進入單元(i，j)的凈對流通量，近似為 <mrow><mfrac> <mo>∂</mo> <mrow><mo>∂</mo><mi>y</mi> </mrow></mfrac><mrow> <mo>(</mo> <msub><mi>v</mi><mi>y</mi> </msub> <mi>C</mi> <mo>)</mo></mrow><mo>≈</mo><mfrac> <mn>1</mn> <mi>Δy</mi></mfrac><mo>{</mo><msub> <mi>v</mi> <mrow><mi>y</mi><mrow> <mo>(</mo> <mi>i</mi> <mo>+</mo> <mn>1</mn> <mo>/</mo> <mn>2</mn> <mo>,</mo> <mi>j</mi> <mo>)</mo></mrow> </mrow></msub><mo>[</mo><mrow> <mo>(</mo> <mn>1</mn> <mo>-</mo> <mi>α</mi> <mo>)</mo></mrow><msub> <mi>C</mi> <mrow><mi>i</mi><mo>,</mo><mi>j</mi> </mrow></msub><mo>+</mo><mi>α</mi><msub> <mi>C</mi> <mrow><mi>i</mi><mo>+</mo><mn>1</mn><mo>,</mo><mi>j</mi> </mrow></msub><mo>]</mo><mo>-</mo><mo>-</mo><mo>-</mo><mo>-</mo><mrow> <mo>(</mo> <mn>8</mn> <mo>)</mo></mrow> </mrow> <mrow><msub> <mi>v</mi> <mrow><mi>y</mi><mrow> <mo>(</mo> <mi>i</mi> <mo>-</mo> <mn>1</mn> <mo>/</mo> <mn>2</mn> <mo>,</mo> <mi>j</mi> <mo>)</mo></mrow> </mrow></msub><mo>[</mo><mrow> <mo>(</mo> <mn>1</mn> <mo>-</mo> <mi>α</mi> <mo>)</mo></mrow><msub> <mi>C</mi> <mrow><mi>i</mi><mo>-</mo><mn>1</mn><mo>,</mo><mi>j</mi> </mrow></msub><mo>+</mo><mi>α</mi><msub> <mi>C</mi> <mrow><mi>i</mi><mo>,</mo><mi>j</mi> </mrow></msub><mo>]</mo><mo>}</mo> </mrow>其中，α為空間加權因子；中心加權，α等于0.5；采用上風加權方案公式(2)右邊第七項是由源/匯引起的流入或流出單元(i，j)凈質(zhì)量通量；公式(2)右邊第八項和第九項是單元(i，j)內(nèi)一級化學反應引起的質(zhì)量損失或增加量；吸附濃度表示為溶解濃度的函數(shù)，具體取決于所用的吸附等溫線；公式(2)的左邊是儲存于單元(i，j)的質(zhì)量變化率，近似為 <mrow><mi>R</mi><mfrac> <mrow><mo>∂</mo><mi>C</mi> </mrow> <mrow><mo>∂</mo><mi>t</mi> </mrow></mfrac><mo>≈</mo><mi>R</mi><mfrac> <mrow><msubsup> <mi>C</mi> <mrow><mi>i</mi><mo>,</mo><mi>j</mi> </mrow> <mrow><mi>n</mi><mo>+</mo><mn>1</mn> </mrow></msubsup><mo>-</mo><msubsup> <mi>C</mi> <mrow><mi>i</mi><mo>,</mo><mi>j</mi> </mrow> <mi>n</mi></msubsup> </mrow> <mi>Δt</mi></mfrac><mo>-</mo><mo>-</mo><mo>-</mo><mrow> <mo>(</mo> <mn>9</mn> <mo>)</mo></mrow> </mrow>其中延遲因子R由吸附等溫線確定；6)對以上各分量差分計算過程詳細實現(xiàn)過程如下①彌散系數(shù)D由邊界流速分量及公式(10)，利用IDL矩陣運算計算彌散系數(shù)分量； <mrow><msub> <mi>D</mi> <mrow><mi>xx</mi><mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>±</mo> <mn>1</mn> <mo>/</mo> <mn>2</mn> <mo>)</mo></mrow> </mrow></msub><mo>=</mo><mfrac> <mrow><msub> <mi>α</mi> <mi>L</mi></msub><msubsup> <mi>v</mi> <mrow><mi>x</mi><mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>±</mo> <mn>1</mn> <mo>/</mo> <mn>2</mn> <mo>)</mo></mrow> </mrow> <mn>2</mn></msubsup> </mrow> <msub><mi>v</mi><mrow> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>±</mo> <mn>1</mn> <mo>/</mo> <mn>2</mn></mrow> </msub></mfrac><mo>+</mo><mfrac> <mrow><msub> <mi>α</mi> <mi>T</mi></msub><msubsup> <mi>v</mi> <mrow><mi>y</mi><mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>±</mo> <mn>1</mn> <mo>/</mo> <mn>2</mn> <mo>)</mo></mrow> </mrow> <mn>2</mn></msubsup> </mrow> <msub><mi>v</mi><mrow> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>±</mo> <mn>1</mn> <mo>/</mo> <mn>2</mn></mrow> </msub></mfrac><mo>+</mo><msup> <mi>D</mi> <mo>*</mo></msup> </mrow> <mrow><msub> <mi>D</mi> <mrow><mi>yy</mi><mrow> <mo>(</mo> <mi>i</mi> <mo>±</mo> <mn>1</mn> <mo>/</mo> <mn>2</mn> <mo>,</mo> <mi>j</mi> <mo>)</mo></mrow> </mrow></msub><mo>=</mo><mfrac> <mrow><msub> <mi>α</mi> <mi>L</mi></msub><msubsup> <mi>v</mi> <mrow><mi>y</mi><mrow> <mo>(</mo> <mi>i</mi> <mo>±</mo> <mn>1</mn> <mo>/</mo> <mn>2</mn> <mo>,</mo> <mi>j</mi> <mo>)</mo></mrow> </mrow> <mn>2</mn></msubsup> </mrow> <msub><mi>v</mi><mrow> <mi>i</mi> <mo>±</mo> <mn>1</mn> <mo>/</mo> <mn>2</mn> <mo>,</mo> <mi>j</mi></mrow> </msub></mfrac><mo>+</mo><mfrac> <mrow><msub> <mi>α</mi> <mi>T</mi></msub><msubsup> <mi>v</mi> <mrow><mi>x</mi><mrow> <mo>(</mo> <mi>i</mi> <mo>±</mo> <mn>1</mn> <mo>/</mo> <mn>2</mn> <mo>,</mo> <mi>j</mi> <mo>)</mo></mrow> </mrow> <mn>2</mn></msubsup> </mrow> <msub><mi>v</mi><mrow> <mi>i</mi> <mo>±</mo> <mn>1</mn> <mo>/</mo> <mn>2</mn> <mo>,</mo> <mi>j</mi></mrow> </msub></mfrac><mo>+</mo><msup> <mi>D</mi> <mo>*</mo></msup><mo>-</mo><mo>-</mo><mo>-</mo><mrow> <mo>(</mo> <mn>10</mn> <mo>)</mo></mrow> </mrow> <mrow><msub> <mi>D</mi> <mrow><mi>xy</mi><mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>±</mo> <mn>1</mn> <mo>/</mo> <mn>2</mn> <mo>)</mo></mrow> </mrow></msub><mo>=</mo><mrow> <mo>(</mo> <msub><mi>α</mi><mi>L</mi> </msub> <mo>-</mo> <msub><mi>α</mi><mi>T</mi> </msub> <mo>)</mo></mrow><mfrac> <mrow><msub> <mi>v</mi> <mrow><mi>x</mi><mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>±</mo> <mn>1</mn> <mo>/</mo> <mn>2</mn> <mo>)</mo></mrow> </mrow></msub><msub> <mi>v</mi> <mrow><mi>y</mi><mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>±</mo> <mn>1</mn> <mo>/</mo> <mn>2</mn> <mo>)</mo></mrow> </mrow></msub> </mrow> <msub><mi>v</mi><mrow> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>±</mo> <mn>1</mn> <mo>/</mo> <mn>2</mn></mrow> </msub></mfrac> </mrow> <mrow><msub> <mi>D</mi> <mrow><mi>yx</mi><mrow> <mo>(</mo> <mi>i</mi> <mo>±</mo> <mn>1</mn> <mo>/</mo> <mn>2</mn> <mo>,</mo> <mi>j</mi> <mo>)</mo></mrow> </mrow></msub><mo>=</mo><mrow> <mo>(</mo> <msub><mi>α</mi><mi>L</mi> </msub> <mo>-</mo> <msub><mi>α</mi><mi>T</mi> </msub> <mo>)</mo></mrow><mfrac> <mrow><msub> <mi>v</mi> <mrow><mi>x</mi><mrow> <mo>(</mo> <mi>i</mi> <mo>±</mo> <mn>1</mn> <mo>/</mo> <mn>2</mn> <mo>,</mo> <mi>j</mi> <mo>)</mo></mrow> </mrow></msub><msub> <mi>v</mi> <mrow><mi>y</mi><mrow> <mo>(</mo> <mi>i</mi> <mo>±</mo> <mn>1</mn> <mo>/</mo> <mn>2</mn> <mo>,</mo> <mi>j</mi> <mo>)</mo></mrow> </mrow></msub> </mrow> <msub><mi>v</mi><mrow> <mi>i</mi> <mo>±</mo> <mn>1</mn> <mo>/</mo> <mn>2</mn> <mo>,</mo> <mi>j</mi></mrow> </msub></mfrac> </mrow>其中D*為有效分子擴散系數(shù)； <mrow><msub> <mi>v</mi> <mrow><mo>*</mo><mo>,</mo><mo>*</mo> </mrow></msub><mo>=</mo><msqrt> <msubsup><mi>v</mi><mrow> <mi>x</mi> <mrow><mo>(</mo><mo>*</mo><mo>,</mo><mo>*</mo><mo>)</mo> </mrow></mrow><mn>2</mn> </msubsup> <mo>+</mo> <msubsup><mi>v</mi><mrow> <mi>y</mi> <mrow><mo>(</mo><mo>*</mo><mo>,</mo><mo>*</mo><mo>)</mo> </mrow></mrow><mn>2</mn> </msubsup></msqrt><mo>-</mo><mo>-</mo><mo>-</mo><mrow> <mo>(</mo> <mn>11</mn> <mo>)</mo></mrow> </mrow>②彌散項設置初始溶質(zhì)濃度場，即確定污染事件發(fā)生位置及該物質(zhì)的飽和溶解度作為模型的初始濃度分布矩陣C；由彌散項計算公式(3)、公式(4)、公式(5)及公式(6)，結合彌散系數(shù)計算結果矩陣，利用IDL矩陣運算計算各彌散項；③對流項由對流項計算公式(7)、公式(8)，結合單元格邊界流速分量、土壤彌散度屬性，利用IDL矩陣運算計算各對流項；④公式(2)的第七項，源匯項由于快速處理模型、源匯項設置為0；⑤公式(2)的第八、第九項，化學反應項根據(jù)各化學物質(zhì)的等溫吸附曲線獲??；⑥單元節(jié)點濃度由公式(9)及各個分量矩陣計算結果，利用IDL矩陣運算計算單元節(jié)點后一時刻濃度場；7)返回步驟5)中的計算過程，即實現(xiàn)一個時間段后，區(qū)域流場及溶質(zhì)濃度場的計算；在步驟5)前加一個總時段的循環(huán)控制語句，即實現(xiàn)各時段區(qū)域流場及溶質(zhì)濃度場的計算；8)結果輸出經(jīng)過上述步驟實時輸出地表高程矩陣matrix1、地下水位矩陣matrix2、地下水埋深矩陣matrixh0、地下水流場、溶質(zhì)濃度場和時間項。FSA00000208012600023.tif,FSA00000208012600039.tif,FSA000002080126000310.tif
2.如權利要求1所述的一種地下水環(huán)境的建模及數(shù)值模擬方法，其特征在于所述步 驟4)中，根據(jù)中心加權方案，從節(jié)點濃度獲得界面濃度
全文摘要
本發(fā)明涉及一種地下水環(huán)境的建模及數(shù)值模擬方法，其包括以下步驟1)應用IDL軟件的空間插值函數(shù)KRIG2D()或TRI_SURF()函數(shù)對地表高程、潛水埋深進行空間網(wǎng)格插值，同時劃分空間有限差分數(shù)值模擬網(wǎng)格，得到兩個相同的差分網(wǎng)格矩陣；設置模擬區(qū)土壤屬性文件和地表高程與地下水水位屬性文件；計算地下水位埋深h0和單元邊界流速場；2)建立非飽和帶一維溶質(zhì)運移模型；3)建立飽和帶溶質(zhì)遷移二維數(shù)值模型4)將上式中的各項式進行近似表示5)利用IDL矩陣運算對以上各分量進行求解。本發(fā)明基于IDL平臺，可高效、快速、簡便地模擬地下水環(huán)境，因此，可廣泛用于地下水環(huán)境突發(fā)污染事件的數(shù)值模擬、預測、應急預測等過程中。
文檔編號G06F19/00GK101908100SQ201010238839
公開日2010年12月8日 申請日期2010年7月26日 優(yōu)先權日2010年7月26日
發(fā)明者吳文強, 李偉峰, 陳求穩(wěn), 馬金鋒, 黃國鮮 申請人:中國科學院生態(tài)環(huán)境研究中心