[發明專利]一種片麻巖土石山區山坡尺度水文過程模擬方法

有效

• 申請號: 201710095169.3
• 申請日: 2017-02-22
• 公開（公告）號: CN106874605B
• 公開（公告）日: 2018-12-28
• 發明（設計）人: 甘永德;賈仰文;劉歡;龔家國;牛存穩;仇亞琴;郝春灃;杜軍凱
• 申請（專利權）人: 中國水利水電科學研究院
• 主分類號: G06F17/50
• 分類號: G06F17/50
• 代理公司: 北京國林貿知識產權代理有限公司 11001
• 代理人: 李富華;李桂玲

1.一種片麻巖土石山區山坡尺度水文過程模擬方法，其特征在于，所述方法的計算過程如下：

山坡計算單元劃分：采用等流時線法將山坡劃分為若干基本計算單元；

計算單元垂直剖面劃分：根據山坡植被、土壤和巖石特性，在基本計算單元內分為4層：植被冠層截留層、地表儲留層、土壤層、基巖層；植被冠層截留層又細分為：高植被截留層和矮植被儲留層、草地層和裸地；土壤層進一步分為均質土壤層、土石二元混合介質層；考慮到片麻巖基巖層內廣泛分布著構造節理，且片麻巖本身具有一定持水性和透水性，進一步將基巖層分為基質流區和優先流區；

計算單元內狀態變量包括：植被冠層截留量、洼地儲留量、枯枝落葉儲留量、土壤含水量；主要參數包括：植被最大截留深、洼地最大儲留深、枯枝落葉干重、土壤導水系數、土壤水分特征曲線、土壤含水量、各層土壤厚度、土石二元混合介質碎石質量比系數、基巖層厚度、坡面糙率；

計算單元水文過程計算：包括水文氣象數據展布、降雨期間產匯流過程計算、非降雨期間產匯流過程計算；

水文氣象數據展布：包括水文氣象過程空間尺度展布和降雨時間降尺度展布：

采用泰森多邊形法和反距離加權平均法進行流域內氣象數據的空間展布，包括降雨、氣溫、風速、空氣濕度、凈輻射，計算公式如下：

式中：D表示待插值點估計值；D_pi表示第pi個參證站點數據；pm表示參證站點個數；λ_pi表示第pi個參證站點數據權重；d_pi表示第pi個參證站點同待插值點的距離；pn表示權重指數；

由于日降雨過程的非穩定性，對日降雨數據進一步進行降尺度展布，具體公式如下：

S＝a·P+b (5)

式中：I為時段t_k內最大降水平均雨強；S表示暴雨參數；t為時間，其中，t_k-1<t≤t_k；t_k為時段，其中，k∈0,1,...,N區間的時間，N為時段數；μ表示暴雨衰減系數；P表示日降雨量；T表示日降雨總歷時；a，b表示參數；

降雨期間產匯流過程計算：降雨期間，土壤蒸散發量較小，予以忽略，計算單元內水文過程主要由降雨→植被截留→入滲產流→匯流過程構成；

植被截留計算：植被截留計算，計算公式如下：

Wr_max＝0.2·Veg·LAI (8)

式中：Veg表示裸地-植被域中植被的面積占計算單元的面積比例；Wr表示植被截留水量；Wr_max表示植被最大截留水量；I為時段t_k內最大降水平均雨強；Rr表示植被冠層流出水量；LAI表示葉面積指數；

洼地儲留計算，計算公式如下：

式中：I’為時段t_k內凈雨強；H_u2為地表儲留；H_umax2為土壤表層最大儲留深；R_u2為土壤表面徑流；f_in為入滲率；

枯枝落葉儲留計算，計算公式如下：

H_umax＝z_maxG_k (11)

式中：G_k為枯枝落葉干重；z_max為枯枝落葉最大持水系數；

土壤及基巖層水分運動過程計算：

片麻巖分布區土壤和基巖均具有導水性和持水性，所有土壤層水分運動過程均采用理查茲公式計算：

式中：h為土壤水吸力；C為容水度；SS為源匯項；z為坐標軸；K(h)為導水系數；t為時間；考慮到土壤和基巖層均廣泛分布著裂隙優先流，土壤水分運動過程中將計算單元分為基質區和優先流區；其中計算單元內基質區所占面積比例為W_m，優先流區所占面積比例為W_f；

根據碎石質量比系數R_v大小，將土壤層分為以下土層：

1)當R_v＝0時，土壤層為均質土壤；

2)當1>R_v>0時，土壤層為土石二元混合介質層；

3)當R_v＝1時，土壤層為基巖層；

此外，假設碎石內和土壤內水勢相等，將基質區土壤水分運動過程進一步修正為：

其中，K_m(h)為基質區內平均導水系數，由基質區內碎石和土壤的導水系數加權平均得到，引入碎石形狀系數ε和質量比系數R_v，得到公式：

K_m(h)＝(1-R_v)K_ss(h)+εR_vK_sr(h) (14)

其中：

式中：w_m為基質區面積比例；h為土壤水吸力；C_ms和C_mr分別為土壤和碎石的容水度；SS為源匯項；Γ為不同區間水量交換量；下標m表示基質區；下標i表示土壤層；K_ss(h)為土壤非飽和導水系數；K_sr(h)為碎石非飽和導水系數；K_ss為土壤飽和導水系數；K_sr為碎石飽和導水系數；α、vn和vm為參數，vm＝1-1/vn；下標1和2分別表示土壤和碎石；z為坐標軸；

土石二元混合介質達到飽和時，h＝0，基于公式(14)，計算單元內飽和導水系數為：

K_m＝(1-R_v)·K_ss+εR_v·K_sr (17)

土壤水分過程計算公式表示如下：

式中：下標f表示優先流區；w_f為優先流區所占面積比例；

上下邊界條件均采用通量邊界條件：

式中：q為各層土壤上下界面處水分通量；

基于Van Genuchten模型分別計算兩者的土壤水分特征曲線：

式中：S_e為飽和度系數；θ為土石二元混合介質土壤含水量；θ_s為土石二元混合介質土壤飽和含水量；θ_r為土石二元混合介質殘余含水量；α、vn和vm為參數，vm＝1-1/vn；h為土壤水吸力；下標1和2分別表示土壤和碎石；

壤中流計算：

R_sub,i＝q_m,i+q_f,i (25)

q_m,i＝w_m(q_m,i-1-q_m,i) (26)

q_f,i＝w_f(q_f,i-1-q_f,i) (27)

式中：Q為壤中過流斷面流量；K_s為土壤飽和導水系數；W為計算單元寬度；Φ為壤中流水位高度，即土水勢；x和z為坐標軸；e為土壤內孔隙度；R_sub為垂向入流量；q為各層土壤上下界面處水分通量；下標m和f分別表示基質區和優先流區；下標i表示土壤層；

地表匯流過程計算，計算公式如下：

連續方法：

運動方程：

S_f＝S₀ (29)

曼寧公式：

式中：Q₀表示地表過流斷面流量；A表示過流斷面面積；R_surf表示地表產流量；S_f表示摩擦坡降；S₀表示計算單元平均地面坡降或河道坡降；R_wr表示過流斷面水力半徑；kn表示曼寧糙率系數；

非降雨期土壤蒸散發過程計算，計算公式如下：

E＝E₁+E₂+E₃ (31)

式中：E表示計算單元總蒸散發量；E₁表示植被截留蒸發量；E₂表示植被蒸騰量；E₃表示裸地土壤蒸發量；

土壤潛在蒸發量由Penman公式計算：

式中：RN為凈輻射量；G為傳入水中的熱通量；Δ為飽和水汽壓對溫度的導數；ρ_a為空氣密度；C_p為空氣的定壓比熱；δ_e為實際水蒸氣壓與飽和水蒸氣壓的差值；r_a為蒸發表面空氣動力學阻抗；λ為水的氣化潛熱；γ為空氣濕度常數；PR為大氣壓；

空氣動力學阻抗計算：其計算公式如下：

式中：r_a為蒸發表面空氣動力學阻抗；hz為氣象站觀測點離地面的高度；hd為置換高度；z_om表為水蒸氣紊流擴散對應的地表粗度；z_ox為地表粗度；κ為von Karman常數；U為風速；h_c為植被高度；

植被截留蒸發量E₁使用Noilhan-Planton模型計算：

E₁＝Veg·δ·E_p (36)

δ＝(Wr/Wr_max)^2/3 (39)

Wr_max＝0.2·Veg·LAI (40)

式中：Veg為裸地-植被域中植被的面積占計算單元的面積比例；δ為濕潤葉面占植被葉面的面積比例；E_p為土壤潛在蒸發量，即最大蒸發量；Wr為植被截留水量；Wr_max為植被最大截留水量；I為時段t_k內最大降水平均雨強；Rr為植被冠層流出水量，即超出最大植被截留水量的部分；LAI表示葉面積指數；

植被蒸騰量E₂采用Penman-Monteith公式計算，土壤各層蒸散量采用雷志棟的根系吸水模型進行計算，公式如下：

E₂＝Veg·(1-δ)·E_PM (41)

式中：E_PM為采用Penman-Monteith公式計算的潛在蒸騰量；G為傳入水中的熱通量；r_c為植被群落阻抗；

根系吸水模型：

T_r＝E₂ (44)

式中：lr表示根系層厚度；T_r為實際植被蒸騰量；SS為源匯項；

植被群落阻抗計算：采用Dickinson公式計算植被群落阻抗：

σ₁^-1＝1-0.0016(25-T_a)² (46)

σ₂^-1＝1-VPD/VPD_c (47)

式中：r_c為植被群落阻抗；r_{s min}為最小氣孔阻抗；LAI為葉面積指數；σ₁為溫度影響函數；σ₂為大氣水蒸氣壓飽和差影響函數；σ₃為光合作用有效放射的影響函數；σ₄為土壤含水量的影響函數；T_a為氣溫；VPD為飽和水蒸氣壓同實測值之間的差；VPD_c為氣孔閉合時的VPD值；r_smax為最大氣孔阻抗；PAR為光和作用有效放射；PAR_c為PAR的臨界值；θ₃為根系層土壤含水量；θ_w為植被凋萎時的土壤含水量；θ_s為土壤飽和含水量；

裸地土壤蒸發量E₃由下式計算：

式中：β為土壤濕潤函數；θ₄為表層土壤含水量；θ_m為土壤分子吸力對應的土壤含水量；θ_h為表層土壤田間持水量；

非降雨期土壤水分運動過程采用修正的Richards公式計算，計算公式同非降雨期；當非降雨期地表有地下水流出時，采用動力波方程進行匯流計算，計算公式同降雨期；非降雨期壤中流采用改進后的動力波方程，計算公式同降雨期。