漳浦縣舊鎮(zhèn)梅宅村海水入＊侵的動力系統(tǒng)模型研究

蔡建堤

(福建省水產研究所,福建 廈門 361012)

漳浦縣舊鎮(zhèn)梅宅村位于福建省南部沿海。降水量一般年份全年降水日數(shù)100多天,1-3月份常有小雨,4-6月份多雷陣雨,7-9月份常有臺風暴雨,10-12月份少雨。鹿溪是舊鎮(zhèn)主要河流,鹿溪在入海處設有水閘,用于排洪水和含蓄淡水,斷開和海水的直接聯(lián)系。

漳浦縣舊鎮(zhèn)梅宅村一帶有上千畝的池塘圍墾養(yǎng)殖,當?shù)仞B(yǎng)殖戶大量抽取地下水進行水產養(yǎng)殖。長期不合理和過量開采地下水,造成海水入侵淡水含水層,養(yǎng)殖區(qū)附近的地下水海水入侵十分嚴重,且入侵有進一步加深的趨勢。因此我們加強對地下水資源的監(jiān)管和建立海水入侵的有效預警機制尤為重要。

海水入侵是世界上海岸含水層中十分普遍和廣泛分布的污染問題,認識、預測、防治以及含水層管理最好的途徑便是進行監(jiān)測和計算機的數(shù)值模擬[1-2,4,8,13]。在海水入侵模擬的數(shù)學模型與數(shù)值方法方面做了大量的研究工作[5,7,9,14-15],數(shù)學模型主要有突變界面模型和過渡帶模型兩種;主要的數(shù)值方法,海水入侵水質模型求解方法以有限元方法居多,為了解決求解對流彌散方程所引起的數(shù)值彌散及過量現(xiàn)象,所采取的方法有迎風(上游加權)有限差分法或有限元法、交替方向隱式差分或有限元法、特征法和變形網格法、人工彌散加權法及歐拉拉格朗日混合解法(即特征有限元法,簡稱ELM法)等。但是,這些方法的推廣還存在一定困難,特別是在三維情況下,而且在程序中對邊界的處理也是十分復雜、費時的。

動力系統(tǒng)模型[16]廣泛用于農業(yè)、金融、醫(yī)療等領域,獲得很大的成功。本文首次推導出海水入侵預警的動力系統(tǒng)模型,獨創(chuàng)模型參數(shù)的求解方法,并對模型參數(shù)進行靈敏性檢驗,建立模型的仿真系統(tǒng)。本模型不需要詳細了解：海水入侵過程中交換Na+、Ca2+的運移行為[6,12],咸淡水界面[3,10-11]和海水的潮水變化的影響等復雜問題。通過歷年的海水入侵監(jiān)測數(shù)據(jù)：地下水開采量和監(jiān)測站點礦化度等信息,確定模型參數(shù),建立模型,預測海水入侵的發(fā)展趨勢,提出科學的防治措施。目前,國內有關動力系統(tǒng)模型研究海水入侵的方法未見有報道。

1 資料和方法

1.1 研究區(qū)位

受國家海洋局和福建海洋與漁業(yè)局的委托,對福建省漳州地區(qū)進行海水入侵監(jiān)測。在漳浦縣舊鎮(zhèn)梅宅村設1個斷面,斷面布設站位3個監(jiān)測站位(海水入侵區(qū),過渡帶,人為入侵區(qū)),站號分別是1號、2號、3號(圖1)。

水樣品中礦化度測定,采用重量法(SL79-1994),樣品的檢測結果,見表1。

圖1 海水入侵監(jiān)測站位示意圖Fig.1 Station locations for seaw ater intrusionmonitoring

表1 海水入侵監(jiān)測數(shù)據(jù)Table 1 Seaw ater intrusion data from themonitoring

1.2 方法

1.2.1 推導模型的數(shù)學表達式

漳浦縣舊鎮(zhèn)梅宅村海水入侵是海灣層狀沉積體孔隙式層狀入侵,此種入侵取決于咸、淡水之間的平衡狀態(tài)。一般來講,阻止海水入侵必須形成一定的淡水水頭壓力,為此在原有平衡基礎上需要有源源不斷的淡水供應。咸、淡水之間平衡的破壞主要表現(xiàn)為地下淡水資源量的明顯減少,即地下水開采量明顯大于地下水在相同時間段內的儲積增量,造成咸、淡水相互作用帶上的地下淡水水位下降、水頭壓力不足,不足抵制海水水頭壓力,導致海水入侵。對基巖海岸海水入侵來講,也應具備類似上述的地下水動力學特征。因此,淡水水頭壓力引起監(jiān)測點淡水水位的變化(?d/?t)應該和地下水的水頭水位(d)和海水水頭的水位(h)有關。淡水水頭水位(d)造成淡水水頭壓力,增大監(jiān)測點淡水的供應量,使得監(jiān)測點淡水水位的變化(?d/?t)變大,從而抵制海水入侵。由于淡水含水層有一定的水容量,因此,淡水和海水存在競爭關系。當?shù)叵滤^水位升高,造成監(jiān)測點淡水水位的變化(?d/?t)變大,使得海水水位(?h/?t)變小。同時,監(jiān)測點淡水水位的變化(?d/?t)和監(jiān)測地點距離淡水水頭的距離(x2)有關,離淡水水頭越近,淡水水位的變化(?d/?t)就越大。

因此,監(jiān)測點淡水水位的變化和監(jiān)測點離淡水水頭地點(x2)的平方成反比和淡水水頭水位d的平方成反比。具有這樣性質的海水入侵的動力系統(tǒng)模型為：

式中,α2為淡水擴散率,與淡水密度和水文地質有關;β為淡水含水層資源限制強度度量,與水文地質有關;x2為離淡水水頭的距離;k2為抽取監(jiān)測點的地下水中淡水水位。

同樣,監(jiān)測點海水水位變化(?h/?t)的動力系統(tǒng)方程為：

式中,α1為海水擴散率,與海水密度和水文地質有關;β為淡水含水層資源限制強度度量,與水文地質有關;x1為監(jiān)測點離海水水頭的距離;k1為抽取監(jiān)測點的地下水中海水水位。

海水入侵的動力系統(tǒng)和淡水含水層的實際水位、水文地質、海水的潮汐變化、降雨量、開采地下水的量、地下水賦存條件等很多因素有關,而有些因素我們無法準確的得知。如果建立模型都考慮這些因素,那么將使得模型建立或者求解變得很困難,甚至無法建立或者求解,這樣的模型毫無意義。因此,我們應該對模型進行簡化。由于海水水量相比淡水水量大的多,監(jiān)測點抽取地下水用于養(yǎng)殖用,對淡水頭水位有較大的影響,對海水頭水位影響忽略不計,同時從時間大的尺度(年)上,不考慮海平面變化和潮汐作用影響。另外,式(1)和式(2)中,當?shù)氐乃牡刭|情況反映在 α1,α2和β參數(shù),因此,不需要知道實際水文地質。這樣簡化的動力系統(tǒng)模型,考慮了海水入侵最基本要素的淡水和海水水頭水位、水文地質、地下水開采量。這樣的模型簡單,且有實際的意義,能預警海水入侵。

1.2.2 模型參數(shù)的求解方法

隨機取一批淡水水頭和海水水頭的水樣,測其礦化度,淡水水頭水樣礦化度平均值為m,海水水頭水樣平均的礦化度為n,根據(jù)公式：監(jiān)測點的礦化度=n×c+m×(1-c),求解出海水比例c和淡水的比例(1-c)。通過監(jiān)測斷面各個監(jiān)測點的觀測和測量,取得監(jiān)測點的礦化度,采水量k1+k2,x1,x2和監(jiān)測點的水位變化(?d/?t+?h/?t)。這樣,只要利用兩個監(jiān)測點的數(shù)據(jù),就可以計算出 ?d/?t和?h/?t,再利用監(jiān)測點的淡海水的比例,計算出該監(jiān)測點的水位大小。根據(jù)監(jiān)測點的水位大小和當?shù)貙嶋H海水開始入侵的時間,估算出淡水水頭水位(d)。根據(jù)Fick化學物體擴散的模型,擴散率和擴散物質的密度有關,假設淡水的擴散率α2=1,則海水的擴散率α1=1.03。這樣就可以計算出所有參數(shù)h,β,k1和 k2。

2 分 析

2.1 漳浦縣舊鎮(zhèn)梅宅村動力系統(tǒng)模型的建立

測量海水礦化度的平均值為34.369 5 g/L,淡水礦化度的平均值為0.179 5 g/L。根據(jù)公式：礦化度=34.369 5×c+0.179 5×(1-c),算出各個監(jiān)測站點海水比例c和淡水的比例(1-c),以及淡水和海水的比(1-c)/c(表2)。

表2 淡水海水的比例Table 2 The proportion of freshwater to seaw ater

1號監(jiān)測站位 x1=2.50 km,x2=2.65 km(圖1)。2008年4月-2009年4月,抽水引起的1號監(jiān)測站位井水水位年變化量為?d/?t+?h/?t=-0.044 6單位,抽水量年變化量為k2+k1=0.682 0單位,其中1單位=5 m。由式(1)+式(2),得

同樣,3號監(jiān)測站位 x1=5.15 km,x2=0.5 km(圖1);2008年4月-2009年4月,抽水引起的3號監(jiān)測站位井水水位年變化量為?d/?t+?h/?h=-0.100 0單位,抽水量年變化量為k2+k1=4.172 5單位。

由式(1),(2),(3),(4)得：

d11表示1號監(jiān)測站位2008年4月淡水水位,h11表示1號監(jiān)測站位2008年4月海水水位,由表3得：

由式(7),(8),(9),(10)得：h11=1單位,d11=3.044 6單位。即2008年監(jiān)測站位的水位為4.044 6單位。根據(jù)Fick化學物體擴散的模型,擴散率和擴散物質的密度有關。假設α2=1,則 α1=1.03,據(jù)此估算淡水源的水位 d=5單位。由(5)和式(6)得h=5.121 8單位,β=0.935 6。由式(1),(2),(5),(6),(7),(8)得,k1=0.191 4單位,k2=0.490 6單位。

因此,1號監(jiān)測站位的動力系統(tǒng)模型方程為：

2.2 海水入侵的趨勢分析

當?h/?t=0的時候,海水停止入侵,假設海水源水位保持不變。如果我們停止開采地下水資源,k1=0,k2=0。由式(11)得出d=5.638 6單位,由式(12)得出 ?d/?t=0.189 1。

初始值h=1單位;d=3.149單位,保持淡水水位為5.638 6單位,可以阻值海水入侵,停止開始地下水,經過t=(5.638 6-1-3.044 6)/0.189 1=8.429 4 a后,1號監(jiān)測站位的井水礦化度淡化為6.243 09 g/L。

當我們停止開采,并讓淡水源水位保持在5.638 6單位,經過大約8 a海水礦化度減輕為6.243 0 g/L。因此一旦發(fā)生海水入侵,要比較長時間減輕海水入侵,且恢復比較緩慢。

假設我們停止開采地下水,并讓淡水水位上升并保持為6.2單位,則?h/?t=-0.621 8。

初始值h=1單位;d=3.044 6單位,最高水位 6.2單位,經過 t=1/0.621 8=1.608 1 a,經過大約2 a海水完全淡化。表明：增加淡水的資源可以加速抵制海水的入侵。

圖2 動力系統(tǒng)模型的向量圖Fig.2 Vector diagram from the dynam icmodel

動力系統(tǒng)模型的向量圖(圖2),可以清晰分析漳浦縣舊鎮(zhèn)梅宅村1號監(jiān)測站位淡水和海水的變化關系。當?shù)此籨保持不變,且h＞1.100 1×d時,海水水源h越大,向量越大,向量向負的方向,說明第2年1號監(jiān)測點淡水水位下降越明顯,海水的水位則上升明顯。再考慮,海水水源h保持不變,且d＞1.068 8×h時,淡水水源水位d變大,向量越大,向量向正的方向,說明第2年1號監(jiān)測點淡水水位出現(xiàn)明顯上升,海水水位下降。向量場在h=d附近,向量大小比較小,說明如果保持淡水水源水位和海水水源水位一樣大時,第2年1號監(jiān)測點的海水入侵變化不明顯,動力系統(tǒng)到達基本的平衡態(tài)。如果增加淡水源水位,使得h＞1.068 8×h,則監(jiān)測點的第2年的礦化度變小,減緩海水入侵。這為我們科學合理的利用地下水資源而幾乎不破壞地下動力系統(tǒng)的平衡,提供理論基礎。

2.3 模型的預警

采用兩種預警方式,靜態(tài)預警方式：根據(jù)當?shù)貙嶋H估算第2年的淡水水頭水位和地下水用水量等參數(shù)(表3)。動態(tài)預警方式：根據(jù)實際的降水量和地下水用水量,動態(tài)調整淡水水頭水位和開采量大小等參數(shù)(表4)。

表3 靜態(tài)預警模型方式Table 3 Static prediction model

表4 動態(tài)預警模型方式Table 4 Dynam ic predictionmodel

兩種預警方式主要誤差原因,2010年3-4月漳浦出現(xiàn)大量的降水,增加了淡水的水位,同時用于養(yǎng)殖的地下水開采量也變小,海水入侵減慢。動態(tài)預警模式估算的礦化度和實際監(jiān)測結果比較接近,誤差在2%以內,預警精確度比較高。而靜態(tài)模型主要優(yōu)點在于能提前對海水入侵趨勢做出初步的預測,為水資源的規(guī)劃和管理提供參考依據(jù),缺點是預測進度比較差,誤差在15%左右。

3 討 論

1)動力系統(tǒng)在建模過程中對一些問題做出假設,我們很少保證這些假設都是完全正確的。因此我們需要考慮模型對每個假設的敏感程度,這就是靈敏性分析,是數(shù)學建模的一個重要面。我們不能對模型中的每個參數(shù)都計算靈敏性系數(shù),只選擇那些有較大不確定性的參數(shù)進行靈敏性分析。本模型的推導過程,估算淡水水頭水位d=5,并假設α2=1,對于這些假設,我們需要對靈敏性分析,確定對其它參數(shù)的敏感程度。

h對d的靈敏性分析,靈敏度分析公式為：

由式(13),當我們在 d=5,h=5.121 8時,由(5)(6)得出：

由式(14)說明：若d增加1.000 0%,將h只增加0.796 6%。說明d的不確定性變化,會導致h的稍微小的變化,但是這樣變化不是很明顯。

β對d的靈敏性分析,靈敏度分析公式為：

由式(15),當我們在 d=5時,β=0.935 6由(5)(6)得出：

由式(16)說明：若d增加1%,β將大約減少1%。說明我們原始假設的d不確定度,對β成反比,而對h影響成正比。

同樣,模型的推導假設α2=1,我們確定與α2相關的主要靈敏性的分析。

β對α2的靈敏性分析,靈敏度分析公式為：

h對α2的靈敏性分析,靈敏度分析公式為：

由式(17),(18)和(19)得：如果 α2增加1%,則 β和h增加均1.043 4%,而 d減少0.960 2%。說明我們原始假設的α2不確定度,對d成反比,而對 β和h影響成正比,且影響的比例一致。

2)應用MATLAB數(shù)學軟件的SIMULINK對模型進行計算機仿真,通過建模、仿真和分析,可以實時修改參數(shù),減少繁雜的數(shù)學運算和推導,對模型的推廣有著重要的作用。如：對式(2),建立仿真模型,組建各個模塊(圖3)。在2個常值信號輸入模塊d,輸入淡水水位d的值;在常值信號輸入模塊h,輸入海水水頭水位h值;在2個常值信號輸入模塊x2,輸入離淡水水頭的距離x2的值;在常值信號輸入模塊k2,輸入抽取監(jiān)測點的地下水中淡水水位k2;在常數(shù)增益模塊a2,輸入淡水擴散率α2值;在常數(shù)增益模塊b,輸入淡水含水層資源限制強度度量β值;動態(tài)仿真結果就可以在顯示輸入模塊Disp lay1上顯示?d/?t的值。圖3中的參數(shù)是按式(12)中的參數(shù)輸入的,容易得出 ?d/?t?-0.324 6。同樣,式(1)的仿真模型和圖3大致相同。

d對α2的靈敏性分析,靈敏度分析公式為：

圖3 動力系統(tǒng)模型仿真Fig.3 Flow chart for the dynam icmodel

3)模型預警需要通過監(jiān)測點的水位、礦化度、采水量,預警海水入侵。除了礦化度數(shù)據(jù)可以精確測量外,水位數(shù)據(jù)是通過當?shù)氐慕邓拷Y合井深水位估算、采水量根據(jù)養(yǎng)殖戶的用水量估算,精確度難以保證。為了獲得更為準確的數(shù)據(jù),需要改進監(jiān)測設備,如：采用YSI等自動測定儀器,連續(xù)跟蹤監(jiān)測井水的鹽度、溫度、水深,可以得到更準確和詳細的井水開采量和礦化度,從而完善動力系統(tǒng)模型,對精確預警海水入侵有著重要的意義。此外,模型參數(shù)在一定程度上反映當?shù)氐膶嶋H地質特征,因此每個地方參數(shù)是不同的,推廣應用在其它海水入侵地方需要重新計算參數(shù)。動力系統(tǒng)模型是一種新的研究海水入侵的方法,還需要在預警實踐中不斷地完善。

4 結 論

海水入侵預警的動力系統(tǒng)模型特點：1)模型建立、求解容易;2)模型參數(shù)簡單：地下水開采量、監(jiān)測站點礦化度、降水量等信息,不需要了解復雜的水文地質條件;3)采用動態(tài)預警方式預警海水入侵,誤差在2%以內,精度高。

漳浦縣舊鎮(zhèn)梅宅村海水入侵的特點：1)要比較長時間減輕目前的海水入侵,大約需要8 a時間,且恢復比較緩慢,礦化度=6.243 0 g/L;2)采用大量增加淡水資源的措施可以加速抵制海水的入侵,如：讓淡水水位上升并保持為6.2個單位,大約只需要2 a,監(jiān)測點完全淡化。3)合理保持淡水水位和海水水位一致,則海水入侵不會加劇。

通過以上綜合分析,模型能有效預警漳浦縣舊鎮(zhèn)梅宅村海水入侵,為我們合理規(guī)劃和管理水資源提供理論依據(jù)。

致謝：參加本項目的人員還有姜雙城、陳財珍、洪云、鐘碩良,謹致謝忱。

[1] 丁玲,李碧英,張樹深.海岸帶海水入侵的研究進展[J].海洋通報,2004,(2)：82-87.

[2] 王潘平,李天科,王兵,等.萊州灣海水入侵原因分析與防治措施[J].山東農業(yè)大學學報：自然科學版 ,2009,(1)：93-97.

[3] 艾康洪,陳崇希.漫尾島咸淡水界面運移剖面二維流水質模型模擬研究[J].勘察科學技術,1994,(6)：3-9.

[4] 安永會,張福存,賈惠穎,等.萊州灣濱海平原咸水入侵分析與趨勢預測—以山東省廣饒縣為例[J].環(huán)境科學與技術,2009,(5)：79-82.

[5] 成建梅,黃丹紅,胡進武.海水入侵模擬理論與方法研究進展[J].水資源保護,2004,(2)：3-10.

[6] 吳吉春,薛禹群.海水入侵過程中水-巖的陽離子交換[J].水文地質工程地質,1996,23(3)：18-19.

[7] 張乃興,鄭新奇,李新運.萊州灣東南沿岸海水入侵發(fā)展規(guī)律與趨勢預測[J].山東師范大學學報：自然科學版,1996,11(4)：72-77.

[8] 李國敏.海水入侵研究現(xiàn)狀與展望[J].地學前緣,1996,(2)：161-168.

[9] 李新運.萊州灣東南沿岸海水入侵相關分析和趨勢預測[J].中國地質災害與防治學報,1994,5(4)：20-23.

[10] 姚秀菊,王洪德.黃河三角洲地區(qū)地下淡水(微咸水)的形成與演化[J].地球學報,2002,23(4)：375-378.

[11] 姜效典,王碩儒.海水入侵地區(qū)咸水與淡水分界面計算-B樣條函數(shù)方法的應用[J].青島海洋大學學報：自然科學版,1995,25(2)：213-220.

[12] 郭篤發(fā).萊州灣東南岸海水入侵區(qū)地下水中若干離子的主成分分析[J].海洋科學,2004,(9)：6-9.

[13] 黃奕普.海水入侵問題研究綜述[J].水文,2003,(3)：10-15.

[14] 蔡祖煌,馬鳳山.海水入侵的基本理論及其在入侵發(fā)展預測中的應用[J].中國地質災害與防治學報,1996,7(3)：1-9

[15] 潘世兵,張福存,安永會.咸水入侵趨勢預測的非確定模型研究[J].勘察科學技術,1999,17(6)：5-8.

[16] M ARK M.Meerschaert.MathematicalModeling[M].北京：機械工業(yè)出版社,2005.