蒸發(fā)條件下寧夏銀北地區(qū)鹽漬土水分運(yùn)移數(shù)值模擬

趙文娟+李春光

摘要：對(duì)鹽漬土中的水鹽動(dòng)態(tài)進(jìn)行定量分析是提出適宜于寧夏銀北地區(qū)鹽漬土調(diào)控措施的基礎(chǔ)工作。蒸發(fā)條件下土壤水分運(yùn)動(dòng)模擬初步研究表明，將有限體積法應(yīng)用到求解蒸發(fā)條件下的土壤水分運(yùn)動(dòng)模擬計(jì)算中，模擬值與實(shí)測(cè)值吻合較好；春季土壤含水率模擬曲線在縱斷面較為一致，地表蒸發(fā)量對(duì)0.2～0.6 m的土壤含水率有較大的影響。

關(guān)鍵詞：土壤水分；蒸發(fā)；數(shù)值模擬；鹽漬土

中圖分類(lèi)號(hào)：S152.7 文獻(xiàn)標(biāo)識(shí)碼：A 文章編號(hào)：0439-8114（2015）15-3638-04

DOI：10.14088/j.cnki.issn0439-8114.2015.15.013

Abstract： Its the basis of the work to use the numerical simulation method for finding appropriate control measures in Yinbei region of Ningxia province. The results showed that there was a good agreement between simulated and measured values when the finite volume method was selected to simulate the movement of soil moisture evaporation. The simulated curve of soil moisture in the horizontal section was consistent with that of the real values. Surface evaporation had a greater impact on the soil moisture of 0.2～0.6 m.

Key words： soil moisture； evaporation； numerical simulation； saline soil

寧夏回族自治區(qū)屬于黃河中上游半干旱-半漠境鹽漬區(qū)，鹽漬化土地面積為20 718 km2，占自治區(qū)總面積的40%。其中尤以銀北地區(qū)分布最廣，鹽漬化面積為17 070 km2，約占全區(qū)鹽漬化土地總面積的81%左右，該地區(qū)含鹽量0.6%～1.0%的重度鹽漬土就占到本區(qū)鹽漬化土壤面積25.66%[1]。由于銀北地區(qū)地形低洼，排水不暢，地下水位高致使水中的鹽分在強(qiáng)烈的蒸發(fā)作用下堆積于土壤內(nèi)，這是形成土壤鹽漬化的主要原因。通過(guò)對(duì)該地區(qū)的鹽化土地實(shí)施水利工程改造和農(nóng)業(yè)生物改良調(diào)控措施，鹽漬化得到了整體控制[2]。但近幾年來(lái)該地區(qū)局部的鹽漬化現(xiàn)象明顯加重，主要是由于打淺井使用高礦化度地下水進(jìn)行灌溉及淋洗，加劇了土壤的鹽害[3]。這就需要提出更為適宜銀北地區(qū)土壤鹽漬化調(diào)控措施，作為該地區(qū)農(nóng)業(yè)可持續(xù)發(fā)展和生態(tài)環(huán)境保護(hù)的基礎(chǔ)。

本研究是通過(guò)建立能夠描述寧夏銀北地區(qū)鹽漬土中水鹽運(yùn)移的數(shù)學(xué)模型，采用數(shù)值計(jì)算方法進(jìn)行模擬計(jì)算，對(duì)鹽漬土中的水鹽動(dòng)態(tài)進(jìn)行定量分析。分析模擬結(jié)果，掌握水分、鹽分在鹽漬土體中的基本運(yùn)動(dòng)規(guī)律，為確定合理的土壤鹽漬化調(diào)控模式提供科學(xué)依據(jù)。

1 材料與方法

1.1 土壤水分運(yùn)動(dòng)數(shù)學(xué)模型

土壤水分運(yùn)移模型是基于對(duì)流擴(kuò)散方程和初邊值條件共同構(gòu)建的模型，模型中的各類(lèi)土壤水動(dòng)力參數(shù)均需要通過(guò)試驗(yàn)來(lái)確定。該模型屬于非線性偏微分方程，當(dāng)前最有效的求解方法是數(shù)值計(jì)算方法。有限差體積法[4]是當(dāng)前求解流動(dòng)問(wèn)題數(shù)值計(jì)算方法中最成功的方法，該方法是將方程中的每一項(xiàng)在所劃分的控制體積進(jìn)行積分，各項(xiàng)積分均有明確的物理意義，構(gòu)成線性方程組進(jìn)行求解計(jì)算。

1.2.1 有限體積法離散過(guò)程 采用有限體積方法對(duì)二維非飽和土壤水分運(yùn)動(dòng)方程進(jìn)行線性離散。

其中，？準(zhǔn)0為初值含水率（cm3/cm3），D（？準(zhǔn)）為非飽和土壤擴(kuò)散率（cm2/min），為非飽和土壤導(dǎo)水率（cm/min），？著為水量交換強(qiáng)度（cm/d）。

對(duì)方程在圖1所示的控制體內(nèi)積分可得：

由于全隱式格式可以在較大時(shí)間步長(zhǎng)保持計(jì)算結(jié)果的穩(wěn)定性[5]，所以采用全隱格式對(duì)擴(kuò)散項(xiàng)進(jìn)行時(shí)間上的積分，設(shè)？駐x=？駐z=？駐d則？駐V=A·？駐x=A·？駐z，得：

1.2.2 邊界條件處理 上邊界條件為地表，與當(dāng)?shù)氐臍庀笠蛩叵嚓P(guān)聯(lián)，即與垂直向的水量交換強(qiáng)度？著（cm/d）有關(guān)。當(dāng)？著為負(fù)數(shù)時(shí)為該時(shí)期的蒸發(fā)強(qiáng)度。根據(jù)監(jiān)測(cè)蒸發(fā)量數(shù)據(jù)即可得到整個(gè)監(jiān)測(cè)區(qū)域內(nèi)的水量交換強(qiáng)度值。

采用上述相同方法給出邊界條件在計(jì)算區(qū)域內(nèi)離散得到各個(gè)端點(diǎn)的離散系數(shù)，見(jiàn)表1。以上系數(shù)通過(guò)編寫(xiě)MATLAB程序進(jìn)行迭代運(yùn)算，即可得到蒸發(fā)條件下的土壤水分運(yùn)動(dòng)曲線。運(yùn)行前，需給出初始值。假設(shè)土壤的含水率隨監(jiān)測(cè)深度的不同而變化，呈現(xiàn)為隨深度變化的曲面，即有：

？準(zhǔn)0（SD，0）=HSL （6）

其中，SD為土壤不同深度（m），以向量表示；HSL為不同深度的土壤含水率（cm3/cm3），以向量表示。

1.2 試驗(yàn)區(qū)情況

試驗(yàn)區(qū)位于寧夏銀川市賀蘭縣南梁臺(tái)子鐵西村。該村地理坐標(biāo)為東經(jīng)105°58，北緯38°41，位于賀蘭縣的西北方向。試驗(yàn)區(qū)為一面積約459 m2的條形區(qū)域。距離地表20 cm試驗(yàn)土壤的平均初始含水率為15.843%，土壤容重為1.6 g/cm3。

2 結(jié)果與分析

利用時(shí)間延遲反射儀對(duì)試驗(yàn)區(qū)的15個(gè)以三行五列形式布置的測(cè)點(diǎn)進(jìn)行土壤含水率和電導(dǎo)率分時(shí)段、分層的測(cè)定工作。其中相關(guān)的水動(dòng)力學(xué)參數(shù)有：

選擇2012年5月16日10點(diǎn)到5月17日10點(diǎn)的時(shí)段對(duì)自西向東的由2號(hào)、5號(hào)、8號(hào)、11號(hào)和14號(hào)、17號(hào)、20號(hào)及23號(hào)組成的縱斷面數(shù)值模擬土壤水的運(yùn)動(dòng)過(guò)程，其他測(cè)點(diǎn)作為邊界條件進(jìn)行模擬計(jì)算。當(dāng)日蒸發(fā)所引起的土壤水通量為6.4 cm/d，蒸發(fā)強(qiáng)烈。模擬結(jié)果表明，同一深度下的土壤含水率自西向東運(yùn)動(dòng)趨勢(shì)基本保持一致，但模擬數(shù)值自西向東逐漸增加，土壤含水率最低點(diǎn)均出現(xiàn)在0.2 m附近，最大值出現(xiàn)在1.0 m附近。模擬值與實(shí)測(cè)值吻合較好，模擬曲線可以基本表達(dá)土壤水分在鹽漬土內(nèi)的運(yùn)動(dòng)趨勢(shì)，但個(gè)別測(cè)點(diǎn)的實(shí)測(cè)值與模擬值的偏移量較大，可能的原因是假設(shè)模擬區(qū)域的土壤質(zhì)地為均質(zhì)土壤，即土壤容重是相同，在模擬過(guò)程出現(xiàn)了偏移現(xiàn)象，另一方面也可能是測(cè)量誤差造成實(shí)測(cè)值與模擬曲線的偏移。

將5月17日的模擬值作為初始值，模擬土壤水分在5 d后的變化過(guò)程，擬合曲線見(jiàn)圖3。由于氣象數(shù)據(jù)的變化，土壤含水率也發(fā)生變化。含水率的最小值出現(xiàn)在了0.4 m附近，土壤監(jiān)測(cè)斷面含水率與5月17日的相比是增加的，且測(cè)定點(diǎn)的模擬數(shù)值自西向東的增量更為明顯。在0.8 m附近土壤含水率逐漸進(jìn)入穩(wěn)定狀態(tài)，且隨著迭代時(shí)間的增長(zhǎng)，數(shù)值解與實(shí)測(cè)值吻合較好。說(shuō)明運(yùn)用有限體積法計(jì)算土壤縱斷面含水率是有效的，模擬曲線也顯示出0.2～0.6 m處土壤含水率較低。

3 結(jié)論

采用數(shù)值模擬的方法研究蒸發(fā)條件下的寧夏銀北鹽漬土區(qū)土壤水分運(yùn)動(dòng)，在整個(gè)模擬過(guò)程中，模擬值與實(shí)測(cè)值吻合較好，將有限體積法應(yīng)用到求解土壤水分運(yùn)動(dòng)方程是可行的。模擬結(jié)果表明，土壤的含水率在春季的變化較為頻繁，在同一深度的土壤含水率從數(shù)值上呈現(xiàn)西低東高的現(xiàn)象；縱斷面上同一時(shí)刻的土壤含水率模擬曲線形態(tài)是基本一致的，土壤蒸發(fā)量使0.2～0.6 m的土壤含水率較小。

參考文獻(xiàn)：

[1] 何文壽，劉陽(yáng)春，何進(jìn)宇.寧夏不同類(lèi)型鹽漬化土壤水溶鹽含量與其電導(dǎo)率的關(guān)系[J].干旱地區(qū)農(nóng)業(yè)研究，2010，28（1）：111-113.

[2] 俞仁培.對(duì)鹽漬土資源開(kāi)發(fā)利用的思考[J].土壤通報(bào)，2001，32S1：139-140.

[3] 史曉杰，萬(wàn) 力，張永庭，等.銀北地區(qū)土壤鹽漬化形成機(jī)理與模擬研究[J].水文地質(zhì)工程地質(zhì)，2006，34（6）：116-117

[4] VERSTEEG H K，MALALASEKERA W.An introduction to computational fluid dynamics：The Finite V Method[M].London：Longman Group Ltd，1995.

[5] 姜啟源，謝金星，葉 俊.數(shù)學(xué)模型[M].北京：高等教育出版社，2011.