應用重力模擬與地下水數值模式於比流出量推估 －以濁水溪沖積扇為例

國 立 交 通 大 學

土 木 工 程 學 系

碩 士 論 文

應用重力模擬與地下水數值模式於比流出量推估

－以濁水溪沖積扇為例

Integrating a Gravity Simulation and Groundwater

Numerical Modeling on the Calibration of Specific

Yield for Choshui Alluvial Fan

研 究 生 ： 姚又瑜

指導教授 ： 張良正 博士

應用重力模擬與地下水數值模式於比流出量推估

－以濁水溪沖積扇為例

Integrating a Gravity Simulation and Groundwater

Numerical Modeling on the Calibration of Specific

Yield for Choshui Alluvial Fan

研 究 生：姚又瑜

Student：Yu-Yu Yao

指導教授：張良正 博士 Advisor：Dr. Liang C. Chang

國 立 交 通 大 學

土 木 工 程 學 系

碩 士 論 文

A Thesis

Submitted to Department of Civil Engineering

National Chiao Tung University

in Partial Fulfillment of Requirements

for the Degree of

Master of Science

in

Civil Engineering

July 2013

Hsinchu, Taiwan, Republic of China

I

應用重力模擬與地下水數值模式於比流出量推估

－以濁水溪沖積扇為例

學生：姚又瑜 指導教授：張良正 博士 國立交通大學土木工程學系

摘要

地下水補注量為地下水永續管理重要的資訊之一，而比出水量或儲水 係數之不確定性則直接影響地下水補注量推估之正確性，然而現地比出水 量與儲水係數僅能藉由複井抽水試驗取得，必需在試驗場址鑽鑿兩口水井， 一用於抽水，另一則用於觀測，因此鑽探成本較單井試驗高，受到鑽井成 本的限制，造成比出水量之現地試驗數量遠低於所需，使區域性地下水補 注量之推估具有高度不確定性。 有鑑於此，本研究將以濁水溪沖積扇為研究區域，結合現地重力量測、 重力模擬與地下水模擬，重力量測為非侵入式之調查方式，以取代傳統複 井抽水試驗，進行比流出量之調查，降低地下水補注量推估之不確定性。 本研究分別於溪州國小、土庫國中與客厝國小進行重力調查，量測 2012 年五月、八月與十月豐枯水期間的重力變化，量測重力值須先校正海潮、 天文潮、大氣壓力與地層下陷等不同機制之影響，校正後之重力殘差即代 表地下水量變化對重力之影響。其次，應用重力積分公式，搭配 MODFLOW 系統水量變化進行重力模擬，比較量測重力殘差與模擬重力變化幅度之比 值，再據以調整比出水率參數值，反覆進行上述過程，直至模擬重力變化 與重力殘差變化幅度吻合，即完成參數調查。 調查結果如下，由於溪州(1)為溪州國小鄰近之地下水觀測站，其曾進

II 行複井抽水試驗，如複井試驗之結果如可與重力調查結果吻合，即表示重 力量測與抽水試驗之高度關聯性。研究結果顯示，抽水試驗結果與本研究 調查結果吻合，充分驗證本研究方法之可行性與正確性。在土庫國中與客 厝國小方面，原量測重力殘差變化幅度對模擬重力變化幅度之比值分別為 1.8 與 50 倍，經將原參數值分別調大 1.8 倍與 50 倍後，模擬重力則與量測 重力殘差吻合。經參數調整後，土庫國中鄰近區域之年淨抽水量由 5.285(百 萬噸)提高至 7.232(百萬噸)，而客厝國小則由 11.827 (百萬噸)改變至-3.584 (百萬噸)。 本研究成果顯示應用重力現地量測，配合地下水模擬與重力模擬，可 取代傳統之複井抽水試驗，取得現地比出水量，相對於傳統需鑿設複井， 花費大量鑽鑿成本，重力量測屬於非侵入式的地球物理量測方法，可節省 大量成本，對於降低地下水參數及相關推估值之不確定，有極大助益。

III

Integrating a Gravity Simulation and Groundwater

Numerical Modeling on the Calibration of Specific

Yield for Choshui Alluvial Fan

Student：Yu-Yu Yao Advisor：Dr. Liang-Cheng Chang Department of Civil Engineering

National Chiao Tung University

Abstract

In Taiwan, groundwater resources play a vital role on the regional supply management. Because the groundwater resources have been used without proper management in decades, several kinds of natural hazards, such as land subsidence, have been occurred. The Choshui alluvial fan is one of the hot spots in Taiwan. For sustainable management, accurately estimation of recharge is the most important information. The accuracy is highly related to the uncertainty of specific yield (Sy). Besides, because the value of Sy should be tested via a multi-well pumping test, the installation cost for the multi-well system limits the number of field tests. Therefore, the low spatial density of field test for Sy makes the estimation of recharge contains high uncertainty.

Gravity measurement is a kind of physical geographical approach. Because gravity is a function of the mass of a material and the inverse square of the distance, the gravity measurement has the ability to observe the mass variation of the shallow groundwater system. Using the groundwater level observation and gravity measurement can be used for the calibration of Sy for a numerical model. Because the approach can observe the mass variation of groundwater without well installation, it can overcome the limitation of well construction for the Sy test. Four steps are used in the study. First, gravity

IV

variations of Si-jhou, Tu-ku and Ke-cuo are observed in May, August and November 2012. The measured values are required to remove the effect of other sources, such as ocean tide and land subsidence, and the removed values are called as gravity residual. Second, a numerical model, MODFLOW 2000, has been applied to simulate and calibrate the groundwater variation basing on the four values of Sy obtained from pumping tests. Although the groundwater levels are well-calibrated in the first version of MODFLOW with low spatial density of field test for Sy, the mass variation of groundwater might contain high uncertainty. Third, integrating the Newton gravity integral with the simulated water mass which obtained from the first version of MODFLOW may simulate the gravity variation of the system. Fourth, comparing the ratio between the variation ranges of gravity residual and the one of simulated gravity, the values of Sy can be modified basing on the ratio and assigned in the second version of MODFLOW.

Because a multi-well pumping test has been applied on Si-jhou station, the values of Sy is 0.216. the gravity residual observed from Si-jhou and the measured Sy can be used to verify the proposed approach. Assigning the field test value of Sy in Si-jhou, the simulated gravity can fit the gravity residual well without parameter calibration. The study of Si-jhou is an evidence to verify the ability of the proposed approach. In Tu-ku and Ke-cuo station, the ratios between the variation ranges of gravity residuals and the ones of simulated gravities approximate 1.8 and 50. The values of Sy are modified as 1.8 and 50 times of the values assigned in the first version of MODFLOW. After the parameter re-assignment, the simulated gravities may fit the gravity residuals better in Tu-ku and Ke-cuo. The comparison of water balance between the two versions of MODFLOW indicates that the quantities of recharge for Tu-ku from -5.285 million tons increase to -7.232 million tons, and for Ke-cuo from 11.827 million tons to -3.584 million tons.

V

The result show that using gravity measurement, groundwater simulation model and the one of simulated gravity, to replace the traditional multi-well pumping test, gravity measurement is a kind of non-intrusion physical

geographical approach. Gravity measurement can saving a lot of cost, and also decrease the groundwater parameter uncertainty.

VI

謝誌

耶！我畢業了！首先，感謝我的指導教授 張良正老師，無論是在研究 領域上的耐心指導，以及做人處事的諄諄教誨，對於我們都有很大的影響。 另外，也要特別感謝 黃金維老師的協助與支持，因為老師的鼓勵使我對於 研究更有信心。 當然，一定不會忘記向最敬愛的陳文哥大大致謝，由於你的全力支持 與教導，促使我不斷前進，在你的身上我學到『只要有心，就會成功』！ 你在美國一定會發光發熱的啦！還有我們 419 研究室的夥伴們，彬哥、祐 誠、雲直、逸儒、逸美、怡瑄、深惠、小肥瑜、阿布、呱呱、ㄅㄅ、國陞、 劉鎧、海倫、奕璋、阿文、曉芸、韋炫、又田、誠胤、修緯…從一開始的 漆彈大戰、一起修課趕作業、參加教育訓練、被計畫追著跑、去美國研討 會大吃大喝灑錢看秀、一起唱歌吃飯…這些點點滴滴都是最美好的回憶， 謝謝你們的陪伴，讓研究生活更加精彩豐盛。在此還要特別感謝黃金維老 師團隊的大力幫忙，鄭博、程博、紫猗、冠宏、小光…謝謝你們總是義不 容辭的解決我在研究上的困難，讓對於重力方面一竅不通的我，能夠快速 吸收理解。 最後，要感謝家人無條件的付出與關愛，無論做任何決定總是在身邊 支持鼓勵，每次回家都會補充滿滿的元氣！喔耶～永遠愛你們。

VII

目錄

摘要 ... I Abstract ... III 謝誌 ... VI 目錄 ... VII 第一章 前言... 1 1.1 研究動機 ... 1 1.2 研究目的 ... 2 1.3 研究流程 ... 2 第二章 文獻回顧 ... 4 2.1 地下水數值模式 ... 4 2.1.1 地下水數值模式 ... 4 2.1.2 參數檢定 ... 6 2.2 重力試驗 ... 8 第三章 理論敘述 ... 10 3.1 地下水數值模擬與參數檢定 ... 10 3.1.1 MODFLOW 控制方程式 ... 10 3.1.2 參數檢定 ... 10 3.2 重力試驗與重力積分 ... 15 3.2.1 重力量測原理、設備 ... 15 3.2.2 應用重力積分技術於地下水負載 ... 17 第四章 結合重力量測於濁水溪沖積扇比流出率檢定 ... 19 4.1 研究區域與相關現地參數概述 ... 19 4.2 地下水數值模式初步建置與檢定 ... 23

VIII 4.2.1 地下水分層架構、邊界條件及網格劃分 ... 23 4.2.2 模式輸入資料 ... 30 4.3 現地重力量測結果 ... 38 4.4 應用重力變化於濁水溪沖積扇比流出率之檢定 ... 45 4.4.1 重力積分半徑影響分析 ... 45 4.4.2 比流出率檢定 ... 47 4.4.3 比出水率檢定成果 ... 58 4.4.4 專家系統參數檢定 ... 64 4.5 水平衡分析 ... 66 第五章 結論與建議 ... 68 5.1 結論 ... 68 5.2 建議 ... 70 參考文獻 ... 71 附錄 A 環境改正 ... 76

IX

表目錄

表 1.3-1 濁水溪沖積扇補注量相關研究比較 ... 5 表 4.1-1 濁水溪沖積扇儲水係數或比流出量現地試驗值 ... 22 表 4.2-1 濁水溪沖積扇淺層觀測井列表(33 口井) ... 35 表 4.3-1 絕對重力儀施測地點 ... 42 表 4.3-2 溪州國小量測結果 ... 42 表 4.3-3 土庫國中量測結果 ... 42 表 4.3-4 客厝國小量測結果 ... 42 表 4.5-1 各區域修正前後淨補注量 ... 66

X

圖目錄

圖 1.1-1 研究流程圖 ... 3

圖 3.1-2 參數檢定流程圖 ... 11

圖 3.2-1 FG5 絕對重力儀與各部件明細[Micro-g Lacoste FG5 Manual] ... 17

圖 4.1-1 濁水溪沖積扇(海園-石榴)水文地質剖面圖 ... 20 圖 4.1-2 濁水溪沖積扇地下水觀測網站井分布圖(繪製於民國 98 年) ... 21 圖 4.2-1 濁水溪沖積扇概念分層 ... 24 圖 4.2-2 濁水溪沖積扇地下水邊界之概念圖 ... 25 圖 4.2-3 濁水溪沖積扇之模擬範圍 ... 27 圖 4.2-4 濁水溪沖積扇模式第 1 分層格網與邊界條件 ... 28 圖 4.2-5 濁水溪沖積扇模式第 2~7 分層格網與邊界條件 ... 29 圖 4.2-6 濁水溪沖積扇模式 AA’段剖面圖 ... 30 圖 4.2-7 含水層一徐昇氏分區 ... 31 圖 4.2-8 含水層二徐昇氏分區 ... 32 圖 4.2-9 含水層三徐昇氏分區 ... 32 圖 4.2-10 淺層含水層種類與已知儲水係數之分佈 ... 34 圖 4.3-1 絕對重力儀(FG5)架設於客厝國小 ... 39 圖 4.3-2 GPS 監測站架設於客厝國小 ... 40 圖 4.3-3 重力測站於濁水溪沖積扇分佈位置 ... 41 圖 4.3-4 GPS 連續監測地表高程變化曲線－客厝國小 ... 43 圖 4.3-5 GPS 連續監測地表高程變化曲線－土庫國中 ... 43 圖 4.3-6 GPS 連續監測地表高程變化曲線－溪州國小 ... 44 圖 4.4-1 影響半徑分析圖 ... 46 圖 4.4-2 加入虛擬觀測站之徐昇式網格分區 ... 49

XI 圖 4.4-3 客厝國小內插水位與鄰近三站觀測水位之變化 ... 50 圖 4.4-4 土庫國中內插水位與鄰近三站觀測水位之變化 ... 50 圖 4.4-5 比流出量或儲水係數分布圖(修改前) ... 52 圖 4.4-6 模擬水位振幅分布圖(修改前) ... 53 圖 4.4-7 模擬重力振幅分布圖(修改前) ... 54 圖 4.4-8 溪州國小實測重力與模擬重力變化圖 ... 56 圖 4.4-9 土庫國中實測重力與模擬重力變化圖 ... 57 圖 4.4-10 客厝國小實測重力與模擬重力變化圖 ... 57 圖 4.4-11 調整範圍示意圖 ... 58 圖 4.4-12 土庫國中實測重力與模擬重力變化圖(修正後) ... 60 圖 4.4-13 客厝國小實測重力與模擬重力變化圖(修正後) ... 60 圖 4.4-14 比流出量或儲水係數分布圖(修改後) ... 61 圖 4.4-15 模擬水位振幅分布圖(修改後) ... 62 圖 4.4-16 模擬重力振幅分布圖(修改後) ... 63 圖 4.4-17 觀測水位與檢定水位檢定關係圖(修改前) ... 64 圖 4.4-18 觀測水位與檢定水位檢定關係圖(修改後) ... 65 圖 4.5-1 土庫國中修正前後淨補注量比較圖 ... 67 圖 4.5-2 客厝國小修正前後淨補注量比較圖 ... 67

第一章 前言

1.1 研究動機

台灣地區蘊含豐富的地下水資源，在區域供水上扮演重要的角色，但是許 多區域地下水因未經妥善管理，過度抽水導致嚴重之環境衝擊，例如地層下陷 或海水入侵等，濁水溪沖積扇為台灣地區最重要的地下水區之一，區域大量使 用地下水，且部分區域已因過度抽水，而有嚴重地層下陷問題。 地下水資源之永續利用需基於收支平衡之觀點，因此在進行地下水資源之 管理，了解年地下水補注量為地下水管理最重要的資訊。傳統上多以地下水數 值模式推估地下水補注量，搭配現地觀測之觀測水位變化，模式藉由模擬地下 水水位變化，進一步推估系統之補注量或抽水量，掌握系統之目前收支狀況。 水文地質參數的給定為地下水模擬最重要的一環，參數不確定性將影響補注量 或抽水量推估之不確定性，其中又以比出水量或儲水係數最為直接。 比出水量或儲水係數可用以描述多孔介質中，地下水水位與系統水量之關 係，故如參數數值如具有不確定性，其直接影響後續之補注量與抽水量之推估。 由於比出水量與儲水係數僅能藉由複井抽水試驗取得，因此需在鄰近區域鑽鑿 兩口水井。以透水係數而言，由於其僅需要單井即可進行試驗，故地下水觀測 站網於設置觀測井時，必藉由抽水試驗取得現地透水係數。對於比流出量或儲 水係數而言，則受到複井試驗鑽探成本之限制，其現地試驗值數量遠低透水係 數。以濁水溪沖積扇為例，合計共建置 203 口觀測井，全數僅有 10 口進行複井 抽水試驗，僅有 6 口非侷限水井有比出水率，分別為合興(1)、溪州(1)、柑園(1)、 二水、烏塗(2)與東光(1)；僅有 4 口拘限水井有儲水係數，分別為線西(1)、海園 (1)、田洋(1)與莿桐(1)觀測井。由於現地試驗值大幅短缺，使得後續補注量與抽 水量之推估，具有高度不確定性，進一步影響地下水之永續管理。 重力量測為地球物理的量測方法，現地重力受到地球中各種物體質量變化

之影響，包含海潮、天文潮、大氣壓力、地殼變形與地下水系統等。重力公式 為物體質量與距離倒數平方之函數，因此重力現地量測可以呈現一定範圍內地 下水系統之質量變化，相對於傳統需藉由侵入式之井體進行複井抽水試驗，重 力量測直接以非侵入式的方式進行調查，可有效克服建置成本之限制。

1.2 研究目的

本研究將建置地下水數值模式，推估區域性淨補注量，並藉由時序重力量 測資訊，降低比出水率之參數不確定性，以提高濁水溪沖積扇補注量與抽水量 推估之精確度。

1.3 研究流程

本研究之研究流程如圖 1.1-1 所示，主要分為四大步驟，首先應用地下水數 值模式(MODFLOW 2000)，搭配少量之比出水率(Sy)現地試驗值，初步建置與模 擬濁水溪沖積扇地下水系統之水量變化；其次以觀測水位搭配專家系統檢定抽 水量，進而模擬地下水系統的蓄水量變化；接著應用重力積分公式，搭配前述 MODFLOW 系統之蓄水量變化進行重力模擬，比較現地重力殘差與模擬重力之 變化幅度；最後依據量測重力與模擬重力變化幅度之差異，調整系統參數，完 成比出水率之檢定。

第二章 文獻回顧

本研究將結合現地重力量測，其原理為運用重力積分公式，呈現一定範圍 內地下水系統之質量變化，藉由非侵入性的方式進行地質參數調查。配合建立 地下水數值模式，以數值試驗的方式，推估地下水系統之比出水率。

2.1 地下水數值模式

2.1.1 地下水數值模式 關於濁水溪沖積扇地下水系統之抽補強度研究甚多，如劉聰桂(1996)利用熱 核爆氚示蹤方法評估地下水補注量；交通大學(2005)以一維垂向之溼地入滲係數 和旱地降雨入滲率，評估地下水補注量；；姜儷安、歐陽湘（1997）、葉文工(1998)、 林再興(1998) 則利用二維多層地下水流數值模式逆推地下水收支；李振誥(1999) 利用基流資料估記法推估濁水溪沖積扇地下水補注量；江崇榮等人(2006)、地調 所及交通大學(2010)、徐年盛、江崇榮、汪中和等(2012)則以水位歷線法推估地 下水補注及抽取；上述利用各種研究方法，推估濁水溪沖積扇歷年來之地下水 補注量，約於 0.35 至 1.26(百萬噸/平方公里)之間，而總補注量約在 3.76 至 23.05(億噸/年)不等，上述重要文獻數據整理至表 1.3-1。

表 1.3-1 濁水溪沖積扇補注量相關研究比較 方法 估算人 單位面積補 注量 (百萬噸/平 方公里) 面積 總補注量 (億噸/年) 碳十四定年與氚示 蹤 劉聰桂 （1996） 0.50 濁水溪沖積扇 /約 1800 平方公里 9 3DFEWA 數值模擬 張信誠、劉振 宇(1996) 0.44 雲林地區 /約 1120 平方公里 4.88 MODFLOW 數值模擬 姜儷安、歐陽 湘 （1997） 0.42 雲林地區 /約 1120 平方公里 4.66 MODFLOW 數值模擬 葉文工 （1998） 0.5 濁水溪沖積扇 /約 1800 平方公里 8.97 非飽和層一維 長期水文模式 劉振宇 （1998） 0.59 至0.61 濁水溪沖積扇 /約 1800 平方公里 10.06 至 10.48 水平衡法 陳進發等 (1998) 0.35 至0.55 彰化地區 /1074平方公 里 3.76至 5.91 MODFLOW 數值模 擬 林再興(1998) 0.97 彰化地區 /約 830 平方公里 8.07 基流資料估計法 李振誥 （1999） 0.52 濁水溪沖積扇 /約 1800 平方公里 9.3 MODFLOW 數值模 擬 張良正、劉振 宇(2002) 1.26 雲林地區 /約 1120 平方公里 14.11 地下水歷線分析法 (未考慮未降雨時期 之邊界地下水補注 和地下水灌溉補注) 江崇榮等 (2006) 0.62 濁水溪沖積扇 /2079 平方公里 12.89 地下水歷線分析法 交通大學 (2010) 0.88 濁水溪沖積扇 /2079 平方公里 18.28

地下水歷線分析法 (同時考慮雨水補 注、河水補注、邊界 地下水補注和地下 水灌溉補注) 徐年盛、江崇 榮、汪中和 (2012) 0.90 濁水溪沖積扇 /2562 平方公里 23.05 2.1.2 參數檢定 如前所述，地下水數值模式常應用在許多地下水管理問題中，惟地下水數 值模式的建置需大量之水文地質參數，如透水係數(Hydraulic conductivity (K)或 儲水係數(Storage coefficient (S)，以及其他抽水量、補注量與邊界條件等資料， 這些參數常難以直接現地量測，或因成本之關係其資料密度常遠低於模式所需， 因此模式建置過程中，常需以參數檢定方式，逆向推求部分參數數值，此一般 稱為參數檢定(Yeh, 1985)。 參數檢定最常用的方式為藉由調整模式參數值而使模式模擬值接近觀測資 料值，而參數調整又可分為人工調整或是藉由演算法由電腦程式調整。人工檢 定最大的問題為費時費工，為了解決這個問題，許多研究利用優選法進行模式 之自動參數檢定(Mazi et al., 2004; Mazi et al., 2000; Hill et al., 1992)。在電腦普及 計算速度突飛猛進的現在，自動參數檢定妥善利用此優點，達到相較於人工率 定省時省力的目的。然而，應用優選法於參數檢定需先確認並建立目標函數與 限制式，亦即需先將如何調整參數的想法結構化成數學方程式，因此相較於人 工參數檢定較缺乏彈性。傳統上，採梯度類型優選法之自動化參數檢定工具， 如 Ucode 等，由於梯度類型之優選方法需以參數與水位之敏感度矩陣訂定搜尋 方向，為使參數檢定模式通用化，在此常以差分方式建立敏感度矩陣。然差分 化敏感度矩陣，需藉由反覆呼叫模擬模式方可求得，故其計算量隨著參數維度 大幅增加，對於高參數維度問題之計算量將大幅增加。且其梯度型演算法受到 初始搜尋位置的影響極劇，若初始解不同常造成不收斂或不一致之結果。 人工參數檢定之優點為模式建立者可藉由參數檢定過程，增加對模式的了

解，修正模式架構上可能的誤差。 Boyle (2000)比較了水文模式之人工參數檢定 及應用優選法之參數檢定，並且合併兩種參數檢定方法以改善兩者之缺點。他 比較人工參數檢定與應用優選法之參數檢定所得到之結果，指出前者較能被水 文學家所接受，因為人工參數檢定之過程不僅可以被檢視，而且能幫助水文學 家進行結果分析。惟人工參數檢定之缺點為參數推估過程複雜而冗長，費時費 力，且其經驗與知識多累積於人，因此如將模式知識傳承轉移給他人，也是應 用人工參數檢定上的一大問題(Madsen et al., 2001; Chau and Chen, 2001; Chau, 2004; Chau, 2006)。

有鑑於人工參數檢定的專業知識轉移困難，以往有許多研究應用專家系統 協助進行參數檢定(Abbott, 1991; Chau and Chen, 2001; Chau and Albermani, 2002, 2003; Chau, 2004; Kim, 2007)，以求兼顧知識透明度與檢定效率。專家系統是一 個仿人類專家之智慧型電腦系統，使用者可以透過與專家系統的互動介面來描 述問題，專家系統再藉由其內部已定義好的知識與推理來回答問題。因此專家 系統的建置，即是在將專家的知識進行結構化及表達，亦即在將專家知識進行 透明化。陳(2008)以專家系統應用於穩態及暫態地下水流模式之參數檢定，其著 重於檢定模式之透水係數(K)及儲水係數(Ss 與 Sy)。王(2010)與張(2011)則延續陳 (2008)之架構，分別應用專家系統於穩態與暫態淨補注量(net Recharge, Q)之檢定 上，兩者除以設計案例驗證系統正確性外，均延伸應用於濁水溪沖積扇上。研 究指出，以專家系統為基礎之參數檢定系統，除可克服因差分帶來之大量計算 量外，對不同初始解而言，均可有類似之檢定結果，顯現該參數檢定系統之強 健性。本研究將延續張(2011)之成果，應用於濁水溪沖積扇之參數檢定上。

2.2 重力試驗

基於牛頓重力理論，萬有引力是為鄰近質量除以距離平方之函數，因此重 力場的時空變化受到自然界中許多機制的影響，包含潮汐、大氣壓力、固體潮、 地層變形與地下水的水文變化等之影響。

對於地下水而言，當地下水體水量的增減，將影響鄰近區域的重力場。Tapley et al. (2004) 與 Wahr et al. (2004) 指出短期的重力變化，主要受到土壤中蓄水量 變化的影響。因此，藉由地下水位之觀測值，可用以推估地下水蓄水量的變化， 進而移除飽和含水層與土壤含水量對重力之影響(Rodell and Famiglietti, 2002; Rodell et al., 2007) 。前述研究以水文觀測資料為基礎，進而移除地下水等水文 變化對重力之影響。Hasan and Troch (2006)應用超導重力儀量測重力場之時序變 化，該研究指出重力殘差變化，與飽和水位、未飽和含水量及降雨量均具有高 度相關，其中又未飽和含水量相關性最高。

在參數推估方面，國外以重力量測應用於儲水係數或比流出量之推估上， Pool and Eychaner (1995)、Pool(2008)與 Gehman(2009)等研究均以無限延伸之布 格平板，來近似模擬現地地下水蓄水質量變化對重力之影響，再以此模型反演 系統之比流出量(specific yield, sy) 等水文地質參數。在布格平板的假設下，計算 上可大幅簡化，但由於布格平板係假設一個無限延伸的均質薄平板，對於實際 地下水問題而言，地下水系統存在高度異質性，且亦地下水系統亦有其邊界， 而非無限延伸之問題，因此在此假設下，可能引入參數檢定上之不確定性。 USGS (1997)亦在美國亞利桑納州進行重力量測，於 1992 年 12 月至 1994 年 1 月共 5 個時間點量測重力值，藉由豐枯水期間水位變化，搭配重力場之變 化，探討系統儲水係數。研究指出若岩性屬於河道沉積物，儲水係數將介於 0.15 ~ 0.34 之間；若岩性較為膠結，儲水係數則介於 0.07~0.18 之間。此外，USGS (2003) 於美國加州再次進行現地注水試驗，藉由注水行為對於地下水水位之影響，以 重力變化與水位變化，研究指出該區域儲水係數約為 0.13，單場注水可使得重

力增加 66( micro-gals)。Christiansen 等人(2011)在研究區域同時量測地下水位與 重力場的時空變化，同時藉由兩種觀測資料，建立系統自動參數檢定系統。該 研究指出，相對於傳統僅以地下水位進行參數檢定，同時考量重力資訊可使參 數檢定結果更敏感與準確。 綜合上述研究，顯示應用重力量測資料推估儲水係數為一可行之方法。因 此，本研究將參考前述文獻，建置地下水數值模式，推估區域性淨補注量，並 藉由時序重力量測資訊，降低比出水率之參數不確定性，以提高濁水溪沖積扇 補注量與抽水量推估之精確度。

第三章 理論敘述

本研究係結合地下水數值模擬與重力量測技術，進行比出水量之檢定，本 章理論概述分別討論地下水數值模擬與參數檢定，以及重力試驗與重力積分兩 者。

3.1 地下水數值模擬與參數檢定

第 3.1.1 小節主要說明 MODFLOW 之地下水流控制方程式，當使用者完成 網格切割、邊界條件與參數設定後，模式即可依控制方程式求得模擬水位。但 部分參數缺乏紀錄或無現地調查，因此難以完成正向模擬，因此於第 3.1.2 說明 應用逆問題的建置方式，進行系統之參數檢定。 3.1.1 MODFLOW 控制方程式 MODFLOW 核心所考量之，三維飽和地下水流可以下列之偏微分方程式來 表示： t h S W z h k z y h k y x h k x xx yy zz s                                     ………...……….(3.1.1) 其中 zz yy xx K K K , , _：沿主軸 Z Y X, , 方向的透水係數(Hydraulic Conductivity)(LT1) h：總水頭(Potentiometric Head)(L）

W：單位體積的體積流率(Volumetric Flux)，代表源匯項 (Sources/Sinks)(T1) Ss：孔隙介質的比儲水量(Specific Storage)（L-1）

t：表時間(T)

3.1.2 參數檢定

地下水流模式之正向模擬是指，使用者需先輸入水文地質參數、抽水量或 補注量等，數值模式即可求得模擬水位，藉由輸入正確參數，以求得模擬水位

之流程，可稱為正向模擬。然部分參數因為缺乏紀錄、無現地調查或參數不確 定性大等原因，如以此殘缺不全或具有高度不確定性之參數值輸入模式，模式 求得之模擬水位將與觀測水位有極大之誤差。 為確保數值模式可正確模擬地下水系統，可以觀測水位為目標，反覆迭代 調整前述參數。如觀測水位與模擬水位間之誤差可降至可接受之檢定範圍內， 該組參數即可視為正確值，此即稱為參數檢定。由於係以觀測水位與模擬水位 之誤差最小為目標，藉此逆向推估參數值，亦可稱為逆問題(inverse problem)。 本研究係採用王雲直(2009)與張弼舜(2010)開發之專家系統參數檢定系統， 系統是以專家系統為檢定核心，擷取參數檢定專家之知識，搭配系統推論法則 進行參數檢定。系統檢定流程如圖 3.1-2 所示，在參數檢定開始後，系統首先初 始化參數值，即淨補注量及抽水量，並讀入相關參數的設定資訊，如分區設定、 觀測井位置、觀測水位等。 開始參數檢定 初始化參數值與 設定模式參數檢定 相關資訊 執行地下水模式 MODFLOW 傳送相關資訊 至專家系統 專家系統對各區 進行參數修正 以修正值更新 MODFLOW輸入檔 各區檢定誤差 均小於容許誤差 是 停止參數檢定 否 圖 3.1-2 參數檢定流程圖

在完成初始參數設定後，接著執行地下水模式 MODFLOW ，MODFLOW 完成地下水流模擬後，將地下水位輸出成 ASCII 格式之檔案。 下一步則為專家系統前處理器判斷各分區之觀測水位與模擬水位之誤差(以 下稱為檢定誤差)是否均小於容許誤差。當有一分區之檢定誤差超過容許誤差， 專家系統前處理器則會取出各分區的相關資訊，包括水位高、淨補注量或抽水 量、模擬水位與觀測水位誤差以及上次參數檢定所得參數修正量，然後傳入專 家系統進行推論分析。 將所需資訊傳送至專家系統後，系統接下來會執行專家系統，專家系統依 據知識庫中之規則集合(rules)，並搭配推理機(inference engine)之運作，而推論出 各區的參數修正量。 當各區待檢定之參數，由專家系統分析出修正量後，專家系統後處理器再 將各修正量依據 MODFLOW 輸入檔格式，對 MODFLOW 輸入檔進行更新。 接著再次執行 MODFLOW，計算出更新參數後的地下水位，再由專家系統 前處理器判斷各分區檢定誤差是否均小於容許誤差，若是則停止參數檢定，若 否則反覆執行上述步驟，直至所有分區檢定誤差均小於容許誤差。  參數檢定資訊系統平台 專家系統乃一模仿人類專家且具有決策能力的智慧型電腦系統，本研究採 用 CLIPS(C Language Integrated Production System)建立地下水模式參數檢定專 家系統，CLIPS 乃一方便採用物件導向概念，用以建立專家系統之工具。在本 研究中，地下水模式以分區為概念進行參數檢定，一個分區即為一個物件，在 CLIPS 中則以建立實例(instance)的方式應用物件導向概念。在每個實例(物件)中， 則有若干屬性(slot)代表各分區的特性，如 K 值、所在分層、淨補注量或抽水量、 模擬水位與觀測水位誤差及各種參數調整依據。研究所建立之專家系統，乃將 熟悉地下水模式參數檢定之人類專家所具有之經驗及知識轉換成規則(Rule)並 將之存放於知識庫(Knowledge Base)，當專家系統對地下水模式進行參數檢定時，

專家系統依據所輸入之情況並透過推理機(Inference Engine)對知識庫內之規則 進行規則推論(Rule Inference)，其推理結果即為所輸入問題之答案。  參數檢定規則與推論 在建立專家系統之前，必須先擷取參數檢定之經驗及知識。由於本研究室 之研究團隊，在地下水檢定方面已有相當多研究，相關經驗豐富，因此與研究 團隊之學長以及指導教授進行訪談及討論，之後再依據討論結果，歸納出以下 幾點概念： 一般地下水流模式中，由於受到流線的影響，使得上游進行調整時，會對 下游影響較大；而下游進行調整時，對上游影響較小。因此模式大多會從上游 開始往下游調整。 1. 在多層地下水模式檢定中，由於表層補注量的大小對於下層含水層抽水 量有顯著的影響，因此進行人工檢定時，通常會先進行表層抽水/(補注)量檢定， 再依序往下層進行抽水量檢定。 2. 在人工檢定地下水模式時，通常會以迭代方式來回調整各分區參數，使 各區之檢定誤差逐漸縮小。 基於上述概念，本研究經由反覆測試，將其歸納成更具體的參數調整原則， 詳細說明如下所示： 1. 在進行參數調整前，各分區應先設定初始修正量(進行參，作為參數調整 幅度之基準。 2. 開始調整時，需判別修正方向。當觀測水位大於模擬水位時，往正方向 調整，即增加水量；當觀測水位小於模擬水位時，往負方向調整，此時應減少 水量。 3. 判別修正方向後，接著應判斷修正量，其判斷方法為比較本次檢定與上 次檢定時之修正方向。若修正方向相同，表示參數修正持續向同方向進行，則 修正量不變；若修正方向相反，表示參數調整震盪，則修正量減半。

4. 經多次調整後，參數修正量會逐漸縮小，然而因多分區同時進行調整， 各分區會受相鄰分區影響，可能使得調整量縮小過快，造成檢定速度趨緩，因 此當檢定速度過慢時，各分區調整量設定為同時放大一定倍數。  參數檢定初始設定 首先，專家系統在進行初次檢定時，會先設定各分區參數調整量Q。之後 每次參數調整，Q都會根據規則進行調整，並儲存起來作為下一次參數調整之 參考。各分區均有自己的Q，根據規則進行調整，不受其它分區影響。此外， 每次進行推論會得出一個Q之外，還會獲得一個正號或負號。正號即代表該區 的實際的水量應比現在要多，因此需增加水量，所以將現在的抽水量(或補注量) n Q 加上所推論出的Qn，而獲得新的抽水量(或補注量) Qn1，如式 3.1-1。反之 亦然，負號即代表該區的水量應比現在要少，因此需減少水量，所以將現在的 抽水量或補注量 n Q 減去所推論出的Qn，而獲得新的抽水量或補注量Qn1，如式 3.1-2。 n n n Q Q Q 1   ……….………(3.1-1) n n n Q Q Q 1   _{………..………(3.1-2)}

3.2 重力試驗與重力積分

本研究藉由絕對重力量測儀器，於研究區域內，分別於豐枯兩季量測地球 重力場變化，搭配豐枯水季間水位變化以及隨之而來的重力變化，檢定研究區 域含水層之儲水係數(S, Storage coeffcience)或比出水量(Sy, Specific yield)。

3.2.1 重力量測原理、設備 重力測量即於某一點位置量取其重力加速度的大小。本文所討論的絕對重 力測量原理，主要是以自由落體法則來測定。假設有某一物體由一高度 h 自由 落下至 h0，則落下距離 h-h0與落下時間 t 的關係式為： h-h0=v0t+(1/2)gt2 即為 h=h0+v0t+(1/2)gt 2 ………(3.2-1) 上式(3.2-1)中 g 為重力加速度值，h0表示於時間 t=0 物體所在位置，h 則是 時間為 t 時刻物體所在位置，v0表示 t=0 時物體的速度。根據上式，只要量測不 同落下時間 ti及與其對應的落下距離 hi-h0，則可求解出重力加速度值 g。由於式 (3.2-1)中有三個未知數(h0、v0與 g)，因此必須測定至少三組的 hi及 ti值以組成聯 立方程式求解出 g 值。本實驗所使用的 FG5 絕對重力儀，即是採用如上觀測原 理。 舉例而言，自由落體在三個位置上的落下時間及距離分別為(t1, h1), (t2,h2)及 (t3,h3)，再假設 d1=h1-h0 d2=h2-h0 d3=h3-h0 依據(3.2-1)式，可得：

d1 = v0t1+(1/2)gt1 2 d2 = v0t2+(1/2)gt2 2 d3 = v0t3+(1/2)gt32 以上三式中，將第二式及第三式均減去第一式，可得： (d2-d1)=(t2-t1)[v0+(1/2)g(t2+t1)] (d3-d1)=(t3-t1)[v0+(1/2)g(t3+t1)] 令 D1=d2-d1，D2=d3-d1，T1=t2-t1，T2=t3-t1，代入上式後，消去 v0，簡化後則 可得 g=2×[(D2/T2)-(D1/T1)] / (T2-T1) ………..………....(3.2-2) 上式中的 g 為重力加速度值，亦即為本計畫欲測定之絕對重力值。以 FG5 絕對重力儀為例，該儀器即是以碘穩頻雷射來測定式(3.2-2)中的落下距離 D，並 以銣原子鐘測定落下時間 T，觀測精確度可達 10-8 ms-2，亦即 μgal 等級。 本實驗使用儀器為內政部採購之絕對重力儀 FG5。絕對重力儀的源起為 Hammond 和 Faller 於 1997 年在美國國家標準局合作下，並得到空軍劍橋研究 室(Air Force Cambridge Research Laboratories, AFCEL)的資助，於 Wesleyan 大學 成功研製了可移動式的絕對重力儀 FG-5，此乃是自十九世紀可移動、可倒擺重 力儀之後的第一台可移動式絕對重力儀，如圖 3.2-1 所示。FG5 主要部份是一台 極精準的碘穩頻雷射干涉儀，用來監測物體自由落下的直角稜角運動，干涉計 所收集到的光學資料就提供了一套非常精確的距離觀測量系統，可以追蹤出絕 對重力波長的標準。簡單來說，其原理就式以及精密的測距系統搭配精密計時 的原子鐘，精準地測定物體落下的重力加速度值，經由大量地收集落下的數據， 最後以電腦軟體程式作業計算得出一平均加速度值，即為絕對重力值。

圖 3.2-1 FG5 絕對重力儀與各部件明細[Micro-g Lacoste FG5 Manual] 3.2.2 應用重力積分技術於地下水負載

前述已經說明重力理論與 FG5 絕對重力儀，對於所量測的重力數據，過去 常以布格平板來解析重力數據與其他地球物理之間的關係(Heiskanen and Moritz 1993)，布格平板係以均質、無限延伸之有限厚度平板來概念描述與簡化地球對 重力量測點之關係，均質且無限延伸之假設，使得此數學問題較易求解分析。 在均質且無限延伸平板之假設下，對於具異質且非無限延伸的問題而言，布格 平板則不敷使用。 目前由於資訊運算技術大幅興盛，各式數值方法蓬勃發展，應用數值積分 技術可無需搭配均質且無限延伸，應用數值積分極適合應用於高度異質性之問 題。Sato 與 Hanada (1984)在其海潮重力負載模式(GOTIC)中，除建立海潮之數 值模擬外，透過數值網格可呈現海潮之空間異質變化，透過數值積分方式，計 算在此海潮變化下之重力負載。

對於地下水系統而言，受到各地沖積條件不同，因此土壤具備高度異質性， 各地區之儲水係數或比流出量均有大幅差異，對於地下水數值模擬而言，可藉

由輸入不同之儲水係數或比流出量，可模擬地下水系統之空間異質變化。地下 水數值模式，可將地下水系統切割為有限個數的矩形區塊，即為數值模式之網 格，重力模擬積分中可採用相同之網格，或是再進一步加密切細。 任一重力積分的矩形區塊均由八個角點所構成，其座標可描述為(xi,yj,zk)， 其中下標 i、j 與 k 等於 1 或 2。另外，重力量測點之座標標示為(x0,y0,z0)，以重 力公式搭配數值積分，該矩形區塊中的質量對於量測點於 z 方向之重力分量， 可寫為式 3.2-3 與 3.2-4 所示。



x y z



d zd yd x z G A x x y y z z z             

  

   2 2 1 2 1 1 3 2 2 2

     







       2 1 2 1 2 1 , , 1 1 1 i j k k j i k j i z y x N    ………...…………...(3.2-3) 其中



x_i y_j z_k



G



x_i



D y_j



 y_j



D x_i



 N , ,  log   log 



k k i k z z x D z     1   tan 2 ….……….……….……...(3.2-4) 其中，令δxi=xi-x0、δyj=yj-y0、δzk=zk-z0與 D=δxi2+δyj2+δzk2。G 代表重力常 數，其數值為 _{( )；ρ 代表區塊內之密度。式 3.2-3 與 3.2-4} 可推估數值模式中的單一網格對於量測點之 z 方向重力分量，因此累計所有網 格對於量測點之重力分量，即為整個系統對於量測點的重力分量。 在評估地下水量變化對模擬重力之貢獻時，可以地下水模式之網格為基礎， 進一步加密切割積分網格，以增加積分精度。由於模式網格較粗，故加密之積 分網格則採模式網格之模擬水位與水量，以此計算網格內之水量密度，再代入 式 3.2-3 與 3.2-4 計算重力貢獻值。

第四章 結合重力量測於濁水溪沖積扇比流出率檢定

由於過去儲水係數現地試驗值十分缺乏，參數不確定性將導致地下水補注 量評估之不確定性。有鑒於此，本研究先應用目前稀少之儲水係數值，搭配 MODFLOW 2000 地下水數值模式，初步建置、檢定與模擬濁水溪沖積扇之地下 水變化。接著，利用重力積分公式，搭配前述 MODFLOW 推估之模擬水量變化 進行重力模擬。最後，依據現地重力量測值與模擬重力值之變化幅度比值，檢 定研究區域內含水層之儲水係數(S, Storage coefficient)或比流出量(Sy, Specific yield)。

4.1 研究區域與相關現地參數概述

濁水溪沖積扇地下水區位於臺灣西部海岸之中段，北起烏溪，南至北港溪 南岸，東以八卦台地及斗六丘陵山脊線為界，並沿兩丘陵間之隘口向東擴及竹 山、名間一帶，西臨臺灣海峽，面積約 1,800 平方公里。本沖積扇最主要之溪流 為濁水溪，其主流發源自中央山脈西翼、合歡山以南與玉山山脈北側之間地區， 向西切穿雪山－玉山帶、麓山帶與八卦－斗六丘陵區，在丘陵區西側形成沖積 扇系統，主流流經沖積扇之中央，向西注入臺灣海峽。其他位於沖積扇之河川， 自北而南尚有舊濁水溪、新虎尾溪、舊虎尾溪及北港溪等。濁水溪沖積扇地勢 平緩，高程約介於海拔 0~100 公尺之間。 一、水文地質架構 中央地調所於民國 88 年完成臺灣地區地下水觀測網第一期計畫濁水溪沖積 扇水文地質調查研究報告，並以濁水溪沖積扇 72 站之地層柱狀圖，完成水文地 質剖面 1 至 12，其深度達 300 公尺左右，以及丘陵及河谷區之 8 站地層柱狀圖 繪製水文地質剖面 13 至 15，其深度約達 250 公尺左右;並以這些剖面為基礎劃 分出濁水溪沖積扇概念分層，包括含水層一(F1)、阻水層一(T1)、含水層二(F2)、 阻水層二(T2)、含水層三(F3)、阻水層三(T3)、含水層四(F4)及阻水層四(T4)。上

述各分層之垂直示意圖如圖 4.1-1 所示。 圖 4.1-1 濁水溪沖積扇(海園-石榴)水文地質剖面圖 二、地下水觀測 濁水溪沖積扇地區上游河川與相關之地質條件不利建置水庫，仰賴地下水 水資源比例偏重，所以於第一階段（81 年度至 87 年度）最先實施。目前計有地 下水觀測站 82 站（203 口），抽水試驗站 28 站（55 口），水文地質調查站 88 站。 圖 4.1-2 為濁水溪沖積扇地下水觀測網分布圖。

圖 4.1-2 濁水溪沖積扇地下水觀測網站井分布圖(繪製於民國 98 年) 三、水文地質參數

濁水溪沖積扇大致可區分為扇頂、扇央及扇尾區，扇頂區位於員林、溪州、 西螺、虎溪、東和連線以東，各含水層間無明顥之阻水層存在，地面水可直接 補注至深層，即所謂濁水溪沖積扇扇頂非受壓含水層區；扇央區位於扇頂區以

西，至好修、趙甲、潭墘、田洋、北港聯線以東，阻水層一覆於含水層一之上， 各含水層間有明顯之阻水層存在，地面水無法直接補注各含水層，唯以地層材 料而言，砂礫材料所佔之比例較大；扇尾區位於扇央區以西至沿海，含水層一 以上亦有阻水層覆蓋，各含水層間亦如扇央區皆有明顯阻水層存在，唯其地層 材料中，砂礫所佔之比例較小，而粉砂、泥、及粘土所佔之比例較大。透水係 數 K 介於 10-3~10-5(公尺/秒)之間，導水係數 T 介於 0.01～4.19 平方公尺/分鐘 之間，比容量 Q/s 介於 0.34~208.33 立方公尺/小時/公尺，含水層水力特性以扇 頂最佳，扇央次之，扇尾最差。 表 4.1-1 為濁水溪沖積扇之儲水係數或比流出量現地試驗值列表，由於儲水 係數或比出水率需經由複井抽水試驗取得，受到複井建置成本之限制，其試驗 數量遠低於所需，直至目前為止僅有 4 口井擁有儲水係數，分別為線西(1)、海 園(1)、田洋(1)與莿桐(1)；僅有 6 口井擁有比出水量，分別為合興(1)、溪州(1)、 柑園(1)、二水(1)、烏塗(2)與東光(1)。 表 4.1-1 濁水溪沖積扇儲水係數或比流出量現地試驗值 井名 X Y 儲水係數 S 比出水量 Sy 線西(1) 195062 2669966 0.00448 - 合興(1) 194000 2643600 - 0.225 溪洲(1) 199206 2639069 - 0.216 柑園(1) 201662 2635696 - 0.12 二水(1) 210242 2634428 - 0.199 海園(1) 165467 2624551 0.000191 - 田洋(1) 178717 2624989 0.00188 - 莿桐(1) 197896 2628388 0.0007074 - 烏塗(2) 208141 2629528 - 0.152 東光(1) 174905 2616870 - 0.289

4.2 地下水數值模式初步建置與檢定

4.2.1 地下水分層架構、邊界條件及網格劃分 一、水文地質與模式分層架構 MODFLOW 模擬的設定上，可將含水層型態設定為拘限含水層或非拘限含 水層，並需輸入各分層之上部和底部之高程值。本研究參考中央地質調查所於 民國 88 年出版之「台灣地區地下水觀測網第一期計畫濁水溪沖積扇水文地質調 查研究報告」，以濁水溪沖積扇 72 站地層柱狀圖，完成平原地區水文地質剖面 一至十二（深度至 300 公尺左右），再依據丘陵及河谷區之 8 站地層柱狀圖，繪 製水文地質剖面十三至十五（深度約達 250 公尺左右），劃分出濁水溪沖積扇概 念分層，依深度分別為含水層一、阻水層一、含水層二、阻水層二、含水層三、 阻水層三及含水層四，濁水溪沖積扇模式共分為第一至七分層，如圖 4.2-1 所示。 本章將含水層一視為表層含水層，含水層二以下含水層視為深層含水層。並由 於含水層四之觀測井數量較少，僅有 11 口主要分布於沖積扇扇尾區域，不足以 提供建模所需資料，因此本計畫不將含水層四納入模擬範圍。

海 第一分層 第二分層 第三分層 第四分層 第五分層 第六分層 第七分層 模式概念分層 濁 水 溪 水 文 地 質 架 構 含水層1 含水層2 含水層3 含水層4 阻水層1 阻水層2 阻水層3 海岸 山麓 圖 4.2-1 濁水溪沖積扇概念分層 二、模式邊界條件與格網劃分 根據中央地質調查所(1999)對於濁水溪沖積扇地下水邊界分析研判，提出概 念性之邊界，如圖 4.2-2 所示。

圖 4.2-2 濁水溪沖積扇地下水邊界之概念圖

a. AB 段

位於車籠埔斷層上，斷層以東除河床表層厚約 10 公尺之河道沉積礫石層有 地下水伏流以外，均為透水及含水不佳之砂頁岩，屬於零流邊界(No flow boundary)。

b. BC 段 本邊界位於濁水溪沖積扇之南側邊緣上，含水層沉積物粒徑及厚度均顯著 變小，然而並未尖滅，因此無實體之阻隔為界，然而從地下水流網分布型態， 顯示 BC 段與地下水流線大致平行，於本身即可認定為一條流線，故亦屬於零流 邊界。海岸附近之洩降錐，其在邊界南北之形狀及大小若相當，則不影響零流 邊界之假設。 c. CD 段 此段為含水層尖滅封閉於阻水層之位置，屬於零流邊界。 d. DEF 段 本段與 BC 段相似，含水層可向北延伸而與台中盆地地下水系統相接，缺乏 實體之阻隔，惟與地下水流方向大致平行，故認定為零流邊界。EF 段位於和美 沖積扇頂上，河水可由河床入滲補注地下水。 e. FG 段 本段邊界大致與地下水等水位線平行，隨時間之不同水位有昇降變化。 f. GA 段 本段位於頭嵙山層之中，為八卦山台地之地下水分水線 (Groundwater divide)。 綜合上述環境特性，本計畫模式格網以 1 公里乘 1 公里之格網建構，南北 方向 85 列與東西方向 75 行，共五層格網。而邊界條件設定於 MODFLOW 模式 中，東部之八卦山地區有八卦山背斜，斗六丘陵有內林背斜及斷層經過，以此 為東邊邊界，且設為零流量之邊界條件。北部邊界為烏溪，南部邊界為北港溪， 各層皆為零流量之邊界，模擬區域如圖 4.2-3 所示。西邊含水層一以海岸線向外 延伸兩公里設為定水頭邊界，如圖 4.2-4 所示，其餘含水層延伸 5 公里後尖滅， 設為零流量之邊界條件，如圖 4.2-5 所示。圖中藍色網格為定水頭邊界，灰色網 格為零流量。圖 4.2-6 所示 AA’為剖面圖，水平向長度為垂向長度 30 倍。

圖 4.2-6 濁水溪沖積扇模式 AA’段剖面圖 4.2.2 模式輸入資料 模式之輸入資料包括儲水係數、透水係數、垂向透水係數、起始地下水水 位 、 補 注 量 及 抽 水 量 等 。 以 下 將 說 明 建 置 濁 水 溪 沖 積 扇 地 下 水 模 擬 模 式 (MODFLOW 模式)，其模式架構所需之地下水力學參數及源匯項輸入資料。各 輸入資料敘述如下： 一、淨補注量 在暫態模式中，水位會隨著不同時刻而改變，因此各時刻之淨補注量亦隨 時間改變。由於台灣地區現地抽水量資料缺乏官方記錄，抽水記錄殘缺，且地 下水天然補注亦屬未知，故無法直接輸入地下水數值模式，難以正向模擬，在 此以逆問題求推估補注量，本研究應用專家系統參數檢定系統，使得模擬水位 與觀測水位吻合，進行淺層淨補注量與深層抽水量之檢定。本研究設定之模擬 時刻長度為月，參數檢定之依據為各觀測井之月平均水位，此外，本研究採用 徐昇氏法進行參數分區，以降低參數維度。 在此受到模式的限制，檢定補注量時，設定分區的方式為每一徐昇氏網格 內需有一個地下水位觀測資料，且進行參數檢定時，每一徐昇氏網格分區內之 淨補注量皆需相同。 二、透水係數 觀測站網建置時，各觀測井必進行單井或複井抽水試驗，因此幾乎各觀測 井於不同含水層均擁有其導水係數（Transmissivity，T 值）現地試驗值，但因 A’ A 水平向每格為 1 公里 第 7 分層 第 5 分層 第 3 分層 第 1 分層 南 投 縣 台 灣 海 峽

MODFLOW 模式所需參數為透水係數(Hydraulic Conductivity，K 值)，因此本研 究以台灣地區地下水觀測網第一期計畫所得之導水係數，再除以含水層厚度（b） 求得透水係數，後續再如同淨補注量分區之方法，以徐昇氏網格為分區方式， 將各井點之參數值作為其對應徐昇氏分區之參數值，圖 4.2-7 至圖 4.2-9 分別為 模式第 1、3 及 5 分層之徐昇氏分區圖。而概念分層中之第 2 及 4 分層為阻水層， 其分布範圍位於濁水溪沖積扇中、下游區域，模式於沖積扇中、下游區域阻水 層之透水係數乃參考 Fundamentals of Ground Water (Franklin W. Schwartz and

Hubao Zhang, 2003)中建議之值域，在此設定為 4×10-3 (公尺/天)。此外，由於抽

水試驗所得結果為水平向之透水係數，故在此假設垂向透水係數為水平向之 1/50。

圖 4.2-8 含水層二徐昇氏分區 圖 4.2-9 含水層三徐昇氏分區 三、儲水係數 儲水係數輸入方式與淨補注量及透水係數相同，皆採用徐昇式網格分區， 將點位型式之資料給定至整個沖積扇。在數值方面，由於儲水係數必須進行複 井抽水試驗才能求得，因此儲水係數現地試驗數量，遠低於透水係數現地試驗 數量，僅有 10 口觀測井進行複井抽水試驗，如圖 4.2-10 所示，其值列於表 4.2-1(以 ※符號標註)。 其餘水井之儲水係數必須藉由現有的試驗值，搭配經驗判斷加以補齊。首 先，針對其餘水井以拘限與非拘限進行分類。其次，水井之儲水係數以同分類 水井之試驗平均值給定，拘限水井包含線西(1)、莿桐(1)、田洋(1)與海園(1)，其 試驗值分布於 4 10 91 . 1   _至4.48103_{之間，其平均值為} 3 10 8146 . 1   ，其餘水井即直 接給予 3 10 8146 . 1   ；非拘限水層也以相同方式推估，但由於東光(1)地處於古河道 流經之處，其地質架構較為特殊，比出水率高達 0.289，遠高於其餘 5 口現地試 驗值，為避免推估之比出水率偏高，故不予採用於計算平均值。因此比出水率

表 4.2-1 濁水溪沖積扇淺層觀測井列表(33 口井) 井名 X Y 儲水係數 S 比出水率 Sy 控制面積 (m2) 透水係數 (m/day) 含水層 種類 全興(1) 199514 2674258 - 0.1824 30620061.65 53.1936 0 線西(1) 195062 2669966 ※0.00448 - 72130899.48 20.491 1 洛津(1) 191220 2661365 - 0.1824 58184684.19 5.6966 0 文昌(1) 190120 2656250 0.0018146 - 130443980.9 17.856 1 花壇(1) 202695 2658255 0.0018146 - 119661405.5 49.356 1 香田(1) 186573.2 2641218 0.0018146 - 133513202.4 19.872 1 合興(1) 194000 2643600 - ※0.225 74315828.31 189.936 0 田尾(1) 201045 2643280 0.0018146 90079238.16 93.1376 1 竹塘(1) 190948 2639823 - 0.1824 45573810.79 154.368 0 溪洲(1) 199206.3 2639069 - ※0.216 32592225.09 50.22 0 柑園(1) 201662 2635696 - ※0.12 36836767.67 48.89 0 田中(1) 207088 2639188 - 0.1824 51189018.09 46.7856 0 二水(1) 210242 2634428 - ※0.199 23168436.14 65 0 豐榮(1) 178954 2632222 0.0018146 - 146221463.5 24.9984 1 西螺(1) 194891 2632723 - 0.1824 46974992.78 42.4832 0 海園(1) 165467 2624551 ※0.000191 - 73753709.91 18.5184 1 田洋(1) 178717 2624989 ※0.00188 - 68125351.65 43.6608 1 莿桐(1) 197896 2628388 ※0.0007074 - 32659495.23 56.73 1 芳草(1) 185400 2624155 - 0.1824 54256101.22 12.3912 0 九隆(1) 191168 2627781 0.0018146 - 70751253.87 77.544 1 六合(2) 204584 2629706 - 0.1824 32859021.1 133.2 0 烏塗(2) 208141 2629528 - ※0.152 23744105.22 61.92 0 箔子(1) 162598 2614898 0.0018146 - 101778284.1 2.9424 1 明德(1) 167485 2617020 - 0.1824 77871839.26 1.54448 0 東光(1) 174905 2616870 - ※0.289 145352763.2 45.4176 0 東和(1) 206145.3 2620189 - 0.1824 22888692.27 2.87 0 宏崙(1) 182680 2620675 0.0018146 - 48026563.67 33.3504 1 虎溪(1) 199331 2624542 0.0018146 - 47206609.82 16.7088 1 石榴(1) 206944 2624360 0.0018146 - 24621186.95 4.8 1 舊庄(1) 188020 2614870 0.0018146 - 79153871.55 20.088 1 溫厝(1) 199400 2617408 0.0018146 - 56160496.8 8.70336 1 古坑(1) 204980 2615932 - 0.1824 14367240.52 4.70592 0

崁腳(1) 202021 2612378 0.0018146 - 13911958.55 3.89952 1 ※表示複井抽水試驗實測值，其餘無實測資料之站以平均值給定

四、阻水層設定 本研究阻水層設定之依據乃藉由檢視水文地質鑽探柱狀圖、屏狀圖以 及同一觀測井之分層觀測水位差，綜合分析後，予以適當之調整。調整方 式說明如下所示： 首先檢視水文地質鑽探柱狀圖及屏狀圖後，若含水層內之某一範圍出 現厚泥層阻隔，本計畫則設定該範圍內之水力參數為阻水層等級之量值； 若含水層間為砂泥互層情形，則再進一步分析該處觀測井之不同分層間的 長期地下水位，若分層間之地下水位差異顯著，則代表該分層間有阻水效 果，故應視為有阻水層，若分層間之水位相近，則此處將不視為有阻水層。 五、起始地下水水位 對於暫態模式而言，初始地下水位代表模式之初始條件，不同初始條 件影響模擬結果極巨，因此應選取適當之初始水位以得到較佳之模擬結果。 本研究建置之地下水系統模擬期距為 2012 年，為月模擬，模擬起始時刻為 西元 2012 年 1 月，因此以西元 2011 年 12 月之平均水位作為起始水位。

4.3 現地重力量測結果

由於重力值會受到高程變化之影響，因此量測時需在相同位置進行量 測，但由於該區域有嚴重之地層下陷問題，故所得之重力值宜藉由觀測下 陷量進行校正，為此本研究於溪州國小(彰化縣溪州鄉)、土庫國中(雲林縣 土庫鎮)及客厝國小(雲林縣元長鄉)等三站，進行重力調查。此外，重力站 亦搭配自行架設之 GPS，藉此觀測各站受地層下陷之狀態，以便後續校正 地層下陷施測之影響，地點座標詳見於表 4.3-1，。 重力量測使用絕對重力儀(FG5)進行現地監測(如圖 4.3-1)，並配合連續 式 GPS 監測站觀測地表高程之變動(如圖 4.3-2)，各重力量測站分佈及所在 之鄉鎮如圖 4.3-3。其中，絕對重力儀器量測值移除及改正海潮、固體潮、 極移、大氣壓力及高程變化等影響，如式 4.3-1 表示，在此稱之為重力殘 差。重力量測值之改正除了高程變化外，其餘皆為環境影響改正。有關各 環境改正內容詳見附錄 A。 重力殘差 觀測量 海潮負載 固體潮 極移改正 大氣效應 高程影響………(4.3-1) 為配合台灣地區豐水期及枯水期地下水位之變化，重力量測時間於 2012 年之 5 月、8 月及 10 月，共計有 8 次重力量測記錄，各監測點分別有 2 至 3 次的量測結果。將取得之 GPS 高程資料以周為單位，但由於資料之 不完整性，再加上避免 GPS 量測時所造成的誤差，故取 7 周移動平均，使 GPS 量測高程起伏較為平緩，結果如圖 4.3-4 至 4.3-6，觀測結果節錄自 101 年度重力基準維護及測量整合服務工作之計畫報告。 Hofmann-Wellenhof 與 Moritz(2006)指出，在布格平板之假設下，如重 力站高程上升 1(cm)，重力量測值將減少 1.9( -Gal)。因此如 GPS 觀測值觀

測到 1(cm)之下陷量，重力值應減少 1.9( -Gal)作為校正。本研究將五月高 程設定為一基準點，結果之重力殘差值如表 4.3-2 至 4.3-4 所示。

表 4.3-1 絕對重力儀施測地點 重力測站 緯度 經度 行政區 溪州國小 SJES 23.85380 120.49820 彰化縣溪州鄉 土庫國中 TKJH 23.68799 120.38982 雲林縣土庫鎮 客厝國小 KTES 23.62661 120.33428 雲林縣元長鄉 表 4.3-2 溪州國小量測結果 量測日期 標準偏差 ( -Gal) GPS 高程 (m) 重力殘差( -Gal) 2012/5/23~24 4.7 56.431 978834152 2012/8/16~17 4.8 56.424 978834175 2012/10/24~25 4.8 56.416 978834163 表 4.3-3 土庫國中量測結果 量測日期 標準偏差 ( -Gal) GPS 高程 (m) 重力殘差( -Gal) 2012/5/22~23 4.8 40.966 978843350 2012/8/17 4.8 40.964 978843372 2012/10/25 4.8 40.947 978843359 表 4.3-4 客厝國小量測結果 量測日期 標準偏差 ( -Gal) GPS 高程 (m) 重力殘差( -Gal) 2012/5/21~22 5.1 36.147 978850626 2012/10/30~31 4.8 36.136 978850635 資料來源：101 年度重力基準維護及測量整合服務工作－期末報告

圖 4.3-4 GPS 連續監測地表高程變化曲線－客厝國小 圖 4.3-5 GPS 連續監測地表高程變化曲線－土庫國中 36.13 36.135 36.14 36.145 36.15 36.155 36.16 36.165 36.17 高程 (m) 日期 客厝國小高程變化 移動平均線 40.94 40.945 40.95 40.955 40.96 40.965 40.97 40.975 40.98 40.985 40.99 高程 (m) 日期 土庫國中高程變化 移動平均線

圖 4.3-6 GPS 連續監測地表高程變化曲線－溪州國小 56.41 56.415 56.42 56.425 56.43 56.435 56.44 56.445 高程 (m) 日期 溪州國小高程變化 移動平均線

4.4 應用重力變化於濁水溪沖積扇比流出率之檢定

由於地下水數值模式之網格為一公里見方之矩形，厚度則依據含水層 厚度而變，約數公尺至數十公尺不等，如直接進行數值積分，過粗之積分 網格則可能產生數值積分誤差。一般而言，數值積分精度隨積分網格之尺 寸越小而相對提升，但是積分網格越小則計算量越高。 基於牛頓重力理論，物體質量越接近量測站，其重力貢獻量越大，因 此較佳之積分網格切割策略應依據重力量測站之位置，逐步往內加密縮 小。 為滿足上述原則，本研究以兩階段混合網格進行數值積分，以重力量 測站所在之 1 公里乘 1 公里之第一層矩形網格進行加密，兩個維度均進一 步加密切割為 20 段，亦即水平方向切割為 50 公尺乘 50 公尺之加密網格， 垂直方向則不予切細。加密區域外，則沿用原始地下水數值模式之網格， 以節省計算量。 在加密區中，單一地下水模式網格會因為加密積分的緣故，其計算量 提升為非加密區之 400 倍，在此採取此兩階段混合網格，可在小幅增加計 算量之情況下，有效提升積分精度。 4.4.1 重力積分半徑影響分析 為利用重力觀測值，以檢定系統儲水係數或比流出量，各重力觀測站 係受到觀測站附近一定半徑內水量之影響，因此僅需調整該半徑內之參數 值，特此應用半徑影響分析決定影響半徑。 由於本系統是以數值方式進行重力積分，因此藉由逐步增加半徑，以 確定積分半徑之影響。選定 8 個地下水觀測水位站，各點位散布於扇頂、 扇央及扇尾，涵蓋侷限含水層與非侷限含水層，以此 8 站所在位置分別進 行重力積分半徑之計算。

接著，逐步增加與調整積分半徑，由 1000 公尺、3000 公尺、4000 公 尺增加至 20000 公尺及無限遠等，得到不同影響範圍內之模擬重力值。最 後，再利用無限遠之重力模擬值( 與各不同範圍內之模擬重力值計算出 重力貢獻百分比( )，如式 4.4-1 所示。 ...(4.4-1) 其中， 代表重力貢獻百分比， 代表無限遠之重力模擬值， 代 表 l 範圍內之重力模擬值。由圖可知，無論各區域水文地質參數之差異， 重力貢獻百分比( )曲線之趨勢相當一致，且多數重力站於半徑為 5000 公 尺時，其重力貢獻百分比多高達 97%以上，顯現該站重力值多受 5 公里內 之參數影響，因此本研究以 5 公里為影響半徑，後續檢定上，則以重力站 為中心，以 10 公里見方之網格為檢定分區。 圖 4.4-1 影響半徑分析圖 96 97 98 99 100 0 5000 10000 15000 20000 25000 30000 重力貢獻百 分比 (%) 距離(m) 東光(1) S=0.289 合興(1) S=0.225 溪洲(1) S=0.216 明德(1) S=0.1824 線西(1) S=0.00448 田洋(1) S=0.00188 舊庄(1) S=0.001815 莿桐(1) S=0.000707

4.4.2 比流出率檢定 前述本研究已應用極少量之現地試驗，其餘無試驗值者，則以平均值 給定，如表 4.2-1 所示，初步建置、檢定與模擬濁水溪沖積扇地下水流水 位與水量變化。由模擬之地下水量變化，並利用重力積分程式，計算各重 力站地下水之重力貢獻值。 以模擬重力貢獻值變化幅度與量測值變化幅度之比值，作為參數調整 之依據，儲水係數與比出水量之參數調整分區則以重力站為中心，則依據 前述積分半徑影響分析之結果，於水平方向建立 10 公里見方之區域，垂直 方向則侷限於第一層含水層作為該參數調整分區。 政府於民國 100 年底完成，土庫國中與客厝國小地層下陷相關監測， 除了地陷井與 GPS 觀測站外，亦於場址建立地下水觀測站，於民國 100 年 完工，尚未開始進行地下水之觀測。本研究以重力觀測作為儲水係數或比 流出量之參數檢定，需再調整淨補注量使得地下水模擬水位與觀測水位吻 合，因此對於參數調整過的區域，需以觀測水位作為淨補注量之檢定目標。 為驗證本研究方法之可行性，在此應用克力金(Kriging)內插演算法推 估土庫國中與客厝國小之水位，以內插水位作為該站觀測水位，作為虛擬 觀測站，後續如兩觀測站已進行觀測，其可直接採用觀測水位。溪洲國小 站雖無建置地下水觀測站，但地下水觀測網已於溪州(1)建置觀測站，兩者 距離約距 1 公里遠，由於原有地下水觀測站距離重力站不遠，故無需以前 述方式推估重力站位置之觀測水位。在加入土庫國中與客厝國小兩站後， 其抽水分區則以新加入虛擬站井，建立徐昇氏網格分區，其分布如圖 4.4-2 所示。 圖 4.4-3 與 4.4-4 分別為客厝國小與土庫國中之內插水位，以及鄰近觀 測站之水位變化圖，以客厝國小為例，其鄰近站井為東光(1)、宏崙(1)與舊

庄(1)，由於濁水溪沖積扇以東方為上游，西方為下游，因此三者依據東西 方之分布逐步遞減，客厝國小與宏崙(1)之上下游區位較為相近，故內插水 位亦與宏崙(1)之水位相類似。以搭配前述徐昇氏分區與內插水位，重新進 行濁水溪沖積扇淨補注量之自動檢定，本研究應用專家系統參數檢定系統 進行自動檢定。

圖 4.4-3 客厝國小內插水位與鄰近三站觀測水位之變化 圖 4.4-4 土庫國中內插水位與鄰近三站觀測水位之變化 藉由專家系統參數檢定系統完成自動檢定後，由於參數本身並無調整， 且增加之虛擬站井亦以內插方式取得，故系統增加虛擬站井與否，對於系 統淨補注量、系統水位與水量並無大幅差異。 將系統水量配合重力積分公式，可以推算地下水水量變化對溪洲國小、 土庫國中與客厝國小等三站之重力貢獻值變化，比較實測之重力殘差變化 -2.0 0.0 2.0 4.0 6.0 8.0 10.0 12.0 14.0 16.0 1 2 3 4 5 6 7 8 9 10 11 12 水位高 (m) 時間(月) 東光(1) 宏崙(1) 舊庄(1) 客厝國小 6.0 7.0 8.0 9.0 10.0 11.0 12.0 13.0 14.0 15.0 16.0 1 2 3 4 5 6 7 8 9 10 11 12 水位高 (m) 時間(月) 芳草(1) 宏崙(1) 舊庄(1) 土庫國中

對模擬重力變化之比值，可做為儲水係數或比出水量檢定之依據。 由少量現地試驗所得之比出水率及儲水係數，透過 MODFLOW 模擬 與檢定完成地下水系統，並配合重力積分公式進行演算後，由 33 個觀測水 位與 2 站虛擬水位站，分別展示分佈於濁水溪沖積扇之比出水率、水位振 幅及重力振幅在空間上的變化，如圖 4.4-5 至 4.4-7 所示。首先，由圖 4.4-5 可看出，扇頂皆為靠近山麓的地區，其地質多為礫石所組成，顆粒較扇央 或扇尾大許多，而扇央區域，在合興及竹塘也有較大的儲水係數。此外， 扇頂區域分別由田中、二水、烏塗、六合及古坑等站所組成，依照中央地 質調查所岩心柱狀圖搭配地質剖面圖判斷可知，該區域為非侷限含水層， 故水位的變化相較於扇央及扇尾來得大；而扇尾區域之東光站位於古河道 上，該區域有局部且特殊的地質架構，因此雖接近扇尾但其水位振幅也較 大，如圖 4.4-6。由圖 4.4-7 可看出，模擬重力變化與模擬水位變化及儲水 係數之間具有正相關。