1、王旭升中國(guó)地質(zhì)大學(xué)(北京)目 錄1. 如何獲取HYDRUS-1D22. 版權(quán)聲明23. 參考資料24. HYDRUS-1D的WINDOWS界面25. 設(shè)計(jì)模型36. 使用HYDRUS-1D創(chuàng)建模型47. 輸入模型控制信息58. 水流模型迭代計(jì)算參數(shù)69. 水流模型土壤水力特性模型710. 水流模型土壤水分特征曲線711. 水流模型邊界條件812. 水流模型定水頭或通量邊界設(shè)置913. 根系吸水吸水模型914. 根系吸水水分脅迫參數(shù)1015. 輸入可變邊界條件的信息1016. 編輯土壤剖面使用圖形界面1117. 編輯土壤剖面使用表格1318. 運(yùn)行模型1419. 察看結(jié)果1420. 輸出結(jié)果14

2、HYDRUS-1D是一個(gè)共享專業(yè)軟件，用于模擬一維變飽和度地下水流、根系吸水、溶質(zhì)運(yùn)移和熱運(yùn)移。本手冊(cè)只介紹應(yīng)用HYDRUS1D模擬垂向剖面水流和根系吸水的操作方法。1. 如何獲取HYDRUS-1DHYDRUS-1D由位于歐盟捷克的PC-Progress工程軟件開(kāi)發(fā)公司發(fā)行，用戶可以登錄該公司首頁(yè): 。為了下載HYDRUS-1D，應(yīng)先注冊(cè)成為用戶，然后下載Hydrus-1D的安裝文件：。這個(gè)文件對(duì)應(yīng)目前HYDRUS-1D的最高版本。2. 版權(quán)聲明HYDRUS-1D的作者為:(1) J. Simunek, Department of Environmental Sciences, Univers

3、ity of California Riverside, Riverside, California, USA.(2) M. Sejna, PC Progress, Prague, Czech Republic.(3) M.Th. van Genuchten, Department of Mechanical Engineering, Federal University of Rio de Janeiro, Rio de Janeiro, Brazil.感謝他們提供了一個(gè)如此精美而又免費(fèi)使用的專業(yè)軟件，幫助我們從事有關(guān)的科學(xué)和教育工作。當(dāng)你運(yùn)行解壓文件后，會(huì)在您的電腦中產(chǎn)生一個(gè)安裝目錄，其中

4、包含Setup.exe可執(zhí)行文件。運(yùn)行這個(gè)文件即可安裝HYDRUS-1D軟件。當(dāng)您安裝HYDRUS-1D時(shí)，象安裝其它軟件一樣，會(huì)出現(xiàn)一個(gè)許可協(xié)議，從中可知本共享軟件也受到美國(guó)法規(guī)的保護(hù)。3. 參考資料 HYDRUS-1D安裝之后，在軟件運(yùn)行目錄下有HYDRS-1D Manual.pdf文件。從這個(gè)文件您可以了解到HYDRUS-1D的一些技術(shù)細(xì)節(jié)，如水流、溶質(zhì)運(yùn)移、熱流的方程、一些處理專門(mén)問(wèn)題的模型、輸入輸出文件等等。有一個(gè)Examples目錄，包含大量的模擬算例可供參考。用戶還可以參考以下文獻(xiàn)：· imnek, J., M. Th. van Genuchten, and M. ej

5、na, Development and applications of the HYDRUS and STANMOD software packages, and related codes, Vadose Zone Journal, doi:10.2136/VZJ2007.0077, Special Issue ”Vadose Zone Modeling”, 7(2), 587-600, 2008. · Jacques, D., J. imnek, D. Mallants, and M. Th. van Genuchten, Modeling coupled hydrologica

6、l and chemical processes: Long-term uranium transport following mineral phosphorus fertilization, Vadose Zone Journal, doi:10.2136/VZJ2007.0084, Special Issue ”Vadose Zone Modeling”, 7(2), 698-711, 2008.· 0074, Special Issue ”Vadose Zone Modeling”, 7(2), 782-797, 2008. 這些文獻(xiàn)都可以從 下載。4. HYDRUS-1D的

7、WINDOWS界面運(yùn)行HYDRUS-1D，可以看到一個(gè)Windows的界面如下：前處理工具后處理工具模擬計(jì)算圖 1 所有的前后處理在界面中一目了然，左邊是前處理工具，右邊是后處理工具。其中前處理的各項(xiàng)功能如下圖所示。模擬內(nèi)容選項(xiàng)幾何形狀參數(shù)及剖面方式時(shí)間信息輸出方式水流 迭代求解控制參數(shù)水流 土壤水分特征模型水流 土壤水分特征曲線參數(shù)水流 邊界條件溶質(zhì)運(yùn)移 一般信息溶質(zhì)運(yùn)移 運(yùn)移參數(shù)溶質(zhì)運(yùn)移 化學(xué)反應(yīng)參數(shù)溶質(zhì)運(yùn)移 邊界條件根系吸水 模型根系吸水 水分脅迫函數(shù)可隨時(shí)間變化的邊界條件土壤剖面 圖形界面土壤剖面 數(shù)據(jù)列表圖 25. 設(shè)計(jì)模型在使用HYDRUS-1D之前，您需要對(duì)飽和-非飽和水流模擬的

8、基本原理有所了解，并設(shè)計(jì)出自己想做的模型，準(zhǔn)備好數(shù)據(jù)。一個(gè)剖面水流模型通常包含以下幾個(gè)要素：(1) 土壤剖面從地面算起的深度，準(zhǔn)備模擬那個(gè)時(shí)間段的水分變化。(2) 土壤分幾層，每層土壤的滲透性參數(shù)和水分特征曲線是怎樣的。(3) 根系是怎么分布的。(4) 是否已經(jīng)確定地面降雨入滲、蒸發(fā)蒸騰的信息，特別是它們隨時(shí)間的變化。(5) 是否已經(jīng)確定剖面底部的狀態(tài)屬于哪種類型的邊界條件。下面是一個(gè)參考模型的設(shè)計(jì)圖:時(shí)間E0降水量根系層30 d初始潛水面隔水底板1 m1.5 m1.5 m1 m圖 36. 使用HYDRUS-1D創(chuàng)建模型打開(kāi)HYDRUS-1D軟件，選擇”File/new”菜單，新建一個(gè)模型。在

9、name一欄中輸入本模型的名稱”test”，更改模型存放的目錄。圖 4需要注意的是，HYDRUS-1D模型本身在計(jì)算機(jī)中就表現(xiàn)為一系列的輸入輸出文件，它們存放在與模型名稱一致的目錄中。本例中，軟件會(huì)自動(dòng)創(chuàng)建一個(gè)名稱為”test”的目錄，而”CATOOLSHYDR1DProjects”中除了test目錄之外，還有一個(gè)test.h1d文件。這是一個(gè)模型項(xiàng)目(project)文件，告訴軟件下次到哪里去尋找模型。 模型創(chuàng)建之后，會(huì)顯示前處理和后處理窗口(圖5)。由于是新模型，還沒(méi)有任何模擬結(jié)果，所以后處理窗口是空白的。圖 57. 輸入模型控制信息首先，在前處理窗口雙擊Main Processes，在彈

10、出的對(duì)話框中輸入模型的描述: a test model. 然后在Simulate 一欄中選中Root Water Uptake，表示想處理根系吸水問(wèn)題。電擊OK之后，前處理窗口將增加處理根系吸水的工具條。圖 6下一步，是輸入模型的幾何信息和土層劃分信息。在前處理窗口雙擊Geometry Information，在彈出的對(duì)話框中輸入如圖7所示的數(shù)據(jù)。模型有2個(gè)土層長(zhǎng)度單位是cm土壤剖面的深度是300 cm圖 7接下來(lái)輸入時(shí)間信息，在前處理窗口雙擊Time Information，會(huì)彈出一個(gè)對(duì)話框（圖8）。時(shí)間單位時(shí)間步長(zhǎng)信息模擬30 d內(nèi)的變化上邊界是隨時(shí)間變化的，每天一組數(shù)據(jù)，共30組數(shù)據(jù)自動(dòng)

11、處理蒸騰量在每天24小時(shí)內(nèi)的變化圖 8這個(gè)對(duì)話框中提供了一些靈活的選項(xiàng)來(lái)處理上邊界條件的變化，下面簡(jiǎn)要加以說(shuō)明：(1) 蒸騰量的每日周期變化 HYDRUS-1D 可以使用一個(gè)經(jīng)驗(yàn)公式來(lái)處理每天24小時(shí)潛在蒸騰量的變化，設(shè)某天的潛在蒸騰量為（例如用Pemman公式獲取的, cm/d），則其中Tp(t)是瞬時(shí)潛在蒸騰量，t為時(shí)間。模型假設(shè)早上6點(diǎn)之前以及晚上18點(diǎn)-24點(diǎn)的蒸騰量總和只占全天蒸騰量的1%。注意本例中蒸騰量的單位是 cm/d。 (2) 降水量的周期變化如果在你的模型中降水量是周期性變化的，HYDRUS-1D也可以用一個(gè)公式來(lái)處理其中是周期Dt內(nèi)的平均降雨量。(3) 使用氣象數(shù)據(jù)也可以

12、在HYDRUS-1D中輸入氣象數(shù)據(jù)，它將自動(dòng)利用這些數(shù)據(jù)計(jì)算潛在蒸散量ETp?？梢赃x擇FAO組織推薦的Penman-Monteith公式，也可以選擇Hargreaves公式。這些公式需要輻射、氣溫、濕度之類的氣象數(shù)據(jù)。模型的另一個(gè)控制信息是對(duì)模擬結(jié)果的輸出如何進(jìn)行設(shè)置。在前處理窗口雙擊Print Information工具條，彈出一個(gè)對(duì)話框。本例中確定輸出30組模擬結(jié)果，每天輸出1組。輸出30組結(jié)果，每天1組可以確定需要輸出哪些時(shí)間點(diǎn)的結(jié)果圖 98. 水流模型迭代計(jì)算參數(shù)HYDRUS-1D是采用迭代法來(lái)處理非線性Richards方程的。在前處理窗口雙擊Water Flow- Iteration

13、 Criteria工具條，彈出一個(gè)設(shè)置迭代參數(shù)的對(duì)話框（圖10）。迭代控制參數(shù)的設(shè)置具有高度的專業(yè)技術(shù)性，除非特別了解，一般可以使用默認(rèn)值。如果模擬結(jié)果出現(xiàn)不收斂的情況，需要對(duì)最大迭代次數(shù)、迭代精度等參數(shù)進(jìn)行調(diào)整，但是在缺乏經(jīng)驗(yàn)的情況下很難操作。最多迭代次數(shù)含水量迭代精度壓力水頭迭代精度增大步長(zhǎng)迭代次數(shù)信號(hào)縮小步長(zhǎng)的迭代次數(shù)信號(hào)步長(zhǎng)增大比例步長(zhǎng)縮小比例最小吸力間距最大吸力間距用于生成水分特征曲線的散點(diǎn)值圖 10HYDRUS-1D 采用自動(dòng)控制時(shí)間步長(zhǎng)的方法來(lái)處理迭代的收斂性。對(duì)于每個(gè)時(shí)步，如果迭代次數(shù)太多，就縮小時(shí)間步長(zhǎng)；如果沒(méi)經(jīng)過(guò)幾次迭代就達(dá)到收斂精度，則適當(dāng)增大時(shí)間步長(zhǎng)。9. 水流模型土壤

14、水力特性模型水分特征曲線是非飽和土壤的重要物理性質(zhì)，HYDRUS-1D提供了幾種方法來(lái)處理與之有關(guān)的參數(shù)。在前處理窗口雙擊Water Flow- Soil Hydraulic Properties工具條，彈出一個(gè)設(shè)置水力特性模型的對(duì)話框（圖11）。單孔介質(zhì)模型雙重介質(zhì)模型雙重孔隙度 雙重滲透性吸濕和疏干滯后過(guò)程模型主要在溶質(zhì)運(yùn)移問(wèn)題中使用圖 11在一般情況下，選擇單孔介質(zhì)模型，并選擇用van Genuchten-Mualem公式處理土壤的水力特性就可以了。如果還要模擬溶質(zhì)運(yùn)移，可能需要考慮雙重介質(zhì)模型。雙重介質(zhì)在同一個(gè)點(diǎn)有兩個(gè)孔隙度或兩個(gè)滲透率，相當(dāng)于兩種介質(zhì)的混雜。雙重介質(zhì)模型能夠模擬這兩種

15、“介質(zhì)”之間的水分和鹽分交換。10. 水流模型土壤水分特征曲線在前處理窗口雙擊Water Flow- Soil Hydraulic Parameters工具條，彈出一個(gè)設(shè)置水分特征曲線參數(shù)的對(duì)話框（圖12）。本例中選擇van Genuchten-Mualem公式處理水分特征曲線，其中a, n, l均為控制因子。HYDRUS-1D軟件中提供了一組土壤經(jīng)驗(yàn)參數(shù)庫(kù)，可供用戶參考。本例中兩層土壤的參數(shù)直接從數(shù)據(jù)庫(kù)中調(diào)出：第1層對(duì)應(yīng)Sandy loam，第二層對(duì)應(yīng)sand。編號(hào) qr qs a (cm-1) n Ks (cm/d) l土壤數(shù)據(jù)庫(kù)細(xì)砂壤土中砂圖 12在輸入?yún)?shù)時(shí)，請(qǐng)注意參數(shù)的單位。11.

16、水流模型邊界條件在前處理窗口雙擊Water Flow- Boundary Conditions工具條，彈出一個(gè)設(shè)置邊界條件的對(duì)話框（圖13）。地面邊界類型定壓力水頭定水分通量大氣邊界, 可積水大氣邊界,產(chǎn)流變水頭變水頭和通量初始條件使用壓力水頭使用含水量下端邊界類型定壓力水頭定水分通量變地下水位變流量自由下滲排水滲出面水平排水大氣邊界基于ETp和LAI消光系數(shù)最大積水厚度圖 13上邊界條件有6種類型，下邊界條件有8種類型。邊界類型的確定需要考慮實(shí)際條件，在本算例中，上邊界選擇大氣邊界條件，在降雨量很大時(shí)地表可以產(chǎn)生積水。植被蒸騰量和土壤蒸發(fā)量分開(kāi)處理，HYDRUS-1D推薦使用一個(gè)經(jīng)驗(yàn)公式來(lái)把

17、潛在蒸散量分割為蒸騰潛力和土壤蒸發(fā)潛力：其中ETp為潛在蒸散量（可以使用Penman-Monteith公式處理氣象數(shù)據(jù)得到, cm/d），Tp為潛在蒸騰量(cm/d)，Ep為土壤潛在蒸發(fā)量(cm/d)，LAI是葉面積指數(shù)，k為消光系數(shù)，取決于太陽(yáng)角度、植被類型及葉片空間分布特征。SCF是一個(gè)中間參數(shù)，即土壤覆蓋度(Soil cover fraction)。在闊葉植被發(fā)育的情況下，消光系數(shù)的經(jīng)驗(yàn)值為k。12. 水流模型定水頭或通量邊界設(shè)置如果邊界條件中包含定水頭或定通量的邊界，則在前處理窗口雙擊Water Flow- ConstantBC工具條，彈出一個(gè)設(shè)置邊界數(shù)據(jù)的對(duì)話框。本算例模型中，下邊界

18、為定流量邊界，實(shí)際上就是隔水邊界，因此直接輸入0即可。13. 根系吸水吸水模型在前處理窗口雙擊Root Water Uptake- Models工具條，彈出一個(gè)處理根系吸水模型的對(duì)話框（圖14）。水分脅迫模型鹽分脅迫模型根系吸鹽模型補(bǔ)償吸水域值（臨界濕潤(rùn)度）圖 14HYDRUS-1D使用水分脅迫和鹽分脅迫模型處理根系的吸水。對(duì)于水分脅迫模型，計(jì)算公式為其中Tp是潛在蒸騰量(cm/d)，Ta是實(shí)際蒸騰量(cm/d)，S(x)是吸水強(qiáng)度函數(shù)(cm/(cm.d)，注意x坐標(biāo)實(shí)際表示深度)，a(h)是水分脅迫函數(shù)，h為土壤壓力水頭(cm)，b(x)是根系吸水分配（密度）函數(shù)，LR為根系層的深度。水分脅

19、迫函數(shù)有2種經(jīng)驗(yàn)表示方法，即Feddes模型和S-Shape模型。Feddes模型是一個(gè)梯形函數(shù)，只需要知道h值。而S-Shape模型把水分脅迫和葉片氣孔的壓力水頭聯(lián)系起來(lái)，需要知道氣孔壓力水頭的數(shù)值hf。土壤的濕潤(rùn)度可以表示為 (2.22)但是如果直接用這種方法來(lái)計(jì)算實(shí)際蒸騰量有一定的問(wèn)題。植被其實(shí)可以調(diào)節(jié)不同深度的水分脅迫響應(yīng)特征；某個(gè)深度土壤干燥吸不上水，植被可以加大在比較濕潤(rùn)的土層的吸水量，以補(bǔ)償不足。這種現(xiàn)象稱為補(bǔ)償吸水。為了模擬根系補(bǔ)償吸水，HYDRUS-1D提供了一種簡(jiǎn)化的模型，即如果濕潤(rùn)度高于某個(gè)臨界值(w>wc)，植被根系可以通過(guò)補(bǔ)償機(jī)制充分吸水達(dá)到潛在蒸騰量。如果濕潤(rùn)

20、度低于這個(gè)臨界值，補(bǔ)償機(jī)制受到抑制，發(fā)生整體的水分脅迫，根系吸水總量將低于潛在蒸騰量，并正比于濕潤(rùn)度。如果不考慮這種補(bǔ)償吸水機(jī)制，可以令wc =1。14. 根系吸水水分脅迫參數(shù)在前處理窗口雙擊Root Water Uptake- Water Stress Reduction工具條，彈出一個(gè)處理水分脅迫參數(shù)的對(duì)話框（圖15）。本算例中直接從數(shù)據(jù)庫(kù)中調(diào)入Wheat的經(jīng)驗(yàn)值。P3 P2L P2H Popt P0壓力水頭(pressure head, cm)a=1a=0r2Hr2L數(shù)據(jù)庫(kù)吸水強(qiáng)度(cm/d)脅迫函數(shù)兩個(gè)臨界吸水強(qiáng)度圖 15 根系水分脅迫Feddes模型參數(shù)15. 輸入可變邊界條件的信息

21、在前處理窗口雙擊Variable Boundary Conditions工具條，彈出一個(gè)處理時(shí)間序列數(shù)據(jù)的對(duì)話框（圖16）。本算例中，在步驟(11)中已經(jīng)把地面處理大氣邊界，同時(shí)又選擇使用消光系數(shù)法劃分植被蒸騰和土面蒸發(fā)，因此需要輸入每天的降水、潛在蒸散量、葉面積指數(shù)等數(shù)據(jù)。還有一個(gè)需要輸入的數(shù)據(jù)是最小壓力水頭值，即地面土壤達(dá)到最干燥狀態(tài)時(shí)的壓力水頭。從理論上講，當(dāng)土壤十分干燥時(shí)，吸力很大，而液態(tài)孔隙水的壓強(qiáng)很小，與空氣濕度保持平衡關(guān)系，因此有其中hA為最小壓力水頭，Hr為空氣絕對(duì)濕度，RT/Mg為空氣的摩爾氣體常數(shù)。空氣濕度雖然可以通過(guò)氣象數(shù)據(jù)得到，但這里公式需要的是近地面的空氣濕度。一般情

22、況下，取飽和水汽濕度是可取的，因?yàn)? cm深度以下土壤空氣的濕度往往都是飽和的，只不過(guò)隨溫度發(fā)生變化。因此，可以根據(jù)近地面氣溫的變化來(lái)推算地表土壤的空氣濕度（飽和水汽濕度），再換算成壓力水頭。HYDRUS-1D中需要輸入的是最小壓力水頭的絕對(duì)值，缺省值為hCritA=|hA|=106 cm=104 m這個(gè)數(shù)值只會(huì)對(duì)土壤蒸發(fā)起作用。HYDRUS-1D建議：hCritA 所對(duì)應(yīng)的土壤含水量應(yīng)該至少比殘余含水量大0.005，在模擬根系吸水的情況下，hA還應(yīng)該低于圖15中的P3。否則（hA>P3），當(dāng)根系吸水的臨界值壓力水頭(P3)和地面蒸發(fā)的最小壓力水頭 (hA) 滿足時(shí),會(huì)導(dǎo)致回流(infl

23、ow)現(xiàn)象，這是不合理的。除非存在特別干燥的情況，模型一般不需要仔細(xì)處理這些問(wèn)題。圖 16表1時(shí)間(d)降水量(cm/d)ETp (cm/d)hCritA (cm)LAI時(shí)間(d)降水量(cm/d)ETp (cm/d)hCritA (cm)LAI101000001601000002010000017010000030100000180100000401000001901000005010000020010000060100000210100000701000002201000008100000230100000901000002410000010010000025100000110100000

24、261000001201000002710000013310000028010000014010000029010000015100000300100000數(shù)據(jù)可以先在Excel中準(zhǔn)備好，如表1。這些數(shù)據(jù)可以拷貝到圖16的電子表格中。這些數(shù)據(jù)顯示的降水量、蒸散潛力和葉面積指數(shù)變化特征如圖16右圖所示。在第17日由于莊稼收割，葉面及指數(shù)大幅度下降。16. 編輯土壤剖面使用圖形界面在前處理窗口雙擊Soil Profile- Graphical Editor工具條，程序?qū)棾鲆粋€(gè)處理土壤剖面的軟件（圖17）。這個(gè)圖形軟件的使用比較簡(jiǎn)單，我們需要注意的是在Conditions菜單下面有很多子菜單，包含

25、處理各種問(wèn)題的功能菜單。首先要做的事情，是確定把土壤剖面離散化為多少個(gè)節(jié)點(diǎn)。本算例土壤模型深度為3 m，我們希望節(jié)點(diǎn)間距達(dá)到 1 cm，因此需要301個(gè)節(jié)點(diǎn)。選擇菜單Conditions/ Profile Discretization，在下拉工具條中把Number修改為301。土壤剖面分成31個(gè)節(jié)點(diǎn)，則節(jié)點(diǎn)間距為 10 cmConditionsProfile DiscretizationMaterial DistributionRoot DistributionScaling factorInitial ConditionsSubregionsObservation PointsConditi

26、ons菜單剖分節(jié)點(diǎn)土壤巖性分層根系分布尺度因子初始條件子區(qū)分布觀察點(diǎn)修改節(jié)點(diǎn)數(shù)圖 17 接下來(lái)，確定土層的分布，本模型有2個(gè)土層。缺省的土層編號(hào)為index=1，就是圖12中的細(xì)砂壤土層。選擇菜單Conditions/Material Distribution，在下拉工具條中使用Edit condition，把下部土壤層設(shè)置index=2，這個(gè)土層編號(hào)為2，實(shí)際上就是圖12中的中砂層。設(shè)置2層土壤巖性編號(hào)1. 點(diǎn)擊工具條2.在剖面底部點(diǎn)擊鼠標(biāo)左鍵再放開(kāi)，從下往上移動(dòng)鼠標(biāo)，選定土層范圍，然后再點(diǎn)擊鼠標(biāo)左鍵，修改土層編號(hào)修改節(jié)點(diǎn)數(shù)右下角是鼠標(biāo)所在的z坐標(biāo)高度ConditionsMaterial D

27、istribution圖 18再就是確定根系隨深度的分布，根系分布函數(shù)b(x)是一個(gè)很特殊的函數(shù)，它滿足以下條件其中LR為根系層厚度, x是深度；Dz是節(jié)點(diǎn)間距，bn是每個(gè)間距中的根系分布函數(shù)值，M是根系層占節(jié)點(diǎn)數(shù)。本算例中根系層的厚度為1 m，假設(shè)根系分布函數(shù)為線性，并有, x£1容易證明上式滿足積分為1的條件。設(shè)置方法為，選擇菜單Conditions/Root Distribution，在在下拉工具條中使用Edit condition，劃定根系層范圍(100 cm)，把頂部數(shù)值設(shè)置為0，底部數(shù)值設(shè)置為2，讓程序自動(dòng)進(jìn)行線性差值。根系分布函數(shù)使用工具條ConditionsRoot Distribution圖 19下一步處理初始條件，選擇菜單Conditions/InitialCondition。本算例中初始地下水位高于底板100 cm，假設(shè)土壤剖面初始狀態(tài)是靜力平衡態(tài)，則模型底部的壓力水頭為+100 cm，地面的壓力水頭為-200 cm。使用Edit Condition工具條，分別設(shè)置頂部和底部的壓力水頭，并讓程序自動(dòng)插值