周燕怡，王旭升

(中國地質(zhì)大學(xué)(北京)水資源與環(huán)境學(xué)院,北京 100083)

潛水蒸發(fā)是指淺層地下水在包氣帶向上運(yùn)移、貢獻(xiàn)給陸面蒸散的過程，屬于地下水的一種排泄方式。潛水蒸發(fā)的機(jī)理涉及液態(tài)水向上運(yùn)移和表土蒸發(fā)2個(gè)過程[1]，后者也可以稱為土體中的水分蒸發(fā)[2]。土體中植被根系吸水產(chǎn)生的蒸騰作用也是陸面蒸散的重要組成部分，屬于土壤-植物-大氣連續(xù)體(SPAC)過程[3-4]。潛水蒸發(fā)可以通過水分和鹽分的輸送影響SPAC過程，因此又被擴(kuò)展為GSPAC過程[5]。在干旱區(qū)，潛水蒸發(fā)是地下水排泄的重要方式，它不僅向地表和大氣輸送水分，而且把地下水循環(huán)過程中攜帶的物質(zhì)(包括潛在污染物)輸送到地表產(chǎn)生環(huán)境影響，例如“水走鹽留”造成土壤次生鹽漬化。潛水蒸發(fā)的強(qiáng)度隨著地下水位的加深而減弱，其關(guān)系可以用簡化公式表示，但屬于地區(qū)經(jīng)驗(yàn)，取決于具體的條件。一般認(rèn)為存在一個(gè)極限埋深[6]，當(dāng)?shù)叵滤簧疃瘸^臨界值時(shí)，潛水蒸發(fā)為零。受地形起伏影響，地下水埋深在很小的空間范圍也可能發(fā)生較大的變化。因此，大面積區(qū)域潛水蒸發(fā)的評估需要考慮不同空間尺度的地形變化。Li[7]等以新疆焉耆盆地為例進(jìn)行分析，發(fā)現(xiàn)區(qū)域潛水蒸發(fā)量的計(jì)算結(jié)果在一定程度上取決于數(shù)字高程的柵格分辨率。地下水鹽度[8]和季節(jié)性凍結(jié)[9]也對潛水蒸發(fā)量有影響。

內(nèi)蒙古西北部的巴丹吉林沙漠屬于典型的極端干旱區(qū)。與干旱氣候形成反差的是，這片沙漠分布大量的湖泊。研究者指出該沙漠湖泊的維持依賴于地下水排泄[10-12]。關(guān)于沙漠地下水的來源，很多研究者提出了一些富有爭議的推測[13-15]。Dong等[16]將這些推測分類為4種假說。導(dǎo)致地下水來源爭議的另一個(gè)焦點(diǎn)問題就是地下水的消耗，即巴丹吉林沙漠每年消耗多少地下水。不少研究者根據(jù)氣象數(shù)據(jù)計(jì)算了沙漠湖泊的蒸發(fā)耗水量[16]，但忽略了干旱區(qū)常見的潛水蒸發(fā)耗水。陳添斐[17]在進(jìn)行巴丹吉林沙漠湖泊群地下水模擬時(shí)，發(fā)現(xiàn)無湖洼地的潛水蒸發(fā)總量可能與湖面蒸發(fā)總量相當(dāng)。這意味著潛水蒸發(fā)在巴丹吉林沙漠具有突出的作用。不過，該模擬結(jié)果是使用線性蒸發(fā)假設(shè)、取極限埋深為1.5 m的情況下得到的，并不一定能夠代表沙漠潛水蒸發(fā)的真實(shí)特征。

為了更加準(zhǔn)確地評估巴丹吉林沙漠地下水的消耗量，必須弄清楚潛水蒸發(fā)強(qiáng)度與地下水埋深之間的關(guān)系。受到環(huán)境條件和儀器設(shè)備的限制，在巴丹吉林沙漠直接觀測不同水位深度的潛水蒸發(fā)非常困難。為此，本文以巴丹吉林沙漠的基本氣候和砂土環(huán)境特征為背景，開展不同地下水埋深條件下的SPAC過程數(shù)值模擬，力圖提供更加可靠的定量依據(jù)用于計(jì)算沙漠潛水蒸發(fā)量。

1 研究區(qū)概況

巴丹吉林沙漠位于內(nèi)蒙古阿拉善高原，面積約5×104km2，其東南部有高大沙山和湖泊交錯(cuò)分布的奇特景觀。沙漠中沒有地表河流，因此地下水是湖泊水分的唯一來源。作為干旱區(qū)，巴丹吉林沙漠多年平均降水量較少，介于50～150 mm，自東南向西北降水量逐漸減少，一年之中降水主要發(fā)生在6—8月。與降水量偏少的特征相反，沙漠周邊氣象站觀測到的水面蒸發(fā)量很大，多年平均蒸發(fā)量均超過2 500 mm(小型蒸發(fā)皿)。在地質(zhì)構(gòu)造上，巴丹吉林沙漠屬于阿拉善地塊，發(fā)育斷陷沉積盆地，厚度超過500 m的孔隙和裂隙含水層廣泛發(fā)育[18]。沙漠表層主要被第四系風(fēng)積沙覆蓋，風(fēng)積砂具有較強(qiáng)分選性，粒徑在0.075 ～2 mm范圍內(nèi)的顆粒占80%以上，粉黏粒含量很低，屬于細(xì)砂[19]，孔隙度0.35～0.42。根據(jù)前人研究[20]，風(fēng)積沙的天然干容重一般為1.48～1.68 g/cm3，內(nèi)摩擦角在31°～39°，壓縮系數(shù)小于0.1 MPa-1。

巴丹吉林沙漠植被主要沿著湖泊分布，在地下水埋藏較淺的洼地也有分布。植被在局部地形上的空間變化反映了潛水蒸發(fā)與地下水埋深之間的關(guān)系。圖1以蘇木吉林湖區(qū)為例，反映了沙丘地下水和植被分布的剖面特征。在靠近湖泊的地帶，地下水埋深變淺甚至形成侵蝕下降泉，流動(dòng)水體減弱了鹽堿化問題，有利于形成草地。但是，在蘇木吉林北湖的南部，地勢平坦，地下水埋藏淺，強(qiáng)烈的潛水蒸發(fā)作用導(dǎo)致了鹽分的累積，形成鹽堿地，僅有喜鹽草類稀疏生長?？傮w上，從湖岸到沙丘高處，植被景觀呈現(xiàn)草地—稠密灌叢—稀疏灌叢或裸沙的分帶變化，與地下水埋深的分帶變化相對應(yīng)。

圖1 蘇木吉林湖區(qū)水土環(huán)境特征典型剖面圖[12]Fig.1 Schematic profile of the soil-water environment across the Sumujilin Lakes[12]

巴丹吉林沙漠潛水蒸發(fā)作用反映在土體含水率隨深度的變化。筆者在蘇木吉林湖區(qū)使用土壤水分傳感器TDR進(jìn)行了夏季不同深度含水率的原位觀測，結(jié)果見圖2。圖中3個(gè)觀測場地代表了地下水埋深在0.5～1.5 m之間不同情況。DS-11和DA-1測點(diǎn)的地下水埋深都小于1.0 m，含水率隨著深度加大而增加。土壤含水率越低，毛細(xì)水吸力越大，在剖面上形成了持續(xù)從潛水面吸收水分到地表蒸發(fā)的水勢梯度。在DW-3測點(diǎn)，地下水埋深接近1.3 m，總體上存在隨深度加大含水率增加的趨勢。但是在剖面深度1.0 m范圍內(nèi)出現(xiàn)了含水率的2個(gè)峰值。第1個(gè)峰值大致位于深度0.1 m處，第2個(gè)峰值大致位于深度0.4 m處，含水率都低于20%。由于風(fēng)積砂比較均勻，這種含水率的變化不太可能是土質(zhì)變化導(dǎo)致的，而是觀測之前若干次強(qiáng)降雨入滲形成的，造成了水勢剖面上的零通量面，峰值點(diǎn)上部水分向上運(yùn)移，峰值點(diǎn)下部水分向下運(yùn)移。這意味著當(dāng)?shù)叵滤裆钶^大時(shí)，潛水蒸發(fā)的作用較弱，強(qiáng)降雨形成入滲可以在某個(gè)時(shí)期內(nèi)超過蒸發(fā)作用導(dǎo)致包氣帶水分向下運(yùn)移。這種蒸發(fā)和入滲的交疊作用增加了潛水蒸發(fā)過程的復(fù)雜性。

圖2 巴丹吉林沙漠地下水淺埋區(qū)土壤含水率隨深度的變化特征[21]Fig.2 Depth-dependent soil moisture content at typical sites with shallow groundwater in the Badain Jaran Desert[21]

圖3 潛水蒸發(fā)一維土柱概念模型Fig.3 1-D conceptual model for estimating groundwater evaporation

2 潛水蒸發(fā)模型設(shè)計(jì)

2.1 概念模型和模擬工具

潛水蒸發(fā)涉及地下水與SPAC過程的相互作用，需要考慮SPAC的各種要素，包括土壤、植物和大氣邊界。SPAC過程以垂向水分交換為主。本文建立一個(gè)受地下水影響的垂向一維SPAC模型，研究不同地下水位埋深情況下的潛水蒸發(fā)問題。模擬工作采用Hydrus-1D軟件[22]。

如圖3所示，模型的底部為潛水面，處理為定水頭邊界，上部是一個(gè)SPAC系統(tǒng)，以大氣邊界為驅(qū)動(dòng)條件。SPAC的包氣帶發(fā)生土壤水分的非飽和滲流運(yùn)動(dòng)，采用一維形式的Richards方程描述[23]：

(1)

式中：z——垂向坐標(biāo)/m, 以向上為正；

t——時(shí)間/d；

h——土壤水的壓力水頭/m；

K(h)——隨壓力水頭變化的滲透系數(shù)/(m·d-1)；

C(h)=?θ/?h——土壤容水度/m-1，反映的是含水率θ隨壓力水頭變化的特征；

θ——隨h的變化采用土水特征曲線表示；

S(z)——隨深度變化的單位體積根系吸水強(qiáng)度/(d-1)。

模型的邊界條件，地面邊界：

(2)

潛水面邊界：

h=0，z=-d，t>0

(3)

式中：P——降水強(qiáng)度/(m·d-1)；

Es——土面蒸發(fā)強(qiáng)度/(m·d-1)；

d——地下水位的埋深/m。

非飽和土壤的壓力水頭為負(fù)值，h=0表示潛水面處砂土飽和、壓力水頭為零。Hydrus-1D模型對蒸散發(fā)的模擬有專門的處理方法。首先，蒸散總量被分解為土面蒸發(fā)和蒸騰2個(gè)部分，并且都與潛在蒸散量成正比[24]：

Ep=ETpe-ηLAI

(4)

Tp=ETp(1-e-ηLAI)

(5)

式中：ETp——潛在蒸散強(qiáng)度/ (m·d-1)，取決于氣象條件；

Ep和Tp——土面蒸發(fā)和植被蒸騰強(qiáng)度的最大可能值/ (m·d-1)；

LAI——葉面積指數(shù)；

η——消光系數(shù)(本文按經(jīng)驗(yàn)取η=0.6)。

在2個(gè)潛在蒸散分量中，Ep控制土面蒸發(fā)Es的計(jì)算，Tp控制根系吸水強(qiáng)度S(z)的計(jì)算。其中，Es的計(jì)算需要考慮表層干砂的最大吸力，即-hmax。按照Hydrus-1D推薦的做法，本文采用土水特征曲線中含水率為θr+0.005時(shí)對應(yīng)的壓力水頭確定最大吸力(θr為砂土殘余含水率)。在Hydrus-1D中，根系吸水采用隨根系密度變化的分配函數(shù)和特定的水分脅迫函數(shù)計(jì)算。本研究采用Hydrus-1D嵌入的S-shaped模型[25]作為根系吸水的水分脅迫函數(shù)，該模型包含P0和P50兩個(gè)參數(shù)。P0為衰減指數(shù)，取Hydrus-1D建議的經(jīng)驗(yàn)值(P0=3)。P50指根系吸水強(qiáng)度衰減到50%對應(yīng)的壓力水頭，取經(jīng)驗(yàn)值(P50=-1 m)。在根系層，根系密度也是隨深度變化的(圖3)，與具體的植被景觀特征有關(guān)。本文將在植被情景設(shè)計(jì)中對此加以處理。模型上邊界的驅(qū)動(dòng)條件采用重復(fù)年周期波動(dòng)的逐日降水量和潛在蒸散量數(shù)據(jù)，因?yàn)楸疚闹饕芯磕承夂蚯榫跋碌亩嗄昶骄芷谛孕袨?取周期長度為365 d)。模型初始條件假設(shè)為壓力水頭與潛水面形成靜水力學(xué)平衡狀態(tài)，在反復(fù)周期模擬中，初始條件的影響將逐步消失，達(dá)到重復(fù)性的周期變化狀態(tài)。

本文模型在垂向上采用均勻的剖分方式，離散節(jié)點(diǎn)的間距為10 cm。測試結(jié)果表明更高的模型分辨率并不會顯著改善模擬結(jié)果的精度，因此10 cm的節(jié)點(diǎn)間距對本文研究是可行的。在時(shí)間剖分上，根據(jù)數(shù)值計(jì)算的收斂特征，采用Hydrus-1D的變步長迭代法，最小時(shí)間步長為0.000 01 d，最大時(shí)間步長為5 d。模型含水率迭代精度取0.001。每一種氣候情景下的模擬時(shí)間為20年，取最后1年的結(jié)果代表多年平均周期狀態(tài)。對比觀察不同時(shí)期的模擬結(jié)果，表明模擬時(shí)間達(dá)到10年以上則初始狀態(tài)不確定性導(dǎo)致的年際變化均小于5%，因此取第20年的結(jié)果足以消除初始條件的影響。

取得多年平均狀態(tài)周期解的逐日模擬數(shù)據(jù)后，計(jì)算多年平均潛水蒸發(fā)量：

(6)

式中：Eg——多年平均潛水蒸發(fā)量/mm；

i——一個(gè)周期年內(nèi)的天數(shù)；

Vi——潛水面處即模型底部的逐日水分通量/(mm·d-1)，向上為正；

δ——階躍函數(shù)，即Vi>0時(shí)δ=1，否則δ=0；

Δti——累加時(shí)間步長，取Δti=1 d。

使用階躍函數(shù)的目標(biāo)，是單純提取出潛水面對蒸散發(fā)的貢獻(xiàn)，而篩除掉降水入滲對地下水的補(bǔ)給。

2.2 氣候情景

巴丹吉林沙漠的周邊分布有多個(gè)氣象站，具備長期氣象觀測數(shù)據(jù)，而沙漠內(nèi)部只在近幾年才開始?xì)庀笥^測。根據(jù)周邊氣象資料，巴丹吉林沙漠降水量具有東南高、西北低的特點(diǎn)，多年平均降水量為40～150 mm。小型蒸發(fā)皿(E20型)水面蒸發(fā)量與降水量分布趨勢相反，東南部為3 000～3 500 mm/a，北部可達(dá)4 000 mm以上。2012年，中國地質(zhì)大學(xué)(北京)在沙漠腹地的蘇木吉林南湖架設(shè)自動(dòng)氣象站[12]，對降水量和湖面蒸發(fā)量(E601型蒸發(fā)皿)進(jìn)行觀測。該自動(dòng)氣象站的數(shù)據(jù)能更準(zhǔn)確地反映沙漠湖泊區(qū)氣象條件，作為本次模擬研究的基礎(chǔ)數(shù)據(jù)，其中2013年完整年的逐日降水量和蒸發(fā)量數(shù)據(jù)見圖4。該年的總降水量約為157.4 mm，總蒸發(fā)量約為1 314 mm。本文取該氣象站E601型蒸發(fā)皿觀測的水面蒸發(fā)量為潛在蒸散量(ETp)，替代其它氣象數(shù)據(jù)計(jì)算潛在蒸散量的簡化方法。

圖4 蘇木吉林湖區(qū)2013年氣象要素逐日變化特征Fig.4 Daily meteorological data in the Sumujilin lake area in 2013

為了反映沙漠氣候環(huán)境的時(shí)空變化，在設(shè)計(jì)模型的氣候情景時(shí)，取參考年的降水量代表沙漠相對濕潤氣候情景(用符號C1表示)。通過去除一些雨強(qiáng)較小的降雨事件，使總降水量降低到100 mm，作為沙漠的平均氣候情景(用符號C2表示)。進(jìn)一步去除更多雨強(qiáng)較小的降雨事件，使總降水量降低到50 mm，作為沙漠相對干燥氣候情景(用符號C3表示)。潛在蒸散量如何隨干濕環(huán)境變化是一個(gè)比較復(fù)雜的問題，目前缺少定論。為便于對比不同地下水埋深對潛水蒸發(fā)的影響，本文假設(shè)潛在蒸散量在各個(gè)氣候情景中保持不變。

2.3 砂土參數(shù)

砂土的非飽和土水特征曲線一般用相對飽和度Se作為參考：

(7)

式中：θr和θs——?dú)堄嗪逝c飽和含水率。

在Hydrus-1D中，推薦使用V-G公式[26]表示土壤含水率與壓力水頭之間的關(guān)系，并以此計(jì)算不同相對飽和度下的滲透系數(shù)，即：

(8)

(9)

式中：α——與土壤進(jìn)氣值有關(guān)的參數(shù)/m-1；

n——無量綱的參數(shù)，m=1-(1/n)；

Ks——飽和滲透系數(shù)/(m·d-1)；

l——孔隙彎曲度的參數(shù)，依據(jù)經(jīng)驗(yàn)可取為0.5。

實(shí)際的土水特征曲線還需要考慮吸濕和脫濕過程的滯后效應(yīng)，對短期行為和微觀機(jī)理研究比較重要。本文主要研究長周期行為，這種滯后效應(yīng)暫不予以考慮。

對于巴丹吉林沙漠，已經(jīng)通過采樣法進(jìn)行了砂土水分特征曲線的研究[19, 21, 27-29]，包括壓力膜儀法和離心機(jī)法測試水分特征曲線獲取的VG公式擬合參數(shù)(θr、θs、α、n)。這些參數(shù)都有一定的變異性，其中θs的變異性最小。此外，采用單環(huán)入滲儀在沙丘表層進(jìn)行了飽和滲透系數(shù)的觀測[30]，發(fā)現(xiàn)風(fēng)積砂的Ks為0.5～75 m/d，有一個(gè)較大的變化范圍。實(shí)際上，巴丹吉林沙漠的第四系沉積物有一定的空間變異性，普遍為細(xì)砂，局部地區(qū)也存在中砂、粗砂沉積。關(guān)于砂土參數(shù)的變異性，前人進(jìn)行了一些定量研究[31-32]。通過收集大量土壤樣品(包括砂、黏土、壤土以及粉土)數(shù)據(jù)，Carsel和Parrish[32]建立了θr、Ks、α和n等參數(shù)的耦合隨機(jī)模型，其中含有砂土的概率統(tǒng)計(jì)參數(shù)，與巴丹吉林沙漠的砂土特征具有相似性。筆者參考該方法對砂土參數(shù)進(jìn)行Monte Carlo模擬，得到5 000個(gè)砂土參數(shù)的隨機(jī)組并分析了包氣帶水下滲的不確定性[33]。本文以研究地下水埋深的影響為重點(diǎn)，在砂土參數(shù)變異性方面做了簡化處理，即設(shè)計(jì)出3個(gè)砂土參數(shù)組合分別代表細(xì)砂、中砂和粗砂的情景。各種情景下的V-G公式參數(shù)見表1。由于不同砂土類型的θs變異性并不大，在此統(tǒng)一取最大值0.42。

表1 砂土特征參數(shù)的三種組合情景

2.4 地下水埋深和植被情景

在巴丹吉林沙漠干旱區(qū)，植被覆蓋特征對地下水埋深具有高度依賴性(圖1)，在模型中的具體表現(xiàn)之一就是葉面積指數(shù)的季節(jié)性變化特征隨地下水埋深而變化。杜占池[34]等對內(nèi)蒙古典型草原地區(qū)沙地植物的葉面積指數(shù)(LAI)進(jìn)行了研究，發(fā)現(xiàn)LAI最大值通常出現(xiàn)于8月份。唐思凌[35]等對典型沙生灌木的研究表明沙蒿生長季為4月底到10月下旬，4月底—7月初LAI為穩(wěn)步上升，達(dá)到最大值的10%～90%，7月初到9月中旬為穩(wěn)定期LAI接近最大值，9月中旬到10月下旬LAI迅速降低。本文參考這些研究結(jié)果，在潛水蒸發(fā)模型中將植被生長季確定為5—10月，各月LAI按照固定的比例周期變化，分別為最大LAI的30%、70%、100%、100%、75%、50%，非生長季的LAI取為零?？紤]到隨著地下水埋深的加大，植被總體上向稀疏和無覆蓋狀態(tài)轉(zhuǎn)變，模型的最大LAI隨水位埋深情景(0.2～3.0 m)的不同而在2.5～0.1之間變化(表2)。

植被根系吸水與地下水埋深有關(guān)。有研究者對沙漠草地根系特征進(jìn)行研究[36-37]，表明沙地草本植物的根系主要分布在深度0.5 m以內(nèi)，蒿類的根系可以達(dá)到1.0 m，根系密度均勻分布或隨深度線性減小，淺根一般均勻分布。前人對油蒿和檸條的研究[38]表明沙漠灌木最大根系深度可達(dá)1.5～3.0 m，最大根系密度一般出現(xiàn)在0.3～0.5 m，以下根系密度隨深度增加近似指數(shù)衰減。結(jié)合已有研究，本文假設(shè)植被根系密度在地下水埋深0.2 m時(shí)為均勻分布，在地下水埋深為0.5～1.0 m時(shí)隨深度線性減小，在地下水埋深1.5～3.0 m時(shí)隨深度呈中部大兩頭小的方式變化。同時(shí)，按照Hydrus-1D的要求，確保根系密度函數(shù)在整個(gè)根系層的積分為1。模型中不同深度段的根系密度見表2。

表2 不同模擬情景的地下水埋深和植被特征

本文假設(shè)根系的深度和密度函數(shù)并不隨季節(jié)變化，只有LAI隨季節(jié)變化，代表特定氣候背景下植被種群穩(wěn)定分布的狀態(tài)。另外，本文暫不考慮鹽分運(yùn)移對潛水蒸發(fā)和植被的影響。地下水埋深淺(d<0.5 m)時(shí)可能發(fā)育植被稀疏的鹽堿地，但這時(shí)的總蒸散量不管有無植被都很強(qiáng)，不會顯著影響對潛水蒸發(fā)的計(jì)算結(jié)果。

本次研究包括3種氣候情景(C1, C2, C3)，3種砂土參數(shù)情景(S1, S2, S3)和6種地下水埋深(最大埋深為3.0 m)情景，一共設(shè)計(jì)有54個(gè)情景模型。實(shí)際上，為了尋找潛水蒸發(fā)的極限埋深，對地下水埋深大于3.0 m的情況也進(jìn)行了模擬，即試探性的將模型深度調(diào)整為3.1 m、3.2 m、……, 直到在年周期的任何時(shí)刻都不存在潛水蒸發(fā)為止。

3 模擬結(jié)果分析

3.1 潛水蒸發(fā)強(qiáng)度的季節(jié)性變化特征

Hydrus-1D能夠計(jì)算逐日的模型底部水分通量，根據(jù)這些數(shù)據(jù)可以觀察潛水蒸發(fā)強(qiáng)度的季節(jié)性變化。圖6給出了平均氣候情景(C2)和細(xì)砂(S1)情景下的模型底部水分通量V(式6)的變化曲線，包括地下水埋深為1 m和2 m兩種情況的模擬結(jié)果。當(dāng)V>0時(shí)，模型底部水分向上運(yùn)移，說明發(fā)生了潛水蒸發(fā)；當(dāng)V<0時(shí)，則水分向下運(yùn)移，包氣帶水下滲補(bǔ)給地下水。從圖6可以看出，當(dāng)?shù)叵滤裆顬? m時(shí)，多數(shù)時(shí)間V>0，即以潛水面蒸發(fā)為主，只有少數(shù)發(fā)生降水的日期及雨后的若干天內(nèi)出現(xiàn)V<0的情況，而且呈現(xiàn)一種顯著的鋸齒狀變化。這說明當(dāng)?shù)叵滤裆畋容^淺的時(shí)候，強(qiáng)降水比較容易轉(zhuǎn)化為地下水補(bǔ)給。對于V>0時(shí)，d=1 m情景的潛水蒸發(fā)強(qiáng)度存在顯著的季節(jié)性振蕩，在7月份達(dá)到峰值，與潛在蒸散量(圖4)的變化趨勢基本一致。相比之下，地下水埋深為2 m的情景顯示出非常平緩的V值變化，似乎只有3次降水事件引發(fā)了潛水面處的下滲補(bǔ)給，但是每次V<0持續(xù)的時(shí)間比d=1 m的情景更長一些，而開始出現(xiàn)V>0的時(shí)刻均落后于強(qiáng)降雨發(fā)生時(shí)刻。這實(shí)際上反映了包氣帶對降水-入滲補(bǔ)給過程的滯后延遲作用，包氣帶厚度越大，滯后效應(yīng)越顯著。不過，僅觀察逐日變化曲線難以判斷地下水的年補(bǔ)給量如何隨潛水面深度而變化。雖然d=1 m時(shí)，響應(yīng)降雨事件的V絕對值要大于d=2 m的情景，似乎意味著地下水獲得了更多的補(bǔ)給，但其潛水蒸發(fā)的強(qiáng)度遠(yuǎn)大于d=2 m的情景，說明消耗的地下水也很多。只有通過累計(jì)年蒸發(fā)量與年降雨量比較才能最終判斷實(shí)際地下水補(bǔ)給量與消耗量之間的關(guān)系。

圖5 平均氣候(C2)和細(xì)砂(S1)情景下底部水分通量(V)以及降水(P)的逐日變化(d表示地下水埋深)Fig.5 Daily variation patterns of the bottom flux (V) in the model and the precipitation (P) for the scenarios of average climate conditions (C2) and fine sands (S1), where d is the depth to water table

3.2 多年平均潛水蒸發(fā)量與地下水埋深的關(guān)系

為了在年尺度上觀察地下水與SPAC系統(tǒng)的水量交換，需要把逐日水分通量的模擬結(jié)果在周期年內(nèi)進(jìn)行累計(jì)，例如多年平均潛水蒸發(fā)量按照式(6)進(jìn)行計(jì)算。采用類似的累計(jì)方法可以計(jì)算出多年平均的土面蒸發(fā)量(Es)和蒸騰量(Ta)，實(shí)際蒸散量為兩者之和，即ETa=Es+Ta。將地下水的凈補(bǔ)給計(jì)算為：

Rg=P-ETa

(10)

其中，Rg代表多年平均狀態(tài)下的地下水凈補(bǔ)給，即超過年蒸散量的部分年降水量轉(zhuǎn)化為地下水補(bǔ)給。如果Rg<0，意味著降水量不足以支撐蒸散量，必須通過潛水蒸發(fā)補(bǔ)充。圖6針對9種典型情景，給出了關(guān)鍵水分通量年累計(jì)值隨地下水埋深的變化特征?？傮w看，當(dāng)潛水面深度小于2 m時(shí)，這些水分通量對地下水埋深的變化都比較敏感，而且不同情景的數(shù)據(jù)顯示出一定的差異；潛水面深度超過2 m時(shí)，水分通量趨于穩(wěn)定，隨地下水埋深的變化只有少量的調(diào)整。

圖6 模型水分通量要素多年平均值隨地下水埋深的變化Fig.6 Variations in the mean annual values of the simulated water flux components

由圖6(a)可知地下水凈補(bǔ)給量(Rg)可正可負(fù)，當(dāng)d<2時(shí)，總是有Rg<0，表明蒸散消耗的水量總是大于降水量，需要潛水蒸發(fā)進(jìn)行補(bǔ)充。在相同的砂土參數(shù)情景和地下水埋深情況下，氣候越干旱(從C1到C3)，Rg越小，即潛水蒸發(fā)的貢獻(xiàn)越大。相比之下，砂土類型變化引起的差異較小。當(dāng)?shù)叵滤裆钶^大時(shí)，可以獲得真正的入滲補(bǔ)給(Rg>0)，但在最好的情景下年補(bǔ)給量也不超過40 mm，說明沙漠的地下水補(bǔ)給受到很大的限制。

土面蒸發(fā)量隨地下水埋深的變化不僅具有非線性特征，而且在d為0.5～1.0 m的區(qū)間會出現(xiàn)峰值(圖6b)，變化曲線整體上呈反S型。這種特征是2個(gè)相反的趨勢疊加造成的。首先，隨著地下水埋深變淺，表土變濕潤，蒸發(fā)強(qiáng)度向最大可能土面蒸發(fā)量Ep靠攏。其次，地下水埋深變淺導(dǎo)致LAI增大產(chǎn)生強(qiáng)烈的遮光效應(yīng)，按照式(4)分配的Ep值變小。這2個(gè)相反變化趨勢的疊加產(chǎn)生了Es的峰值。然而，這種疊加效應(yīng)對蒸騰量的變化曲線(圖6c)幾乎沒有影響，而且不管情景如何，蒸騰量的變化都近似為同一條曲線，似乎說明根系吸水的自我調(diào)節(jié)作用使得實(shí)際蒸騰量只隨最大葉面積指數(shù)變化。不過，需要意識到本文模型沒有考慮發(fā)生鹽堿化的情況，因此當(dāng)d<0.5 m時(shí)，單獨(dú)的土面蒸發(fā)量或蒸騰量模擬結(jié)果都不能代表實(shí)際可能發(fā)生的情況，只有兩者的總和(即ETa)才有意義。

多年平均潛水蒸發(fā)量(Eg)隨地下水埋深的變化顯示在圖6(d)中，可以看出Eg隨地下水埋深增大而減小的趨勢具有顯著的非線性特征，并且存在使得Eg=0的極限埋深?？傮w而言，各種情景下的極限埋深都大于3 m，而且在d=3 m時(shí)潛水蒸發(fā)量都小于最大值的5%。當(dāng)d在0.5～1.5 m區(qū)間，Eg對地下水埋深的變化十分敏感，反映了砂土毛細(xì)水平均上升高度一般小于1.0 m的特征。當(dāng)d=1.0 m時(shí)，Eg對砂土類型的變化也非常敏感，細(xì)砂(S1)的潛水蒸發(fā)量遠(yuǎn)遠(yuǎn)大于粗砂(S2)的情景，相比之下，干濕氣候情景變化對潛水蒸發(fā)量的影響是微弱的。然而，在d<0.5 m和d>1.5 m的區(qū)間，潛水蒸發(fā)量對地下水埋深、氣候情景和砂土類型的變化都不太敏感。

3.3 經(jīng)驗(yàn)公式和極限埋深

圖6(d)潛水蒸發(fā)量隨地下水埋深的變化曲線近似呈S型，與土壤含水率隨壓力水頭的關(guān)系曲線有較高的相似度。筆者試圖采用線性、拋物線型或指數(shù)型的經(jīng)驗(yàn)公式對曲線進(jìn)行擬合，效果都不理想，會產(chǎn)生很大的誤差?？紤]到V-G公式成功地刻畫了S型土水特征曲線，本文提出采用一個(gè)類似V-G公式的關(guān)系式擬合潛水蒸發(fā)曲線：

(11)

其中，β和k都是無量綱參數(shù)，且k>1。dmax為極限埋深，極限埋深的數(shù)值根據(jù)模擬結(jié)果直接確定，因此對于每種情景只需要優(yōu)化出β和k兩個(gè)參數(shù)的取值即可。式(11)可以很好地?cái)M合各種情況下的潛水蒸發(fā)曲線，誤差均小于0.1，而且可以比較準(zhǔn)確地反映在d/dmax處于0.2～0.4區(qū)間時(shí)曲線的急劇變化特征。有關(guān)參數(shù)見表3，典型的擬合曲線特征見圖7。

表3 不同情景下潛水蒸發(fā)特征經(jīng)驗(yàn)式(11)的參數(shù)

圖7 式(11)擬合潛水蒸發(fā)曲線的效果對比圖Fig.7 Comparison between the fitting curves by using Equation (11) and the modeling data points

從表3可以看出，不同情景的極限埋深在3.3～3.8 m變化，平均值接近于3.5 m，極限埋深總體上隨著土質(zhì)變粗(從S1到S3)和氣候情景向干燥方向遷移均不是單調(diào)的增減變化，其機(jī)理尚不清楚。參數(shù)β的變化比較有規(guī)律，即隨著砂土顆粒變粗，其數(shù)值增大，這與表1中α參數(shù)的變化特征類似。氣候情景對參數(shù)β的影響相對較小。參數(shù)k為4.8～10.7，總體上隨著氣候情景向干旱化方向發(fā)展而增大，但隨著砂土顆粒變粗并不顯示出單調(diào)的增減變化，似乎在中砂的情景中其數(shù)值最小。

4 成果應(yīng)用與討論

本文模擬研究得到的結(jié)果集中體現(xiàn)在式(11)和表3中，可以用于對巴丹吉林沙漠的潛水蒸發(fā)量進(jìn)行重新評估。為了便于對比，本文計(jì)算潛水蒸發(fā)量的空間范圍與陳添斐[17]所取的范圍相同，屬于30 km×23.6 km的矩形區(qū)域。該區(qū)域發(fā)育了15個(gè)大小不同的湖泊，湖泊總面積約為8.9 km2。根據(jù)已知的200 m分辨率數(shù)字高程數(shù)據(jù)和地下水位柵格數(shù)據(jù)，可以計(jì)算得到地下水埋深的柵格數(shù)據(jù)。然后，考慮不同的情景，采用式(11)計(jì)算每個(gè)非湖泊柵格的多年平均潛水蒸發(fā)量，累加得到評價(jià)區(qū)的潛水蒸發(fā)總量，與湖面蒸發(fā)量進(jìn)行對比。

在干旱氣候(C3)和細(xì)砂(S1)情景下，計(jì)算得到評價(jià)區(qū)的潛水蒸發(fā)總量為2.9×107m3/a，湖面蒸發(fā)總量為1.1×107m3/a。這意味著評價(jià)區(qū)的潛水蒸發(fā)量比湖面蒸發(fā)量大1.6倍，二者比值接近2.6∶1。在濕潤氣候(C1)和粗砂(S3)情景下，潛水蒸發(fā)總量為2.6×107m3/a，與湖面蒸發(fā)量的比值接近2.5∶1。另外，地下水位動(dòng)態(tài)觀測結(jié)果顯示巴丹吉林沙漠地下水的水位非常穩(wěn)定[39-40]，年動(dòng)態(tài)變幅小于0.3 m，波動(dòng)輕微。在平均水位的基礎(chǔ)上增加0.15 m或減小0.15 m計(jì)算得到的潛水蒸發(fā)量，與采用水位平均值得到結(jié)果相比，相對偏移量均小于3%?？傊u價(jià)區(qū)的潛水蒸發(fā)耗水顯著大于湖面蒸發(fā)耗水。

陳添斐[17]模擬計(jì)算發(fā)現(xiàn)潛水蒸發(fā)總量與湖面蒸發(fā)總量的比值接近1∶1，這是極限埋深取1.5 m且進(jìn)行線性插值得到的結(jié)果。從圖6d可以看出，這樣會低估相同埋深下的Eg值，從而造成計(jì)算結(jié)果偏小?？紤]到實(shí)際極限埋深大于3.0 m以及Eg呈現(xiàn)的非線性變化，本文計(jì)算的結(jié)果更加可靠。巴丹吉林沙漠有100多個(gè)湖泊，總面積接近20 km2，湖面蒸發(fā)總量接近2.3×107m3/a。若按照本文的比例進(jìn)行估計(jì)，則巴丹吉林沙漠無湖洼地的潛水蒸發(fā)總量至少達(dá)到5.8×107m3/a，兩者總和確定的地下水排泄消耗量至少為8.1×107m3/a。以往對巴丹吉林沙漠水均衡的研究容易忽略潛水蒸發(fā)，從水均衡分析上來講顯然是不妥當(dāng)?shù)?，需要重新加以考慮。

本文就一些簡單情景下的潛水蒸發(fā)進(jìn)行了數(shù)值模擬，僅取得了一些初步結(jié)論，還存在著許多問題有待進(jìn)一步研究。首先，目前還缺乏用來校驗(yàn)?zāi)Ｐ偷膶?shí)測潛水蒸發(fā)數(shù)據(jù)。建議下一步工作中對巴丹吉林沙漠的潛水蒸發(fā)進(jìn)行長期監(jiān)測，確定各種典型情景下的實(shí)際潛水蒸發(fā)情況，為進(jìn)一步校準(zhǔn)和完善潛水蒸發(fā)模型提供定量依據(jù)。其次，本模型只考慮了液態(tài)水的運(yùn)移，實(shí)際上，在沙漠包氣帶，汽流和溫度對水分運(yùn)移的作用也很重要[41-42]。特別是在凍結(jié)期，溫度對土壤水的運(yùn)動(dòng)起主導(dǎo)作用[43]。因此，在下一步的研究中，需要進(jìn)一步考慮水、汽、熱耦合遷移以及凍結(jié)期的潛水蒸發(fā)問題。第三，本文模型中采用了給定地下水埋深的方式進(jìn)行模擬，實(shí)際上地下水位與蒸發(fā)存在相互作用并可能造成地下水位的變動(dòng)，從而對模擬結(jié)果產(chǎn)生影響[44]。為了精確地模擬潛水蒸發(fā)量，在下一步的研究中還需要考慮潛水蒸發(fā)與地下水位的耦合效應(yīng)。此外，在沙漠中，湖泊和洼地并非均勻分布。因此，本次研究獲得的潛水蒸發(fā)量和湖面蒸發(fā)量的比值在空間范圍上具有局限性，只能對應(yīng)于本文所選取的計(jì)算范圍。若是將計(jì)算范圍擴(kuò)大到沙漠西北部湖泊較少的區(qū)域，可能該比值會更大，這有待于在未來的工作中加以探討。

5 結(jié)論

(1)隨地下水埋深的加大，多年平均潛水蒸發(fā)量呈現(xiàn)非線性的降低趨勢，極限埋深為3.3～3.8 m，平均極限埋深約3.5 m。

(2)在地下水埋深0.5～1.5 m的區(qū)間，潛水蒸發(fā)量隨地下水埋深的增大而迅速減小。當(dāng)?shù)叵滤裆畹扔?.0 m時(shí)，潛水蒸發(fā)量隨砂土類型的變化也迅速變化。在埋深小于0.5 m和埋深大于1.5 m的區(qū)間，隨著地下水埋深、氣候情景和砂土類型的變化，潛水蒸發(fā)量只發(fā)生輕微變化。

(3)新提出的類VG經(jīng)驗(yàn)公式可以準(zhǔn)確地?cái)M合各情景下的潛水蒸發(fā)曲線，尤其能夠準(zhǔn)確地反映在d/dmax為0.2～0.4時(shí)曲線的急劇變化特征。經(jīng)驗(yàn)公式中的參數(shù)β隨砂土顆粒變粗而增大，對氣候情景的變化并不敏感。經(jīng)驗(yàn)公式中的另一個(gè)參數(shù)k隨著氣候變干而逐漸增大，對砂土特征變化不太敏感。

(4)巴丹吉林沙漠湖泊集中區(qū)的潛水蒸發(fā)總量明顯大于湖水蒸發(fā)總量，在沙漠水均衡分析中不能忽略潛水蒸發(fā)的作用。