fgModflow计算地下水流程式原理

1 1 Modflow 計算地下水流程式原理計算地下水流程式原理 運用 PMWIN processing modflow 程式中之 Modflow 來地下水 模擬 模式採用有線差分法 FDM 來求得近似數值解 模擬之前需 將模擬區域劃分成網格 Cell 每一網格假設為均勻物質 以網格 節點或網格中心之水頭代表網格水頭 含水層厚度在水平方向上可 以變化 但在數值上會有較大的誤差 在 Modflow 中使用的固定密度 Constant Density 地下水流偏微 分方程式 McDonald and Harbaugh 1988 p 2 1 如下所示 式 1 1 t h SW z h K zy h K yx h K x szzyyxx 式中 分別為沿著 x y z 座標軸的水力傳導係數 xx K yy K zz K h 壓力水頭 Potentiometric Head W 地下水每單位體積的體積流量 Volumetic Flux 的補注源 Sources 或輸出源 Sinks 當 W0 0 表示流進地下水系 統 Ss 孔隙介質之儲水率 t 時間 T 當式 1 1 結合邊界及初始條件時 描述在非均質性及非等向性 介質 提供主軸方向的水力傳導係數 中的暫態三維地下水流 左側 代表含水層在不同方向之入滲與水頭變化及地下水補給與抽水造成 之三向度地下水流動模式 在均向性穩定流含水層系統中 因為流 率為一常數 多孔介質之儲水量為定值 此式即為 Laplace 方程式 可由給定之含水層流況 初始水位狀況與邊界條件獲得正確解 惟 在大多數非穩定流系統中 含水層之水位變化在整個模式中是個變 數 因此需將水層系統分隔成小區塊 以求得滿足各小區塊間之解 與整個系統之近似解 在此假設環境下每一小區塊在單位時間內體 積流率即可設定是一常數 McDonald and Harbaugh 1988 而由連續方程式 Continuity Equation 來發展有限差分形式之地 下水流方程式 亦表示網格流入及流出的總和要等於網格中儲存量 的變化率 在地下水密度固定的假設下 連續方程式表示網格中地 下水流平衡方程式為式 1 2 所示 式 1 2 V t h SQ si 式中 流入及流出網格之地下水流量之總和 i Q 在有限差分法方程式中的比儲水係數 比儲水係數的定 s S 義為每單位含水層結構體積及每單位壓力水頭變化之抽出或補助地 下水的變化 網格體積大小 V 在單位時間間隔中 壓力水頭的變化 h t 在式 1 2 等號右邊的項表示在已知的水頭變化及的時間間h t 隔中 網格中儲水體積的變化率 因此 式 1 2 等號的左邊表示流 入及流出網格的項及體積總和 而外流力的作用如井水抽水 補注 蒸發 排水 河流入滲等 對網格之流量率影響亦可由下式 1 3 表示之 McDonald and Harbaugh 1988 式 1 3 nkjikjinkjinkji QhPQs 式中 由外界流入網格 i j k 之總體積流量率 nkji Qs 外流力介質之導水係數 nkji P 網格中之水位 kji h 直接注入網格之流量率 nkji Q 外流力數 因此整個網格 i j k 之三向度水流有限差分連續方程式 即可由 網格間之流量率變化與外流力作用比表示成式 1 4 t h vrcSsQsQQQQQQ kji kjikjikjikjikjikjikjikjikji 2 1 2 1 2 1 2 1 2 1 2 1 式 1 4 其中右側之為不同時間之水位變化 亦可以有限差分法表 t h kji 示成式 1 5 式 1 1 1 mm mkjimkjikji tt hh t h 5 而將單位網格整體地下水流有限差分方程統整如下式 1 6 1 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 mm mkjimkji kjikjikji m kjihhvhhv hhrhhrhhchhc tt hh vrcSsQhPTT TTTT kji m kji m kjikji m kji m kjikji m kji m kjikji m kji m kjikji m kji m kjikji m kji m kjikji 式 1 6 式 1 6 即為 McDonald 與 Harbaugh 1988 發展之三向度地下水 流有限差分連續方程式 可由起始水頭與邊界條件以反覆迭代 Interative Method 的方式 來得到每個時間階段 Time Step 之有限 差分方程式系統的解 在每個時間階段的計算期間 初始壓力水頭 分佈會在迭代法的迭代過程中 連續地產生許多中間的壓力水頭分 佈 1 2 地層下陷分析方法地層下陷分析方法 因地下水位下降造成之地層下陷 其下陷量大小與土壤種類有 關 砂質土層的沉陷速度較快但沉陷量較小 而黏性土層之壓縮則 較具延滯性且沉陷量較大 因地下水位下降所引致的地表壓密沉陷之分析方法大致可分為 兩類 一為基於 Terzaghi 壓密理論之 分離式方法 decoupled approach 與另一為基於 Biot 理論之 耦合式方法 coupled approach 前者是先求出孔隙水壓力之分佈 再應用有效應力觀念 計算土層之應變及地表壓密沉陷量 後者是基於孔隙水與介質之間 存在某種互制關係 以介質位移和孔隙水壓力為基本變數之耦合壓 密理論 Terzaghi 壓密理論 本模式以 Terzaghi 1943 之單向壓密理論分析地層下陷問題 該理論係利用飽和土壤內孔隙水在穩定層流情況下之連續條件 及 有效應力原理和土壤本身之應力 應變關係推得計算式 其基本之 假設為 1 土壤是均質且於飽和狀態 2 土壤顆粒與水的壓縮可忽略 3 孔隙水流符合達西定律 Darcy s law 4 土壤之壓縮性與滲透性 在受壓過程中保持不變 5 小載重增量作用基本上沒有產生厚度的 變化 亦即小應變 且壓力不影響土壤特性 即壓縮係數及水力傳導 係數於壓密過程中均為常數 依據 Terzaghi 建議 正常壓密土壤之極限沉陷量之計算 如下 式 1 7 所示 式 1 7 log 1 0 00 0 P P H e C S vc u 式中 極限沉陷量 u S 壓縮係數 c C 初始孔隙比 0 e 初始壓密壓力 0 P 有效應力增加量 0v 土層厚度H 若土壤之壓密性質或孔隙比隨深度有很大的變化 或截然為不 同土層 則總沉陷量 S 為各土層沉陷量之總和 極限沉陷量求得之後 利用以下壓密方程式 式 1 8 則可計算 歷時性之沉陷量 S t 式 1 8 2 2 z u C t u v 式中 u 孔隙水壓 t 時間 Cv 壓密係數 Z 深度 Taylor 1948 對上式壓密方程式 提出一個以 Fourier 級數展開的方 程式來表示的數學精確解 如式 1 9 所示 式 1 9 v n T n n tU 4 12 exp 12 8 1 22 0 22 式中 平均壓密度 時間因數 tU v T 1 3 地層下陷機制地層下陷機制 在含水層系統中儲水量變化的主要因素有三 1 水位下降使孔 隙間排水 2 地下水位下之含水層 受壓水層 有效應力增加 導致 含水層骨架壓縮 3 孔隙水壓下降導致水膨脹 Jacob 1940 其 中儲水量變化並可由沉積物骨架之比儲水係數與含水層單位時間水 頭變量表示 由於地層下陷是因儲存於多孔介質中之水被排除所造成 當 排水時土壤架構將產生應變而使儲水能力降低 根據 Terzaghi 壓密 理論與有效應力觀點 進行抽水時 儲存於多孔介質中之水被移除 將造成土壤有效應力增加 若此時有效應力增量未造成含水層中細 顆粒沉積物夾層重新排列時 沉積物的壓縮行為屬彈性壓縮 此時 若地下水位上升即可使沉積物骨架回復 然而 當有效應力超過彈 性壓縮界限時 造成永久壓縮 此時即使有效應力降低 地下水位上 升 沉積物亦僅產生部分回脹現象 永久壓縮行為亦會在有效應力 再次超過最近之最大應力階段後再發生 且非彈性壓縮範圍內之有 效應力增量造成的單位壓縮量較彈性範圍為大 Jorgenson 1980 Riley 1969 1 4 地層下陷量地層下陷量 根據 Riley 1969 與 Helm 1975 研究指出 自由水層沉積物骨架 之彈性與非彈性壓縮量或回脹量 與沉積物骨架之單位體積比儲水 係數成正比 其關係式為式 1 10 與式 1 11 式 1 100 bS p b ske 式 1 110 bS p b skv 式中 彈性壓縮量 L b 沉積物骨架之彈性比儲水係數 L 1 ske S 非彈性壓縮量 L b 沉積物骨架之非彈性比儲水係數 L 1 skv S 夾層之厚度 L 0 b 在受壓水層內因總應力為常數 期由有效應力增量造成受壓水層沉積 物骨架之壓縮量亦可由式 1 12 與式 1 13 之關係式表示 其中有效 應力增量以受壓水位變量質替代 可以由下式 1 12 與 1 13 表示之 式 1 12 0 bhSb ske 式 1 13 0 bhSb skv 式中 水頭變化量 L h 其中沉積物骨架之彈性與非彈性比儲水係數 可由地下水流模式 中沉積物骨架之體積流量率及水頭變量 Leake Prudic 1998 如 下式 1 14 式 1 14 t h Sq ski 式中 地下水流入或流出可壓縮性夾層儲水量之單位體積流量率 i q 夾層之比儲水係數 sk S 上式 1 14 表示細顆粒夾層中地下水流量率與比儲水係數之關係 其中之值相對於彈性與非彈性比儲水係數值 將取決於相對水 sk S 位與沉積層預壓密水位之大小 其中之預壓密水位對應於最大有效 應力 即最小壓密水位 當水位大於預壓密水位時 為彈性比儲 sk S 水係數 而若水位較預壓密水位低時 為非彈性比儲水係數 兩者 將因地下水之水位變化而改變 其中沉積物股價之比儲水係數與沉 積物厚度之乘積即為沉積物骨架之儲水係數 ke S 細顆粒沉積層之比儲水係數亦可由土壤壓密試驗之壓縮指數性 質與有效應力的關係獲得 Jorgenson 1980 式 1 15 1 434 0 0 ep C S r ske 式 1 16 1 434 0 0 ep C S c skv 其中 細顆粒沉積物再壓縮指數 r C 細顆粒沉積物壓縮指數 c C 有效應力 p 初使孔隙比 0 e 根據 Jorgenson 1980 研究指出 在一多夾層且不同儲水係數 與導水係數之含水層系統內 系統整體之儲水係數與導水係數可由 每單一夾層之儲水係數與導水係數累加決定之 設一系統具有 個 細顆粒夾層之比儲水係數 Ss1 Ss2 Ssn 滲透係數為 k1 k2 kn 每一夾層厚度分別為 b1 b2 bn 則系統之儲水 係數與導水係數可由下式 1 17 1 18 獲得 式 1 17 nsnsssystem bSbSbSS 2211 式 1 18 nnsystem bkbkbkT 2211 地層下陷有限差分方程式 地層下陷模式之有限差分方程式係以隱性差分法 imlicit method 從 每一時階 time step 完成後之水頭與預壓密的關係式來計算沉積物 骨架之體積流率方程式 Leake Prudic 1988 如下式 1 19 式 1 19 111 mm m skmm m skm i hH t S Hh t S q 式 1 20 1 1 mm skv mm ske m sk HhS HhS S 式中 在 m 時階時水流入可壓縮夾層之單位體積流率 m i q m 時刻之比儲水係數 m sk S 時間階段t 在 m 時階完成後網格之水位 m h 在 m 時階完