李天浩臺灣大學：土木工程學研究所李文生Lee, Wen SenWen SenLee2007-11-252018-07-092007-11-252018-07-092005http://ntur.lib.ntu.edu.tw//handle/246246/50285本研究以混合型理查氏方程式變飽和度地下水流方程式，建構能同時模擬飽和與未飽和地下水流，確保模式能達成局部與全域質量守恆的數值模式。首先，探討一維模式在模擬時遭逢的不收斂或不穩定問題，研發解決的方式，以及利用多層次網格細化方法，提升模式計算效率；其次，建立垂直二維並考慮與坡面正交水平方向坡寬的擬三維變飽和度地下水流模式，測試與驗證滲流面邊界條件的處理方法，改進模式解矩陣的效率，探討坡地形狀對於飽和地下水分佈與滲流面的影響。 對於乾燥土壤給予定流量邊界條件，容易發生數值不穩定問題，本研究比較測試三種入滲邊界數值處理方法：(1)以質量平衡原則迭代水頭邊界條件取代流量邊界條件的方法，(2)利用解析解作為迭代求解水頭的初估值，與(3)使用邊界半差分元素控制體積質量守恆法。研究結果發現，方法(1)在所有狀況下，都可以達到收斂與質量守恆，唯迭代計算需花費較多的時間；方法(2)必須利用土壤保水和導水試驗資料推估解析法專用的特性曲線，可以改善收斂性，唯效率與土壤特性曲線與解析法專用的特性曲線的類似程度有關；方法(3)在大部分狀況下可以收斂，也是三種方法中效率最高者，唯在極端狀況下有不能收斂的可能性。 當定流量上邊界大於土壤飽和水力傳導係數時，初期可全部入滲，當地表土層飽和後開始積水，由於外部邊界條件為給定流量邊界條件，在地表土壤飽和後無法以外部給定的定流量入滲。在此狀況下，本研究採用的邊界條件處理方法，是將定流量上邊界轉換為定水頭，在不考慮逕流下，以積水深度與入滲量的和必須等於降雨量（定流量）的方式，迭代決定積水深度（定水頭邊界條件），以及模式是否達到收斂標準。式模擬試驗證實，在重新處理過邊界條件後，不同降雨條件下，均能夠有效率且確定模式收斂的情形下完成模擬。 在空氣可以自由從側邊排除假設下的一維變飽和數值模擬問題中，地表飽和的入滲濕鋒段與地下水位接觸時，土層瞬間飽和將使土壤水壓於短時間內發生急劇變化，可能導致數值不穩定。本研究對於此問題的數值處理方法，是採用土層達到完全飽和前的入滲率，與飽和前後全域水質量差異，推估時間差分步幅的方式處理，可以滿足質量守恆的要求，並且跳過全域壓力水頭應變數大幅驟變的數值處理困難。 將理查氏方程式比擬為運移－擴散微分方程式，由運移－擴散方程式的係數u/D比值，計算佩克勒數(Peclet number)，Pe，數值試驗結果顯示，當Pe值越大時，數值擴散情形越明顯。本研究利用多層次細化網格增加模式計算效率的做法，是選擇Pe值門檻，當Pe的數值大於某個增加網格的Pe值門檻時，以網格尺寸減半的比率進行網格細化、降低Pe值；（在濕鋒後端）當Pe的數值小於某個刪除網格的Pe值門檻時，刪除局部最細網格，使網格尺寸倍增減少不必要的計算。利用多層次網格模式自動產生與去除、並且相鄰網格大小比率不超過2:1和1:2範圍的設計，可以使模擬區域的Pe值均維持在一定的數值以下。數值模擬顯示，多層次網格模式與全域採用細網格模擬的含水量及壓力水頭都幾乎重合，局部網格細化所得的數值準確性與全域採用細網格的結果相當。多層次網格細化法，除了可減少數值擴散問題，改善入滲濕鋒位置與梯度的模擬外；也增加網格配置上的彈性，利用數量較少的網格，達成全域使用細網格模擬相同的計算精度，縮短計算所需時間，提升模式的執行效率；同時，亦可以適用於分層異質土壤問題。 由於多層次網格細化模式建立在質量守恆模式的基礎上，數值模擬可以達到完全的質量守恆；但是在增減網格時，因為保水曲線曲率的凹凸，可能增加或減少水量，造成質量不守恆問題。另外，由於濕鋒前端預備細化網格的配置是採用已知時刻的Pe值進行研判，因此可能在模擬下一時間步的水流時，濕鋒範圍超過細網格範圍的問題，尚待進一步解決。 擬三維數值模式之處理方式是擴充一維模式，加上沿坡面方向的坡向維度，以及與坡面正交水平方向的坡寬變化，以考量寬度變化的質量守恆式搭配達西定律重新推導得擬三維坡地地下水流方程式。在邊界的處理上則多了滲流面邊界的迭代求解過程。求解單次迭代過程中線性矩陣的方式是採用預設式共軛梯度法，能有效的減少運算資源，提高計算效率。由擬三維模式與實驗資料的驗證結果可發現，擬三維模式對於二維與楔形實驗土槽之土壤水分與滲流面模擬具有相當不錯的掌握能力。 　　利用擬三維數值模式模擬穩態與非穩態降雨和坡地滲流，結果顯示，坡地滲流問題必須採用能模擬未飽和層水的流動與勢能分布的變飽和度數值模式，才能正確銜接土壤－大氣接觸的邊界條件與飽和地下水位之間的物理作用機制，對於滲流面分佈的估計具有相當重要的影響。穩態數值試驗顯示，採用杜布假設、以飽和地下水流模式模擬坡地水流，無法正確表現近下游出口處的地下水位曲率和滲流面，並且可能低估實際流量。在非穩態坡地水流問題中，必須具備未飽和土壤水流數值模式，才能描述水自地表入滲到進入飽和地下水面的流動行為，和降雨停止後的退水行為；說明飽和地下水流模式對於坡地水文行為描述能力有限。 　　由於未飽和層水流主要在垂直方向移動，飽和地下水則橫向移動顯著，欲完整描述坡地之水文反應，坡地水流模式至少需具備垂直二維的模擬能力。擬三維模式利用準單向（沿坡地方向）流動理論與質量平衡觀念，推導出具有束縮或擴張因子的擬三維變飽和度水流控制方程式，可配合應用於卡氏座標與圓柱座標系統問題。透過合理的垂直二維與寬度函數結合，模式可應用於三維坡地之水流模擬，並可節省模式運算資源。 　　坡面的形狀影響地下水位在坡地土壤中的變化，向下束縮的坡面容易出現滲流面，同時在靠近坡面下游處，單寬流量和壓力水頭梯度較大；向下擴展的坡面則較少出現滲流面，在坡面下游擴展處壓力水頭梯度較小，上游處寬度較窄，單寬流量和壓力水頭梯度較大。坡面上下游水深差異亦會影響滲流面的高度，相同的坡地長度與土壤厚度，當上下游水深差異較大時，下游較容易出現滲流面，即使在寬度向下游擴展的坡面亦是。 　　坡面飽和水力傳導係數會影響入滲量的多寡，在降雨強度小於飽和水力傳導係數時，降雨可以全部入滲，反之可能只部分入滲。坡面飽和水力傳導係數較降雨強度小者，短延時降雨的累積入滲量不多，加上入滲速度較慢，因此形成飽和地下水位的時間也較慢，只能形成淺薄的飽和地下水層，過程中土壤有效應力均保持較大，坡地較為穩定。飽和水力傳導係數較大的坡面，降雨會快速入滲，雨停後也能夠快速的排水。在降雨強度大、降雨延時比入滲鋒抵達飽和地下水位時間長的狀況下，坡面可能出現滲流面。The purpose of this study is to develop a Finite-Analytic numerical model to calculate moisture and heat fluxes across the soil-atmosphere interface and the vertical distribution of soil moisture and temperature for the bare soil in the unsaturated zone. With more accurate flux estimates, micrometeorology related and natural ground water recharge predictions are likely to be improved. The study contains the following features: (a) It utilizes the equations for moisture and heat transfer following the works of Chen (1995), and reforms the heat flux at the surface boundary condition for evaporation. (b) It accounts for the temperature effect in the metric head and hydraulic conductivity and uses the experimental results of Hopmans and Dane (1985) to calculate the metric head and the water capacity as functions of temperature. (c) It uses the method of de Vries (1975) to calculate the thermal conductivity as well as the correction factor for the average macroscopic temperature gradient. Both are decided by water content and the volumetric fraction of quartz and clay. (d) It applies a Finite-Analytic Numerical Model for efficient solution of strongly non-linear problems for abrupt changes in atmosphere and soil conditions. (e) The results are tested against field measurements, Jackson (1973) and Jackson et al. (1973). Including surface energy-balance components, water content and soil temperature. According to the results, the simulated and observational trends of the soil water content are difference at begging. With few hours ago, the effect of initial condition is reductive and than the simulated and observational trends of the soil water content is closer to closer. The simulated values of the soil temperature also can response the reasonable physical phenomenon. In the process to calculate the metric head and soil temperature, the convergence and the calculating time of metric head is harder and longer than soil temperature. Otherwise, hydraulic conductivity, soil surface roughness length, and sky cover will influence the varies of soil moisture and soil temperature. The hydraulic conductivity controls and regulates the water movement upward to the surface or downward to the groundwater table where the evaporation or infiltration occurs. The surface roughness length controls the value of mass transfer and sensible heat transfer coefficients, which affects the evaporation rate, latent heat fluxes, and sensible heat flux. The sky cover determines the quantity of solar energy, which can produce significant differences in water content, soil temperature and evaporation rate.目 錄 摘要 i 目錄 iv 表目錄 vi 圖目錄 vii 第一章 緒論 1.1問題概述 1-1 1.2研究目標 1-2 1.3文獻回顧 1-2 1.4研究方法 1-8 1.5論文架構 1-10 第二章 坡地變飽和度地下水流數學模式 2.1問題概述 2-1 2.1.1未飽和層現象 2-1 2.1.2坡地水文探討 2-5 2.2變飽和度水流控制方程式 2-8 2.2.1變飽和度地下水流方程式 2-8 2.2.2坡面水流數學模式 2-12 2.2.3擬三維坡地水流方程式 2-13 2.3初始條件 2-15 2.4邊界條件 2-16 2.4.1地表邊界條件 2-16 2.4.2下邊界條件 2-17 2.4.3側邊界條件 2-17 2.4.4滲流面邊界 2-17 第三章 一維變飽和度地下水流模式與應用 3.0模式應用價值 3-1 3.1數值問題和解決方法 3-1 3.2質量守恆問題與處理 3-2 3.2.1數值模式理論 3-2 3.2.2數值離散與線性化 3-4 3.3初始條件與邊界條件處理 3-6 3.3.1初始條件處理 3-6 3.3.2邊界條件處理 3-7 3.4模式演算流程 3-11 3.5瞬間飽和問題與處理 3-14 3.6多層次網格細化邏輯 3-21 3.6.1理查氏方程之數學物理特性 3-21 3.6.2網格細化之探討及處理方式說明 3-22 3.6.3模式應用與分析 3-25 3.7小結 3-43 第四章 坡地變飽和地下水流模式與應用 4.1質量守恆問題與處理 4-1 4.2邊界條件處理 4-5 4.2.1定水頭邊界條件 4-5 4.2.2定流量邊界條件 4-5 4.2.3滲流面邊界 4-11 4.3求解線性聯立方程組方法 4-12 4.3.1共軛梯度法 4-13 4.3.2預設式共軛梯度法 4-14 4.4模式驗證 4-17 4.4.1垂向入滲問題模擬 4-18 4.4.2穩態楔型水槽滲流模擬 4-19 4.4.3非穩態變飽和度入滲模擬 4-21 4.4.4非穩態變飽和度滲流模擬 4-23 4.5小結 4-25 第五章 坡地滲流問題數值試驗與分析 5.1坡面形狀分類與分析案例概述 5-1 5.2上下游不同水深差異之近穩態水位剖面 5-2 5.2.1 模擬條件說明 5-2 5.2.2 矩形坡面之滲流面與流量分析 5-3 5.2.3不同坡面形狀之滲流面與出流量分析 5-4 5.3非穩態坡地降雨之水分分布分析 5-9 5.3.1模擬條件說明 5-9 5.3.2初始坡地水分分布 5-10 5.3.3坡面土壤透水性較降雨強度大之土壤水分分析 5-11 5.3.4坡面土壤透水性較降雨強度小之土壤水分分析 5-14 5.3.5最大滲流面長度等值分布圖 5-18 第六章 結論與建議 6.1研究結論 6-1 6.2研究建議 6-3 參考文獻 7-1 附錄A A-1 附錄B B-1 表目錄 表3.1 Touma和Vauclin（1986）實驗土壤資料表 3-14 表3.2上邊界定流量8.3cm/hr各模擬案例之比較列表 3-28 表3.3上邊界定流量20cm/hr各模擬案例之比較列表 3-33 表3.4上邊界定水頭2.3cm各模擬案例之比較列表 3-36 表3.5網格1cm和5cm之零流量穩態水分分佈質量表 3-42 表4.1 Haverkamp等人(1977)土壤參數表 4-18 表4.2 Hall(1955)之土壤參數表 4-19 表4.3 Vauclin(1979)之土壤參數表 4-21 表5.1模擬坡面條件列表 5-2 表5.2 S1~S9九種坡面形狀在上下游水頭差為50公分時 (下游150公分，上游200公分)的出口流量比較 5-5 表5.3 S1~S9九種坡面形狀在上下游水頭差為150公分時 (下游50公分，上游200公分)的出口流量比較 5-7 表5.4九種坡面形狀在不同上下游水深差異下，近穩態之 滲流面長度表 5-8 圖目錄 圖1.1研究流程圖 1-8 圖2.1土壤孔隙示意圖 2-1 圖2.2毛細管壓力分布示意圖 2-3 圖2.3土壤進水、退水保水曲線差異示意圖 2-3 圖2.4土壤遲滯現象之物理機制示意圖 2-4 圖2.5遲滯現象中，保水曲線示意圖 2-5 圖2.6降雨時期，變源面積之變化圖 2-6 圖2.7石門水庫上游秀巒站單元集水區 2-7 圖2.8不同坡地型態最大可蓄水容積示意圖 2-7 圖2.9坡地土壤之微小元素水頭分配關係 2-12 圖2.10擬三維斜坡地之微小元素之通量示意圖 2-14 圖3.1各節點權重係數簡圖 3-6 圖3.2模擬區域控制體積示意圖 3-8 圖3.3邊界格點之控制體積示意圖 3-10 圖3.4模式運算流程圖 3-12 圖3.5 Touma和Vauclin（1986）實驗土壤水力特性曲線圖 3-15 圖3.6定流量上邊界qub=8.3cm/hr持續入滲使地下水位上漲 示意圖 3-15 圖3.7模擬降雨強度8.3cm/hr，延時3.6小時之土壤入滲 含水量(左圖)，與壓力水頭(右圖)剖面分布圖 3-17 圖3.8模擬降雨強度8.3cm/hr，延時3.6小時之地表入滲率 歷線 3-17 圖3.9模擬降雨強度8.3cm/hr，延時3.6小時之地表水頭 歷線(上圖)，與地下水位歷線(下圖) 3-18 圖3.10模擬降雨強度20cm/hr，延時1.5小時之土壤入滲 含水量(左圖)，與壓力水頭(右圖)剖面分布圖 3-19 圖3.11模擬降雨強度20cm/hr，延時1.5小時之地表入滲 率歷線 3-19 圖3.12模擬降雨強度20 cm/hr，延時1.5小時之地表水頭 歷線(上圖)，與地下水位歷線(下圖) 3-20 圖3.13模式計算含水量垂直分佈與對應之 比值關係圖 3-22 圖3.14單一濕鋒之u/D值分佈與多層次網格判斷示意圖 3-25 圖3.15定流量8.3cm/hr之多層次網格細化模組與固定網 格模式模擬比較圖 3-27 圖3.16多層次網格模式隨入滲濕鋒移動之細網格配置 3-28 圖3.17定流量20.0cm/hr之多層次網格細化模式與固定網 格模式模擬比較圖 3-31 圖3.18原始網格1cm多層次網格細化模組之u/D值與 貝克勒數分佈 3-32 圖3.19上邊界定流量20cm/hr下，固定網格0.1cm，固定 網格10cm，與多層次網格原始網格10cm之計算 地表水頭時間序列 3-32 圖3.20定水頭2.3cm之多層次網格細化模式與固定網格 模式模擬比較圖 3-35 圖3.21定水頭2.3cm固定網格0.1cm，固定網格10cm， 及多層次網格模式之計算地表入滲率時間序列 3-36 圖3.22多層次網格細化模組在模擬歷程中，因為增刪網 格點造成水質量差異的累積值 3-38 圖3.23多層次網格模式模擬上邊界定流量8.3cm/hr入滲， 對網格重新進行分配的過程中，增加網格點與刪除 網格對水質量造成的影響 3-39 圖3.24多層次網格模式模擬上邊界定流量20cm/hr入滲， 對網格重新進行分配的過程中，增加網格點與刪除 網格對水質量造成的影響 3-40 圖3.25多層次網格模式模擬上邊界定水頭2.3cm入滲， 對網格重新進行分配的過程中，增加網格點與刪除 網格對水質量造成的影響 3-41 圖3.26均勻格點1cm和5cm下，零流量穩態條件下的含水 量分佈 3-42 圖4.1二維坡地土壤空間離散格點示意圖 4-2 圖4.2 x-z平面的模擬區域中，邊界控制體積示意圖 4-6 圖4.3利用擬三維模式模擬一維垂向入滲水流之模擬結果 4-19 圖4.4利用擬三維模式模擬楔型土槽穩態滲流面高度之 模擬結果 4-20 圖4.5利用擬三維模式模擬二維具補注地下水之飽和水面 與未飽和水分分布之模擬結果 4-23 圖4.6利用擬三維模式模擬二維非穩態滲流面之模擬結果 4-24 圖5.1坡面形狀分類圖 5-1 圖5.2矩形坡面(S7)在上下游水頭固定情形下，24小時後 的壓力水頭分布圖 5-3 圖5.3矩形坡面(S7)在上下游水頭固定情形下，24小時內 的入出流量歷線圖 5-3 圖5.4兩端水深差異50公分(下游150公分，上游200公分)， 近穩態時各坡面之壓力水頭分布圖 5-5 圖5.5兩端水深差異150公分(下游50公分，上游200公分)， 近穩態時各坡面之壓力水頭分布圖 5-7 圖5.6設計降雨之降雨組體圖 5-10 圖5.7初始壓力水頭分布圖 5-10 圖5.8 Imax/Ks=0.25時，矩形坡面S7在模擬24小時內的 壓力水頭分布圖 5-12 圖5.9 Imax/Ks=0.25時，矩形坡面S2_1在模擬24小時內的 壓力水頭分布圖 5-13 圖5.10 Imax/Ks=0.25時，矩形坡面S5_1在模擬24小時內的 壓力水頭分布圖 5-14 圖5.11 Imax/Ks=4時，矩形坡面S7在模擬24小時內的 壓力水頭分布圖 5-16 圖5.12 Imax/Ks=4時，寬度向下束縮坡面S2_1在模擬24 小時內的壓力水頭分布圖 5-17 圖5.13 Imax/Ks=4時，寬度向下擴展坡面S5_1在模擬24 小時內的壓力水頭分布圖 5-18 圖5.14最大滲流面位置等值圖 5-19 圖B.1最陡坡降法延負梯度方向收斂到唯一解的迭代過程 B-3 圖B.2利用正交向量修正迭代解直接求解之示意圖 B-4 圖B.3 Gram-Schmidt共軛對概念示意圖 B-53970848 bytesapplication/pdfen-US理查氏方程式變飽和度網格細化擬三維滲流面variably saturated flowgrid refinementquasi-3Dseepage facethesishttp://ntur.lib.ntu.edu.tw/bitstream/246246/50285/1/ntu-94-D88521015-1.pdf