fgModflow计算地下水流程式原理.doc

1.1 Modflow計算地下水流程式原理
　　運用PMWIN(processing modflow)程式中之Modflow來地下水模擬，模式採用有線差分法(FDM)來求得近似數值解。模擬之前需將模擬區域劃分成網格(Cell)，每一網格假設為均勻物質，以網格節點或網格中心之水頭代表網格水頭，含水層厚度在水平方向上可以變化，但在數值上會有較大的誤差。
　　在Modflow中使用的固定密度(Constant Density)地下水流偏微分方程式(McDonald and Harbaugh，1988，p.2-1)如下所示：
　　　　　　式1-1
式中、、分別為沿著x、y、z座標軸的水力傳導係數。
　　h=壓力水頭(Potentiometric Head)
　　W=地下水每單位體積的體積流量(Volumetic Flux)的補注源(Sources)或輸出源(Sinks)。當W<0.0表示流出地下水系統(Out Flow of Ground-Water System)，反之，W>0.0表示流進地下水系統。
　　Ss=孔隙介質之儲水率。
　　t = 時間(T)。
　　當式1-1結合邊界及初始條件時，描述在非均質性及非等向性介質(提供主軸方向的水力傳導係數)中的暫態三維地下水流。左側代表含水層在不同方向之入滲與水頭變化及地下水補給與抽水造成之三向度地下水流動模式；在均向性穩定流含水層系統中，因為流率為一常數，多孔介質之儲水量為定值，此式即為Laplace方程式，可由給定之含水層流況、初始水位狀況與邊界條件獲得正確解。惟在大多數非穩定流系統中，含水層之水位變化在整個模式中是個變數，因此需將水層系統分隔成小區塊，以求得滿足各小區塊間之解與整個系統之近似解；在此假設環境下每一小區塊在單位時間內體積流率即可設定是一常數(McDonald and Harbaugh，1988)。
　　而由連續方程式(Continuity Equation)來發展有限差分形式之地下水流方程式，亦表示網格流入及流出的總和要等於網格中儲存量的變化率。在地下水密度固定的假設下，連續方程式表示網格中地下水流平衡方程式為式1-2所示：
　　　　　　　　式1-2
式中：流入及流出網格之地下水流量之總和
　　　：在有限差分法方程式中的比儲水係數。比儲水係數的定義為每單位含水層結構體積及每單位壓力水頭變化之抽出或補助地下水的變化。
　　 ：網格體積大小。
　　 ：在單位時間間隔中，壓力水頭的變化。
　　在式1-2等號右邊的項表示在已知的水頭變化及的時間間隔中，網格中儲水體積的變化率。因此，式1-2等號的左邊表示流入及流出網格的項及體積總和。
　　而外流力的作用如井水抽水、補注、蒸發、排水、河流入滲等，對網格之流量率影響亦可由下式1-3表示之(McDonald and Harbaugh，1988)：
　　　　　　式1-3
式中：由外界流入網格(i,j,k)之總體積流量率。
　　 ：外流力介質之導水係數。
　　 ：網格中之水位。
　　 ：直接注入網格之流量率。
　　 ｎ　：外流力數。
　　因此整個網格(i,j,k)之三向度水流有限差分連續方程式，即可由網格間之流量率變化與外流力作用比表示成式1-4：
　　　　　　　　　　　　　　　　　　　　　　　　　　　式1-4
　　其中右側之為不同時間之水位變化，亦可以有限差分法表示成式1-5：
　　　　　　　　　　　　　　　　　　　　　式1-5

　　而將單位網格整體地下水流有限差分方程統整如下式1-6
　　　　　　　　　　　　　　　　　　　　　　　　　式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
式中，：極限沉陷量
　　　 ：壓縮係數
　　　 ：初始孔隙比
　　　：初始壓密壓力
　　：有效應力增加量
　　　：土層厚度
　　若土壤之壓密性質或孔隙比隨深度有很大的變化，或截然為不同土層，則總沉陷量S為各土層沉陷量之總和。
　　極限沉陷量求得之後，利用以下壓密方程式，式1-8則可計算歷時性之沉陷量S(t)。
　　　　　　　　　　　　　　　式1-8
式中u：孔隙水壓
　　 t：時間
　　Cv：壓密係數
　　Z：深度
Taylor(1948)對上式壓密方程式，提出一個以Fourier級數展開的方程式來表示的數學精確解，如式1-9所示：
　　　　　　　　　　式1-9
式中：平均壓密度　　　　　　：時間因數

1.3地層下陷機制
　　在含水層系統中儲水量變化的主要因素有三：(1)水位下降使孔隙間排水；(2)地下水位下之含水層(受壓水層)有效應力增加，導致含水層骨架壓縮；(3)孔隙水壓下降導致水膨脹(Jacob，1940)。其中儲水量變化並可由沉積物骨架之比儲水係數與含水層單位時間水頭變量表示。
　　由於地層下陷是因儲存於多孔介質中之水被排除所造成，當排水時土壤架構將產生應變而使儲水能力降低。根據Terzaghi壓密理論與有效應力觀點，進行抽水時，儲存於多孔介質中之水被移除，將造成土壤有效應力增加，若此時有效應力增量未造成含水層中細顆粒沉積物夾層重新排列時，沉積物的壓縮行為屬彈性壓縮，此時若地下水位上升即可使沉積物骨架回復；然而，當有效應力超過彈性壓縮界限時，造成永久壓縮，此時即使有效應力降低(地下水位上升)，沉積物亦僅產生部分回脹現象。永久壓縮行為亦會在有效應力再次超過最近之最大應力階段後再發生，且非彈性壓縮範圍內之有效應力增量造成的單位壓縮量較彈性範圍為大(Jorgenson，1980；Riley，1969)。

1.4地層下陷量
根據Riley(1969)與Helm(1975)研究指出，自由水層沉積物骨架之彈性與非彈性壓縮量或回脹量，與沉積物骨架之單位體積比儲水係數成正比，其關係式為式1-10與式1-11：
　　　　式1-10
式1-11
式中：彈性壓縮量(L)
：沉積物骨架之彈性比儲水係數(L-1)
：非彈性壓縮量(L)
沉積物骨架之非彈性比儲水係數(L-1)
：夾層之厚度(L)
在受壓水層內因總應力為常數,期由有效應力增量造成受壓水層沉積物骨架之壓縮量亦可由式1-12與式1-13之關係式表示，其中有效應力增量以受壓水位變量質替代，可以由下式1-12與1-13表示之：
式1-12
式1-13
式中：水頭變化量(L)
其中沉積物骨架之彈性與非彈性比儲水係數，可由地下水流模式中沉積物骨架之體積流量率及水頭變量(Leake ＆ Prudic，1998)如下式1-14
　　　　　　　 式1-14
式中：地下水流入或流出可壓縮性夾層儲水量之單位體積流量率
　　：夾層之比儲水係數
　　上式1-14表示細顆粒夾層中地下水流量率與比儲水係數之關係，其中之值相對於彈性與非彈性比儲水係數值，將取決於相對水位與沉積層預壓密水位之大小。其中之預壓密水位對應於最大有效應力(即最小壓密水位)，當水位大於預壓密水位時，為彈性比儲水係數；而若水位較預壓密水位低時，為非彈性比儲水係數；兩者將因地下水之水位變化而改變。其中沉積物股價之比儲水係數與沉積物厚度之乘積即為沉積物骨架之儲水係數。
　　細顆粒沉積層之比儲水係數亦可由土壤壓密試驗之壓縮指數性質與有效應力的關係獲得(Jorgenson，1980)：
式1-15
式1-16
其中：細顆粒沉積物再壓縮指數
　　：細顆粒沉積物壓縮指數
　　：有效應力
　　：初使孔隙比
　　根據(Jorgenson，1980)研究指出，在一多夾層且不同儲水係數與導水係數之含水層系統內，系統整體之儲水係數與導水係數可由每單一夾層之儲水係數與導水係數累加決定之。設一系統具有ｎ個細顆粒夾層之比儲水係數Ss1、Ss2、...Ssn，滲透係數為k1、k2、...kn，每一夾層厚度分別為b1、b2、...bn，則系統之儲水係數與導水係數可由下式1-17、1-18獲得：
式1-17
式1-18
地層下陷有限差分方程式
　　地層下陷模式之有限差分方程式係以隱性差分法(imlicit method)從每一時階(time step)完成後之水頭與預壓密的關係式來計算沉積物骨架之體積流率方程式(Leake ＆ Prudic，1988)，如下式1-19：
式1-19
式1-20
式中：在m時階時水流入可壓縮夾層之單位體積流率
　　：m時刻之比儲水係數
　　：時間階段
　　：在m時階完成後網格之水位
　　：在m時階完成後網格之預壓密水位
　　：在m-1時階(上一時階)完成後網格之水位
　　：在m-1時階(上一時階)完成後網格之預壓密水位
　　在式1-20中，若時階完成後之水位上小於預壓密水位時，則選擇非彈性比儲水係數作為下一模擬時階進行計算；相反的，若水位大於預壓密水位時，則選擇彈性比儲水係數進行下一時階計算。
　　上述之含水層夾層壓縮理論，亦可用來計算厚之侷限層或難透水層因水位變化產生的壓縮量；模擬方法系將侷限層或難透水層細分為多層不同夾層計算即可(Leake ＆ Prudic，1988)