任紅蕾,陶月贊,林 飛,韋 婷
(合肥工業(yè)大學(xué) 土木與水利工程學(xué)院,安徽 合肥 230009)
在地下水滲流力學(xué)中,受第一類(lèi)邊界(如河渠水位)條件控制的一維潛水非穩(wěn)定滲流模型,是最基本的經(jīng)典模型之一[1-4],是灌排溝渠工程設(shè)計(jì)[5]、水文過(guò)程分析與水資源量評(píng)價(jià)[6-9]、排水工程附近地下水動(dòng)態(tài)過(guò)程分析[10-11]等研究的重要理論工具。近年來(lái),一維潛水非穩(wěn)定滲流模型進(jìn)一步拓展到其它研究中,如水岸生態(tài)工程設(shè)計(jì)中的地表水與地下水之間的水量交換過(guò)程評(píng)價(jià)[12-13]、濱海含水層保護(hù)工程設(shè)計(jì)中的地下水水位對(duì)海水位及潮汐的響應(yīng)分析[14]。
關(guān)于這類(lèi)經(jīng)典模型解析解問(wèn)題的研究,多關(guān)注含水層接受補(bǔ)給影響的河渠類(lèi)水位邊界[15-18],對(duì)僅能接受含水層排泄的排水溝渠類(lèi)邊界多側(cè)重?cái)?shù)值解[11,19],而少見(jiàn)有解析解研究文獻(xiàn)。在河渠類(lèi)水位邊界附近非穩(wěn)定滲流模型解析解問(wèn)題研究中,多未考慮降水入滲補(bǔ)給、灌溉回歸和潛水蒸發(fā)等潛水含水層垂向水量補(bǔ)排作用[1-4,15-16],并在假設(shè)河渠邊界水位瞬時(shí)上升ΔH且保持不變[1-4]或在對(duì)變化水位適當(dāng)處理[15-16]的條件下,討論模型的解及其應(yīng)用。在同時(shí)考慮邊界與垂向水量補(bǔ)排作用的文獻(xiàn)中,文獻(xiàn)[17]對(duì)于時(shí)變的潛水含水層垂向水量交換強(qiáng)度ε(t),依據(jù)ε(t)為常數(shù)時(shí)模型的解,對(duì)ε(t)進(jìn)行離散,采用疊加方法獲得時(shí)變問(wèn)題的解;文獻(xiàn)[18]討論了依據(jù)解求算模型參數(shù)的方法。
潛水位動(dòng)態(tài)的計(jì)算與預(yù)測(cè),是農(nóng)田灌排水渠系設(shè)計(jì)等實(shí)際問(wèn)題研究的理論基礎(chǔ)。潛水具有自由水面,通過(guò)包氣帶可與外環(huán)境直接相聯(lián)通,潛水蒸發(fā)、降水入滲補(bǔ)給、灌溉回歸等潛水含水層的垂向水量補(bǔ)排作用在潛水水位動(dòng)態(tài)研究問(wèn)題中往往難以回避。
排水期間溝渠水位特征有別于河渠,為保證排水效率,排水期間溝渠水位一般不會(huì)有較大幅度上升,如淮北平原沿淮地帶,潛水水位埋深較淺,普遍設(shè)置排水大溝調(diào)控地下水水位,排水大溝設(shè)置有滾水壩[19]以控制排水水位。在潛水蒸發(fā)和潛水向溝排泄的共同作用下,農(nóng)田區(qū)的潛水水位緩慢下降,但溝中水位在滾水壩控制下基本維持不變。這一過(guò)程,可能持續(xù)數(shù)天甚至數(shù)十天,此時(shí)溝渠水位顯然是穩(wěn)定水頭(水位高度為滾水壩壩頂高度)邊界。
數(shù)學(xué)模型的解析解,是問(wèn)題機(jī)理研究、變量間關(guān)系定量分析的重要工具。就第一類(lèi)邊界附近的一維潛水滲流模型而言,隨著源匯項(xiàng)和邊界條件函數(shù)形式的變化,模型求解方法與過(guò)程、解的表達(dá)形式與應(yīng)用方法都有很大變化[15-18]。因此,準(zhǔn)確概化模型的邊界條件和初始條件,對(duì)解的適用性有著重要作用。
結(jié)合上述實(shí)際,根據(jù)現(xiàn)行降水量、蒸發(fā)量等觀測(cè)制度,將垂向水量交換強(qiáng)度ε(t)離散成逐日變化的階梯函數(shù),建立溝渠水位“穩(wěn)定不變”的邊界條件下、含階梯函數(shù)型源匯項(xiàng)ε(t)的潛水非穩(wěn)定滲流模型。對(duì)模型中Boussinesq 方程,采用第一線性化方法,通過(guò)Laplace 變換得到模型解析解,并對(duì)解析解進(jìn)行數(shù)值驗(yàn)證,討論其物理意義。根據(jù)解析解反問(wèn)題研究思路,建立模型參數(shù)求算方法,相對(duì)文獻(xiàn)[17-18],邊界條件更加符合實(shí)際,垂向水量交換強(qiáng)度處理方法及模型求解過(guò)程更為嚴(yán)謹(jǐn),模型的解也相對(duì)簡(jiǎn)潔、應(yīng)用方便。
2.1 數(shù)學(xué)模型一順直溝渠,如圖1,其所處地段的水文地質(zhì)條件,可概括為:
圖1 溝渠附近潛水滲流場(chǎng)示意
①一側(cè)為溝渠邊界的孔隙潛水含水層均質(zhì)且各向同性、隔水底板水平,在水平方向上無(wú)限延展;
②溝渠在邊界處完整切割潛水含水層,潛水水流可視為一維流;
③潛水初始水位h(x,0)與溝渠水位一致,呈水平狀態(tài);
④溝渠水位在研究期內(nèi)保持不變,可視為一類(lèi)邊界;
⑤垂向水量交換強(qiáng)度ε(t),在研究區(qū)域內(nèi)各處相等。
該問(wèn)題的數(shù)學(xué)模型可寫(xiě)成模型(Ⅰ):
式中:μ為潛水含水層的給水度;K為潛水含水層的滲透系數(shù),m/d;h為潛水水位,m;ε(t)為垂向水量交換強(qiáng)度,m/d;x為計(jì)算點(diǎn)距邊界的距離,m。
在經(jīng)典的J.G.Ferris 模型[18]中,上述水文地質(zhì)條件中的條件④為“河渠水位迅速升高△H后、水位保持穩(wěn)定不變”,而模型(Ⅰ)是“溝渠水位在研究期內(nèi)保持不變”;另外,模型(Ⅰ)在J.G.Ferris 模型基礎(chǔ)上,增加了條件⑤。
2.2 垂向水量交換強(qiáng)度的離散農(nóng)田或邊坡所在地段的潛水含水層,多在野外直接出露于外環(huán)境,可形成的垂向水量交換項(xiàng)有灌溉回歸、降水入滲補(bǔ)給、潛水蒸發(fā)等,相關(guān)觀測(cè)和統(tǒng)計(jì)的變量,對(duì)應(yīng)有灌溉量、降水量、水面蒸發(fā)量,這在實(shí)際工作中都是按時(shí)間步長(zhǎng)為一日進(jìn)行觀測(cè)統(tǒng)計(jì)的,因此,依據(jù)觀測(cè)資料給出的垂向水量交換強(qiáng)度ε(t)一般也是逐日量,即,ε(t)是一個(gè)時(shí)間步長(zhǎng)為一日的階梯函數(shù)(也稱(chēng)分段常值函數(shù))。對(duì)于評(píng)價(jià)期內(nèi)隨時(shí)間變化的垂向水量交換強(qiáng)度ε(t),令評(píng)價(jià)期為n段,第i段時(shí)間為ti-1~ti,對(duì)應(yīng)的垂向水量交換強(qiáng)度為εi。εi為正,表示補(bǔ)給含水層;εi為負(fù),表示排泄含水層。采用“疊加法”將ε(t)寫(xiě)成階梯函數(shù)形式:
式中:ε1為第1 個(gè)時(shí)段(即t0~t1時(shí)段)內(nèi)對(duì)應(yīng)的垂向水量交換強(qiáng)度,m/d;H(t-ti-1)是Heaviside 函數(shù),當(dāng)t<ti-1時(shí)、H(t-ti-1)=0,當(dāng)t≥ti-1時(shí)、H(t-ti-1)=1。
2.3 模型求解在模型(Ⅰ)中,當(dāng)h(x,t)-h(x,0)≤0.1hm(hm為潛水含水層的平均厚度,這在解決實(shí)際問(wèn)題的過(guò)程中基本可以滿(mǎn)足[15-18])時(shí),可以采用Boussinesq 方程第一線性化方法,令u(x,t)=h(x,t)-h(x,0)。
對(duì)模型(Ⅰ)求關(guān)于t的Laplace 變換,注意階梯函數(shù)Laplace 變換的性質(zhì)L[ε(t)]=ε(t)/s[20],L[ε(t)]是關(guān)于ε(t)的Laplace 算符,可得模型(Ⅱ):
式中:a=Khm/μ,a為導(dǎo)壓系數(shù),m2/d;s為L(zhǎng)aplace 算子;uˉ為u關(guān)于t的Laplace 變換過(guò)程中的象函數(shù)。
模型(Ⅱ)中的式(6)是二階常微分方程,由該二階常微分方程的通解,結(jié)合邊界條件式(7)、式(8),可獲得模型(Ⅱ)的特定解:
對(duì)式(9)進(jìn)行Laplace 逆變換,[f(t)]/s項(xiàng)求逆變換時(shí),注意Laplace 變換中的“積分性質(zhì)”,即:
式中er(fz)為誤差函數(shù),。
因?yàn)閡(x,t)=h(x,t)-h(x,0),并將式(5)代入式(11),結(jié)合Heaviside 函數(shù)的性質(zhì),可得:
式(12)即為模型(Ⅰ)的解,也即水位“穩(wěn)定不變”溝渠邊界條件下、含時(shí)變垂向水量交換項(xiàng)ε(t)的潛水非穩(wěn)定滲流模型的解,同時(shí)也是文獻(xiàn)[18]的解在△H=0 時(shí)的特例。
2.4 解析解的數(shù)值驗(yàn)證在研究河渠附近一維潛水非穩(wěn)定滲流問(wèn)題時(shí),針對(duì)不同的實(shí)際問(wèn)題,Zissis和Bansal 等通過(guò)Laplace 變換方法,得到了一維線性化形式的Boussinesq 方程的解析解,利用MacCor?mack 顯式差分方案對(duì)解析解進(jìn)行了數(shù)值驗(yàn)證[21-23]。MacCormack 格式是一種求解可壓縮流體流動(dòng)問(wèn)題的顯式有限差分格式,在時(shí)間和空間上具有二階精度[21-29]。
本文利用Laplace 變換法求解了線性化的Boussinesq 方程,得到了含時(shí)變垂向水量交換項(xiàng)ε(t),水位“穩(wěn)定不變”溝渠邊界條件下潛水非穩(wěn)定滲流模型的解,為了驗(yàn)證線性化方法的有效性和解析解的可靠性,通過(guò)上述研究中MacCormack 格式的有限差分法計(jì)算相應(yīng)非線性方程式(1)的數(shù)值解,并將解析解與數(shù)值解進(jìn)行比較。假設(shè)n時(shí)刻各節(jié)點(diǎn)的水位hk,n均已知,利用MacCormack 格式的預(yù)測(cè)校正兩步法進(jìn)行時(shí)間推進(jìn),求解出n+1 時(shí)刻各個(gè)節(jié)點(diǎn)的水位hk,n+1。首先,在預(yù)測(cè)中對(duì)空間導(dǎo)數(shù)進(jìn)行向前差分,得到用正向差分代替時(shí)空導(dǎo)數(shù)得到h的預(yù)測(cè)方案[21-23],如下:
第二步,對(duì)x的偏導(dǎo)數(shù)用后向差分近似代替,而對(duì)于t的偏導(dǎo)數(shù)用正向差分近似代替,方程式(1)可寫(xiě)為[21-23]:
最終,hk,n+1的修正值為和的算術(shù)平均值[21-23]:
上述式(13)—(15)構(gòu)成了非線性方程式(1)的MacCormack 格式的有限差分法計(jì)算方案,其穩(wěn)定準(zhǔn)則為[21-23]:
根據(jù)式(12)—(16),參考淮北平原區(qū)域的水文地質(zhì)特征,模型選取粗砂、中砂、細(xì)砂三種含水層介質(zhì),給水度分別設(shè)為0.30、0.22 和0.17[30],滲透系數(shù)分別設(shè)為30 m/d、20 m/d 和7 m/d[30],初始潛水位取27.000 m,潛水含水層平均厚度取3.5 m。計(jì)算期為5 天,假設(shè):第一天、第二天為降水入滲階段,對(duì)應(yīng)的降水入滲補(bǔ)給強(qiáng)度分別設(shè)為0.050 m/d 和0.030 m/d;后三天為潛水蒸發(fā)階段,潛水蒸發(fā)強(qiáng)度均設(shè)為0.001 m/d。利用以上參數(shù),24 h和48 h時(shí)的潛水位及相對(duì)誤差計(jì)算結(jié)果分別如圖2、圖3所示。
由圖2 和圖3,在不同的水文地質(zhì)參數(shù)條件下解析解與數(shù)值解的吻合度都比較好,相對(duì)誤差最大不超過(guò)0.15%,表明推導(dǎo)解析解過(guò)程中的線性化方法是有效的,并驗(yàn)證了模型解析解的可靠性,同時(shí)也說(shuō)明了將垂向水量交換強(qiáng)度ε(t)離散成時(shí)間步長(zhǎng)為一日的階梯函數(shù)的方法是可行的。
圖2 t=24h 時(shí)的潛水位解析解與數(shù)值解計(jì)算結(jié)果及相對(duì)誤差
圖3 t=48h 時(shí)的潛水位解析解與數(shù)值解計(jì)算結(jié)果及相對(duì)誤差
在離溝渠邊界30 m 和300 m 處虛擬兩個(gè)潛水位觀測(cè)點(diǎn),分別研究其潛水位變化過(guò)程,計(jì)算結(jié)果如圖4所示。溝渠邊界與垂向水量交換(前兩天的降水入滲及之后的潛水蒸發(fā))共同作用下,降水入滲階段潛水位上升,導(dǎo)致潛水位高于溝渠水位而向溝渠排泄,與遠(yuǎn)處水位觀測(cè)點(diǎn)(x=300 m)相比,近溝渠邊界觀測(cè)點(diǎn)(x=30 m)的潛水位上升速度較慢;潛水蒸發(fā)階段,近溝渠邊界點(diǎn)的潛水位下降速度較快,這表明,溝渠邊界對(duì)潛水位的影響作用隨x增大而減小。由圖4 及式(12),對(duì)于計(jì)算域內(nèi)任意點(diǎn)x,潛水位h始終是隨時(shí)間t變化的,計(jì)算期間計(jì)算域內(nèi)的潛水位始終保持非穩(wěn)定狀態(tài),這符合潛水非穩(wěn)定滲流理論的基本特點(diǎn)。
圖4 距溝渠邊界的距離分別為30m 和300m 處的潛水位變化過(guò)程
2.5 特定解及其物理意義
(1)當(dāng)x→0 時(shí)
此時(shí),對(duì)應(yīng)的是邊界處的潛水位與溝渠水位等高,這符合地下水滲流力學(xué)中的基本原理。
(2)當(dāng)x→∞時(shí)
則有:
此條件下,對(duì)應(yīng)的是距離邊界無(wú)窮遠(yuǎn)處的潛水水位與溝渠邊界無(wú)關(guān),即模型(Ⅰ)中x→∞時(shí)的邊界條件。實(shí)際中,一般離邊界較遠(yuǎn)距離時(shí),溝渠對(duì)潛水水位的影響就基本可以忽略。
(3)當(dāng)n=1 時(shí)
此條件下,計(jì)算期內(nèi)ε1為恒定常數(shù),這也是文獻(xiàn)[17]在ΔH=0 條件下這類(lèi)模型的基本解。
根據(jù)地下水位監(jiān)測(cè)數(shù)據(jù)求解含水層參數(shù),是溝渠附近潛水非穩(wěn)定滲流理論和模型研究的重要應(yīng)用之一。直接體現(xiàn)在式(12)中的參數(shù)有a、μ,通過(guò)a和μ之間的關(guān)系,還隱含著滲透系數(shù)K=μa/hm。
通常情況下,根據(jù)野外常規(guī)試驗(yàn),可以得到參數(shù)μ的值,利用非穩(wěn)定滲流模型,目的是為了求解參數(shù)a。關(guān)于正問(wèn)題研究方法,一般在給定參數(shù)初值的條件下,利用解計(jì)算出水位或水位變動(dòng)速度值,將計(jì)算值與實(shí)際值擬合,通過(guò)不斷試算,以確定參數(shù)值,這一試算過(guò)程顯然比較繁瑣。因此,可利用反問(wèn)題研究思路,建立相應(yīng)的求解方法。
令:
利用n=1 時(shí)段的水位動(dòng)態(tài)數(shù)據(jù),由式(12),得:
式中φ(x,t)為潛水位變動(dòng)速度。
3.1 拐點(diǎn)法由式(21),可得:
由式(22),得:
當(dāng)然,只有n=1 時(shí)段足夠長(zhǎng),才能出現(xiàn)拐點(diǎn)。實(shí)例研究表明,出現(xiàn)拐點(diǎn)的時(shí)間小于24 h,這在實(shí)踐中可以滿(mǎn)足。
3.2 配線法由式(21),對(duì)于一個(gè)到邊界直線距離為x的觀測(cè)孔,z=x/2(at)1/2,首先建立不同a值對(duì)應(yīng)的er(fz)~t曲線圖作為理論曲線圖族。再由實(shí)測(cè)潛水水位,制出曲線,該曲線應(yīng)該可以同erf(z)~t理論曲線圖族中的某條曲線完全重合,僅相差ε/μ倍,此時(shí),兩條曲線所對(duì)應(yīng)的a值相等。因此,根據(jù)該觀測(cè)孔潛水水位實(shí)際觀測(cè)數(shù)據(jù),通過(guò)配線法,可直接確定含水層參數(shù)a的值(如圖6)。
安徽省境內(nèi)淮北平原近淮河地段,潛水含水砂層埋深為5 ~ 6 m,底部為不完全連續(xù)黏性土層?,F(xiàn)狀條件下,該地段潛水埋深多為1.5 ~ 3.0 m,潛水蒸發(fā)(該地段潛水極限蒸發(fā)深度為4.8 m[31])、降水入滲等垂向水量交換作用對(duì)潛水位動(dòng)態(tài)影響明顯。研究區(qū)內(nèi)明溝排水系統(tǒng)中,排水大溝設(shè)計(jì)深度為5 m,基本完全切割潛水含水層,相鄰大溝之間的距離約為2 km[32-33]。受渠首滾水壩控制,排水過(guò)程中溝渠內(nèi)的水位基本保持不變。2009年7月11日和12日,產(chǎn)生了一次時(shí)程分配比較均勻的連續(xù)降雨,在此之前,已存在很長(zhǎng)一段時(shí)間內(nèi)連續(xù)無(wú)降雨的情況,11日、12日的日降雨量分別為86 mm、74 mm,降雨前后,溝渠的水位沒(méi)有變化。在離溝渠直線距離65 m 處有一口地下水位觀測(cè)井(自記式),潛水位隨時(shí)間的變化過(guò)程見(jiàn)表1。平均降水強(qiáng)度為80 mm/d,根據(jù)區(qū)域研究成果,降水入滲補(bǔ)給系數(shù)取0.20,降水入滲補(bǔ)給強(qiáng)度為16 mm/d。
(1)拐點(diǎn)法求導(dǎo)壓系數(shù)a。由表1 中的地下水位監(jiān)測(cè)數(shù)據(jù),用向前插值法計(jì)算和,按時(shí)間間隔為6 h 進(jìn)行摘錄,曲線在18~21 h 之間出現(xiàn)拐點(diǎn),因而,該時(shí)段內(nèi)摘錄的時(shí)間間隔加密至1 h。相關(guān)數(shù)據(jù)和求解過(guò)程見(jiàn)表1,計(jì)算得到的曲線變化過(guò)程見(jiàn)圖5。
圖5 拐點(diǎn)法求tg的過(guò)程
據(jù)上述有關(guān)基礎(chǔ)數(shù)據(jù)與拐點(diǎn)時(shí)間:x=65 m,tg=19.5 h,由式(24),可求出a=865 m2/d。
(2)配線法求導(dǎo)壓系數(shù)a。計(jì)算過(guò)程中,μ通過(guò)野外實(shí)測(cè)數(shù)據(jù)獲得,μ的值為0.035。由表1 中的數(shù)據(jù),計(jì)算出的φ(x,t)~t曲線與x=65 m 處的er(fz)~t理論曲線圖中的曲線組配線。根據(jù)觀測(cè)數(shù)據(jù)得到的φ(x,t)~t曲線處于a值為800 m2/d 和1000 m2/d 的理論曲線之間,并靠近a值為800 m2/d 的理論曲線(見(jiàn)圖6)。由此,含水層參數(shù)a可取為850 m2/d,這同拐點(diǎn)法的計(jì)算結(jié)果僅相差1.73%。
圖6 配線法求a 的過(guò)程
表1 x=65m 處地下水位動(dòng)態(tài)監(jiān)測(cè)數(shù)據(jù)(2009.7.11—2009.7.12)與整理
理論曲線族繪制過(guò)程中,可根據(jù)實(shí)際計(jì)算精度要求對(duì)理論曲線進(jìn)行加密,從而提高配線法求解導(dǎo)壓系數(shù)的精度。
(3)參數(shù)反演結(jié)果的合理性分析。通過(guò)數(shù)值法驗(yàn)證了模型解析解式(12)的可靠性,因此可將上述計(jì)算得到的導(dǎo)壓系數(shù)a的值代入式(12),計(jì)算x=65 m 處不同時(shí)刻的潛水位,并與實(shí)測(cè)值進(jìn)行比較,如圖7所示。
圖7 x=65m 處潛水位解析解與實(shí)測(cè)值對(duì)比結(jié)果
對(duì)比分析,當(dāng)a=865 m2/d 和a=850 m2/d 時(shí),x=65 m 處潛水位的解析解與實(shí)測(cè)值都非常接近,表明拐點(diǎn)法和配線法求算導(dǎo)壓系數(shù)a的方法是合理、可靠的。
在類(lèi)似淮北平原沿淮潛水位埋深較淺的區(qū)域,溝渠附近的地下水水位變動(dòng),不僅與灌溉回歸、降水入滲補(bǔ)給、潛水蒸發(fā)等垂向水量交換作用有關(guān),而且受溝渠等排水工程影響明顯。在研究垂向水量補(bǔ)排作用下溝渠附近潛水非穩(wěn)定滲流模型的過(guò)程中,形成以下主要結(jié)論:
(1)在潛水位埋深較淺的區(qū)域,用于控制潛水位的排水溝渠,對(duì)潛水水位消退變化的影響顯然有控制性作用,通常潛水蒸發(fā)、降水入滲補(bǔ)給等垂向水量補(bǔ)排作用也不可忽視,此條件下,宜將這兩者統(tǒng)一納入數(shù)學(xué)模型中。
(2)對(duì)于有滾水壩調(diào)節(jié)控制的溝渠,在自由排水期間,溝渠內(nèi)的水位一般保持穩(wěn)定不變。根據(jù)實(shí)際觀測(cè)統(tǒng)計(jì)制度,潛水蒸發(fā)、降水入滲補(bǔ)給等垂向水量補(bǔ)排作用強(qiáng)度,宜概化為逐日變化的階梯函數(shù)。
(3)利用模型解反問(wèn)題研究思路,根據(jù)潛水位變動(dòng)速度隨時(shí)間的變化規(guī)律,建立拐點(diǎn)法、配線法求模型參數(shù)的方法,直接簡(jiǎn)明;利用文中模型的解,通過(guò)反問(wèn)題研究方法,也可建立依據(jù)潛水位變動(dòng)速度計(jì)算垂向水量交換強(qiáng)度的算式。
(4)值得指出的是,實(shí)際工作中,大范圍農(nóng)田排水往往設(shè)置有多條平行排水溝渠,本文方法適用于溝渠相距足夠遠(yuǎn)的情形;當(dāng)溝渠相距較近并對(duì)地下水水位有疊加影響時(shí),屬兩側(cè)都是一類(lèi)水位邊界控制下的潛水非穩(wěn)定流問(wèn)題,此情形下問(wèn)題求解方法不同于本文方法,值得進(jìn)一步研究。