matlab液体湿润模拟,【水文模型】01 三水源新安江模型
目錄
一、概述
1.1 流域概況
1.2 產(chǎn)流方式論證
1.3 設(shè)計(jì)基本任務(wù)
1.4 設(shè)計(jì)原始資料
二、新安江模型
2.1 流域劃分
2.2 產(chǎn)匯流計(jì)算
2.2.1 蒸散發(fā)計(jì)算
2.2.2 產(chǎn)流計(jì)算
2.2.3 水源劃分
2.2.4 匯流計(jì)算
2.3 模型參數(shù)
2.3.1 參數(shù)物理意義
2.3.2 參數(shù)率定結(jié)果
2.4 模擬結(jié)果
2.4.1 精度統(tǒng)計(jì)表
2.4.2 計(jì)算與實(shí)測流量過程線比較
三、精度統(tǒng)計(jì)與誤差分析
3.1 精度統(tǒng)計(jì)
3.2 誤差分析
3.2.1 模型參數(shù)的影響
3.2.2 人類活動(dòng)的影響
3.2.3 資料代表性的影響
四、參考文獻(xiàn)
五、程序代碼
數(shù)據(jù)文件
一、概述
本博客基于河海大學(xué)水文水資源學(xué)院水文預(yù)報(bào)課程設(shè)計(jì)編寫。采用MATLAB語言編寫三水源新安江模型，并用人機(jī)交互率定模型參數(shù)。模型運(yùn)行所需文件請從新安江模型數(shù)據(jù)下載。
1.1 流域概況
呈村流域控制面積為290 km2，地勢南高北低，相對高差較大，平均海拔高程為583m，流域河道平均坡度為0.95%，最大匯流路徑長度為36 km。流域內(nèi)植被良好，雨量充沛，多年平均降雨量約為2100 mm，流域降水在年內(nèi)年際分配極不均勻，為典型的濕潤流域。該流域植被類型主要包括常綠針葉林、落葉闊葉林、混合林、森林地、林地草原、牧草地與作物地，土壤類型主要為黏壤土。
1.2 產(chǎn)流方式論證
根據(jù)水文氣象條件、下墊面條件及流量過程線分析：
多年平均降雨量約為2100 mm>1000 mm；
流域降水在年內(nèi)年際分配極不均勻，為典型的濕潤流域；
土壤類型主要為黏壤土，土質(zhì)疏松，不易超滲；
綜上，該流域符合蓄滿產(chǎn)流模式。
1.3 設(shè)計(jì)基本任務(wù)
編制預(yù)報(bào)方案，根據(jù)流域概況和原始資料編制預(yù)報(bào)方案；
編寫新安江模型，時(shí)間尺度：日模型(24h)，水源劃分：兩水源(可選)和三水源(可選)；
根據(jù)已給的呈村流域資料，利用編制的新安江模型進(jìn)行日徑流模擬，率定新安江模型參數(shù)；
分析日模型模擬結(jié)果，精度評定時(shí)，日模型采用徑流深相對誤差與確定性系數(shù)。
1.4 設(shè)計(jì)原始資料
逐日資料，包括呈村站流量資料、蒸發(fā)資料，呈村、汪村、樟源口、棣甸、董坑塢、用功城、左龍、馮村、田里與大連站雨量資料。資料起止時(shí)間見表1-1。
參數(shù)建議取值見表1-2。
各雨量站控制面積及比例，見表1-3。
二、新安江模型
2.1 流域劃分
根據(jù)流域地形、地貌條件及布設(shè)的雨量站網(wǎng)，用泰森多邊形法將呈村流域劃分為10塊單元面積，各單元面積的權(quán)重見表1-3。
2.2 產(chǎn)匯流計(jì)算
對每塊單元面積采用三水源新安江模型分別進(jìn)行蒸散發(fā)計(jì)算、產(chǎn)流計(jì)算、水源劃分和匯流計(jì)算，得到單元面積出流過程；由于本流域面積較小，單元面積的出流過程不采用馬斯京根分段連續(xù)演算法進(jìn)行出口以下的河道洪水演算，而是通過調(diào)整滯時(shí) L 求得單元面積在流域出口的流量過程線；將每個(gè)單元面積在流域出口的流量過程線線性疊加，即為呈村流域的出口流量過程。
2.2.1 蒸散發(fā)計(jì)算
(1)當(dāng)W U + P ≥ E P WU+P≥E_PWU+P≥EP?時(shí)
E U = E P ， E L = 0 ， E D = 0 E_U=E_P，E_L=0，E_D=0EU?=EP?，EL?=0，ED?=0
(2)當(dāng)W U + P < E P ， W L ≥ C ? W L M WU+PWU+P
E U = W U + P ， E L = ( E P ? E U ) W L / W L M ， E D = 0 E_U=WU+P，E_L=(E_P-E_U)WL/WLM，E_D=0EU?=WU+P，EL?=(EP??EU?)WL/WLM，ED?=0
(3)當(dāng)W U + P < E P ， C ( E P ? E U ) ≤ W L < C ? W L M WU+PWU+P
E U = W U + P ， E L = C ( E P ? E U ) ， E D = 0 E_U=WU+P，E_L=C(E_P-E_U)，E_D=0EU?=WU+P，EL?=C(EP??EU?)，ED?=0
(4)當(dāng)W U + P < E P ， W L < C ( E P ? E U ) WU+PWU+P
E U = W U + P ， E L = W L ， E D = C ( E P ? E U ) ? E L E_U=WU+P，E_L=WL，E_D=C(E_P-E_U )-E_LEU?=WU+P，EL?=WL，ED?=C(EP??EU?)?EL?
式中：
W U WUWU——上層土壤含水量，mm；
W L WLWL——下層土壤含水量，mm；
P PP——降雨量，mm；
E P E_PEP?——流域蒸發(fā)能力，mm；
E U E_UEU?——上層蒸發(fā)量，mm；
E L E_LEL?——下層蒸發(fā)量，mm；
E D E_DED?——深層蒸發(fā)量，mm；
C CC——蒸發(fā)擴(kuò)散系數(shù)；
W L M WLMWLM——下層土壤含水容量，mm。
2.2.2 產(chǎn)流計(jì)算
流域平均土壤水蓄積容量關(guān)系
W M = W M M 1 + b WM=\frac{WMM}{1+b}WM=1+bWMM?
a = W M M [ 1 ? ( 1 ? W W M ) 1 1 + b ] a=WMM[1-(1-\frac{W}{WM})^{\frac{1}{1+b}} ]a=WMM[1?(1?WMW?)1+b1?]
(1)當(dāng)P E ≤ 0 PE≤0PE≤0時(shí)
R = 0 R=0R=0
(2)當(dāng)P E > 0 ， a + P E ≤ W M M PE>0，a+PE≤WMMPE>0，a+PE≤WMM時(shí)
R = P E + W ? W M + W M ( 1 ? P E + a W M M ) b + 1 R=PE+W-WM+WM(1-\frac{PE+a}{WMM})^{b+1}R=PE+W?WM+WM(1?WMMPE+a?)b+1
(3)當(dāng)P E > 0 ， a + P E > W M M PE>0，a+PE>WMMPE>0，a+PE>WMM時(shí)
R = P E + W ? W M R=PE+W-WMR=PE+W?WM
式中
W M WMWM——流域平均蓄水容量，mm；
W M M WMMWMM——流域最大蓄水容量，mm；
a aa——流域最大初始土壤含水量，mm；
b bb——常數(shù)，反映流域包氣帶蓄水容量分布的不均勻性
W WW——流域平均初始土壤含水量，mm；
P E PEPE——扣除雨期蒸發(fā)后的降雨量，mm；
R RR——產(chǎn)流量，mm。
2.2.3 水源劃分
(1)二水源劃分
α = R P E α=\frac{R}{PE}α=PER?
當(dāng)PE>FC時(shí)
R G = α ? F C RG=α?FCRG=α?FC
R S = α ? ( P E ? F C ) RS=α?(PE-FC)RS=α?(PE?FC)
當(dāng)P E ≤ F C PE≤FCPE≤FC時(shí)
R G = R RG=RRG=R
R S = 0 RS=0RS=0
式中
α \alphaα——產(chǎn)流面積比例；
P E PEPE——扣除雨期蒸發(fā)后的降雨量，mm；
R RR——產(chǎn)流量，mm；
F C FCFC——流域的穩(wěn)定下滲量，mm；
R G RGRG——地下徑流，mm；
R S RSRS——地面徑流，mm。
(2)三水源劃分
流域平均自由水蓄積容量關(guān)系
S m = S m m 1 + E X S_m=\frac{S_{mm}}{1+EX}Sm?=1+EXSmm??
A U = S m m [ 1 ? ( 1 ? S 1 F R 1 F R S m ) 1 1 + E X ] AU=S_{mm} [1-(1-\frac{S_1\frac{FR_1}{FR}}{S_m})^{\frac{1}{1+EX}}]AU=Smm?[1?(1?Sm?S1?FRFR1???)1+EX1?]
當(dāng)P E + A U < S m m PE+AUPE+AU
R S = F R [ P E + S 1 ? F R 1 F R ? S m + S m ( 1 ? P E + A U S m m ) E X + 1 ] RS=FR[PE+\frac{S_1?FR_1}{FR}-S_m+S_m (1-\frac{PE+AU}{S_{mm}} )^{EX+1} ]RS=FR[PE+FRS1??FR1???Sm?+Sm?(1?Smm?PE+AU?)EX+1]
當(dāng)P E + A U ≥ S m m PE+AU≥S_{mm}PE+AU≥Smm?時(shí)
R S = F R ( P E + S 1 ? F R 1 F R ? S m ) RS=FR(PE+\frac{S_1?FR_1}{FR}-S_m )RS=FR(PE+FRS1??FR1???Sm?)
本時(shí)段的自由水蓄量為
S = S 1 ? F R 1 F R + R ? R S F R S=\frac{S_1?FR_1}{FR}+\frac{R-RS}{FR}S=FRS1??FR1??+FRR?RS?
相應(yīng)的壤中流和地下徑流為
R I = K I ? S ? F R RI=KI?S?FRRI=KI?S?FR
R G = K G ? S ? F R RG=KG?S?FRRG=KG?S?FR
本時(shí)段末的自由水蓄量變?yōu)?/p>
S 1 = S ( 1 ? K I ? K G ) S_1=S(1-KI-KG)S1?=S(1?KI?KG)
式中
F R 1 、 F R FR_1 、FRFR1?、FR——上一時(shí)段和本時(shí)段的產(chǎn)流面積比例；
S m S_mSm?——流域平均自由水蓄量，mm；
S m m S_{mm}Smm?——流域最大自由水蓄量，mm；
A U AUAU——流域最大初始自由水蓄量，mm；
E X EXEX——自由水蓄量分布曲線指數(shù)
P E PEPE——扣除雨期蒸發(fā)后的降雨量，mm；
R G RGRG——地下徑流，mm；
R S RSRS——地面徑流，mm；
R I RIRI——壤中流，mm；
R RR——產(chǎn)流量，mm；
S 1 S_1S1?——時(shí)段初始自由水蓄量，mm；
S SS——本時(shí)段的自由水蓄量，mm。
2.2.4 匯流計(jì)算
(1)地面徑流的坡地匯流
地面徑流的坡地匯流時(shí)間不計(jì)，直接進(jìn)入河網(wǎng)，計(jì)算公式為：
Q S ( I ) = R S ( I ) × U QS(I)=RS(I)×UQS(I)=RS(I)×U
式中：
Q S QSQS——地面徑流，m3/s；
U UU——單位轉(zhuǎn)換系數(shù)；
R S RSRS——地面徑流量，mm。
(2)壤中流匯流
表層自由水以K I KIKI側(cè)向出流后成為表層壤中流，進(jìn)入河網(wǎng)。但如土層較厚，表層自由水尚可滲入深層土，經(jīng)過深層土的調(diào)蓄作用，才進(jìn)入河網(wǎng)。深層自由水用線性水庫模擬，其消退系數(shù)為C I CICI，計(jì)算公式為：
Q I ( I ) = C I × Q I ( I ? 1 ) + ( 1 ? C I ) × R I ( I ) × U QI(I)=CI×QI(I-1)+(1-CI)×RI(I)×UQI(I)=CI×QI(I?1)+(1?CI)×RI(I)×U
式中：
Q I QIQI——壤中流，m3/s；
C I CICI——壤中流消退系數(shù)；
R I RIRI——壤中流徑流量，mm。
(3)地下徑流匯流
地下徑流匯流用線性水庫模擬,其消退系數(shù)為C G CGCG，出流進(jìn)入河網(wǎng)。表層自由水以K G KGKG向下出流后，再向地下水庫匯流的時(shí)間不另計(jì)，包括在C G CGCG之內(nèi)，計(jì)算公式為：
Q G ( I ) = C G × Q G ( I ? 1 ) + ( 1 ? C G ) × R G ( I ) × U QG(I)=CG×QG(I-1)+(1-CG)×RG(I)×UQG(I)=CG×QG(I?1)+(1?CG)×RG(I)×U
式中：
Q G QGQG——地下徑流，m3/s；
C G CGCG——地下水消退系數(shù)；
R G RGRG——地下徑流徑流量，mm。
(4)單元面積河網(wǎng)匯流
單元面積的河網(wǎng)匯流用滯后演算法，參數(shù)有滯后量L LL與消退系數(shù)C S CSCS，計(jì)算公式為：
Q T ( I ) = Q S ( I ) + Q I ( I ) + Q G ( I ) QT(I)=QS(I)+QI(I)+QG(I)QT(I)=QS(I)+QI(I)+QG(I)
Q ( I ) = C S × Q ( I ? 1 ) + ( 1 ? C S ) × Q T ( I ? L ) Q(I)=CS×Q(I-1)+(1-CS)×QT(I-L)Q(I)=CS×Q(I?1)+(1?CS)×QT(I?L)
式中：
Q T QTQT——總徑流，m3/s；
C S CSCS——河網(wǎng)水流消退系數(shù)；
Q QQ——單元面積河網(wǎng)總?cè)肓?#xff0c;m3/s；
L LL——河網(wǎng)匯流滯時(shí)，d。
2.3 模型參數(shù)
新安江模型是一個(gè)通過長期實(shí)踐和對水文規(guī)律認(rèn)識基礎(chǔ)上建立起來的一個(gè)概念性水文模型。模型大多數(shù)參數(shù)都具有明確的物理意義，它們在一定程度上反映了流域的基本水文特征和降雨徑流形成的物理過程。因此，原則上可以按其物理意義通過實(shí)測、實(shí)驗(yàn)、比擬等方法確定。本課程設(shè)計(jì)采用參數(shù)的概念分析方法，即先按實(shí)測值或參數(shù)的物理意義初定參數(shù)初值范圍；然后根據(jù)輸入，通過模型計(jì)算輸出；再將輸出過程與實(shí)測過程進(jìn)行比較，作優(yōu)化調(diào)試；根據(jù)特定的目標(biāo)函數(shù)——年徑流相對誤差和確定性系數(shù)確定參數(shù)的最優(yōu)值。
2.3.1 參數(shù)物理意義
(1)蒸散發(fā)能力折算系數(shù)K KK
K是影響產(chǎn)流量計(jì)算最為重要和敏感的參數(shù)，產(chǎn)流計(jì)算中K控制著水量平衡。K主要反映流域平均高程和蒸發(fā)站高程之間差別的影響和蒸發(fā)皿蒸散發(fā)與陸面蒸散發(fā)間差別的影響。蒸散發(fā)能力的地區(qū)分布大體上反映了氣候和自然地理?xiàng)l件的影響，具有較為明顯的區(qū)域性規(guī)律。
(2)自由水蓄水容量S M SMSM
SM反映表土蓄水能力，其值受降雨資料時(shí)段均化的影響明顯。當(dāng)以日作為時(shí)段長時(shí)，在土層很薄的山區(qū)，其值為10mm或更小一些；而在土深林茂透水性很強(qiáng)的流域，其值可取50mm或更大一些；一般流域在10~20mm之間。SM對地面徑流和地下徑流的比重起著決定性作用；SM大，則地下徑流所占比重相對大，地面徑流所占比重相對小，洪峰流量相對小；反之，SM小，則地下徑流所占比重相對小，地面徑流所占比重相對大，洪峰流量相對大。
(3)自由水蓄水庫對地下水和壤中流的日出流系數(shù)K G 、 K I KG、KIKG、KI
KG的大小反映基巖和深層土壤的滲透性，KI的大小反映表層土的滲透性。KG+KI代表出流的快慢，KG/KI代表地下徑流與壤中流的比，對于一個(gè)特定流域它們都是常數(shù)，1-(KG+KI)為消退系數(shù)，它決定了直接徑流的退水快慢。
(4)地下水消退系數(shù)C G CGCG
CG可根據(jù)枯季地下徑流的退水規(guī)律來推求，C G = Q t + ? t Q t CG=\frac{Q_{t+?t}}{Q_t}CG=Qt?Qt+?t??。若以日作為計(jì)算時(shí)段長，則CG=0.950-0.998，大致相當(dāng)于消退歷時(shí)為20~500d。
(5)壤中流消退系數(shù)C I CICI
若無壤中流則CI→0，若壤中流豐富則CI→0.9，相當(dāng)于匯流時(shí)間為10d。
(6)河網(wǎng)水流消退系數(shù)C S CSCS
CS代表坦化作用，其值取決于河網(wǎng)的地貌條件，可通過河網(wǎng)地貌推求。因與時(shí)段長短有關(guān)，其值應(yīng)視河道特征和洪水特性而定。
(7)河網(wǎng)匯流滯時(shí)L LL
L代表平移作用，其值取決于河網(wǎng)的地貌條件，可通過河網(wǎng)地貌推求。
2.3.2 參數(shù)率定結(jié)果
原始資料中給定的參數(shù)值如表2-1所示。
利用人機(jī)交互率定方法得到的參數(shù)率定結(jié)果如下表2-2所示。
2.4 模擬結(jié)果
2.4.1 精度統(tǒng)計(jì)表
呈村流域日模型模擬結(jié)果及精度統(tǒng)計(jì)見表2-3。
2.4.2 計(jì)算與實(shí)測流量過程線比較
其余年份比較圖見附圖。另因計(jì)算時(shí)間尺度長，僅附1989年逐日計(jì)算數(shù)據(jù)，見附表。
三、精度統(tǒng)計(jì)與誤差分析
3.1 精度統(tǒng)計(jì)
從表2-3中可見，年產(chǎn)流量絕對誤差小于100mm的有7次，占總數(shù)的87.5%，絕對誤差小于45mm的有4次，占總數(shù)的50%；所有年份產(chǎn)流量的相對誤差均小于10%；確定性系數(shù)最大為0.9006，最小為0.6872，平均為0.8317。
精度統(tǒng)計(jì)表明，率定的模型參數(shù)基本上是合理的。
3.2 誤差分析
影響流域降雨徑流過程的因素很多，三水源新安江模型的結(jié)構(gòu)和參數(shù)能夠反映濕潤地區(qū)降雨徑流過程的主要規(guī)律和特點(diǎn)，因而能獲得較好的精度。但是，模型本身以及模型計(jì)算中有許多的概化，會(huì)造成誤差。造成呈村流域日模型預(yù)報(bào)方案誤差來源主要有以下三個(gè)方面：
3.2.1 模型參數(shù)的影響
模型參數(shù)是根據(jù)輸入，通過模型計(jì)算輸出，再將輸出過程與實(shí)測過程進(jìn)行比較，用人工調(diào)試的方法進(jìn)行優(yōu)化的，所率定出的參數(shù)可能不是最優(yōu)的。
3.2.2 人類活動(dòng)的影響
隨著經(jīng)濟(jì)建設(shè)的發(fā)展，人類活動(dòng)的影響加劇，流域內(nèi)可能建有塘壩、橡膠壩等蓄水水利工程。這些水利工程沒有固定的調(diào)度規(guī)則和調(diào)洪方式，干旱季節(jié)或汛初時(shí)，流域內(nèi)的水利工程均會(huì)蓄水，會(huì)使實(shí)測徑流量偏小。
3.2.3 資料代表性的影響
本次課程設(shè)計(jì)采用的樣本僅有8年實(shí)測水文系列資料，由于水文自然規(guī)律的隨時(shí)間變化，使觀測到的水文氣象資料代表性不夠，難免造成預(yù)報(bào)方案的誤差。
四、參考文獻(xiàn)
[1] 包為民.水文預(yù)報(bào).北京：中國水利水電出版社，2017.
[2] 李巧玲.水文預(yù)報(bào)課程設(shè)計(jì)任務(wù)書.南京：河海大學(xué)水文水資源學(xué)院，2019．
[3] 翟家瑞.常用水文預(yù)報(bào)算法和計(jì)算程序. 河南：黃河水利出版社，1995
五、程序代碼
clear;
% data=xlsread('呈村流域資料.xls');
data=xlsread('1995參數(shù)檢驗(yàn).xls');
data1=xlsread('呈村流域資料.xls',2); %載入原始資料數(shù)據(jù)
QZ=data(:,2); %QZ表示實(shí)測流量
%%%%%%%%%%%%%需率定的參數(shù)%%%%%%%%%%%%%
par=load('input.txt');
K=par(1); %蒸散發(fā)折算系數(shù) 范圍：0.8～1.2
SM=par(2); %自由水容量 范圍：10～30
KG=par(3); %地下水出流系數(shù) 范圍: 0.01～0.69
CG=par(4); %地下水消退系數(shù) 范圍：0.98～0.998
CI=par(5); %壤中流消退系數(shù) 范圍：0.01～0.9
CS=par(6); %河網(wǎng)水流消退系數(shù) 范圍：0.1～0.9
L=1; %河網(wǎng)匯流滯時(shí) 范圍：0～2
%%%%%%%%%%%%%需率定的參數(shù)%%%%%%%%%%%%%
E0=data(:,3); %E0表示呈村流域器皿蒸發(fā)能力過程,位于原始數(shù)據(jù)表第三列
EP=K*E0; %計(jì)算流域蒸發(fā)能力
F=data1(:,3); %讀取各子流域面積
Fr=data1(:,4); %讀取各子流域面積比例
U=F/(3.6*24); %U為單位轉(zhuǎn)換系數(shù)
WM=120;WUM=20;WLM=60;WDM=WM-WUM-WLM;C=0.18;B=0.3;FC=2.5;FR0=0;PTT=0;%PTT表示流域年降雨量 %常數(shù)
EX=1.5;IM=0.01;KI=0.7-KG;WMM=(1+B)*WM;SMM=(1+EX)*SM;SZ=length(E0); %常數(shù)
S1=zeros(SZ,1);S1(1)=SM; %為S1預(yù)分配內(nèi)存；S1表示本時(shí)段初產(chǎn)流面積上的平均自由水深
FR=zeros(SZ,1); %為FR預(yù)分配內(nèi)存；FR表示本時(shí)段產(chǎn)流面積比例
EU=zeros(SZ,1); %為EU預(yù)分配內(nèi)存
EL=zeros(SZ,1); %為EL預(yù)分配內(nèi)存
ED=zeros(SZ,1); %為ED預(yù)分配內(nèi)存
E=zeros(SZ,1); %為E預(yù)分配內(nèi)存
WU=zeros(SZ,1); %為WU預(yù)分配內(nèi)存
WL=zeros(SZ,1); %為WL預(yù)分配內(nèi)存
WD=zeros(SZ,1); %為WD預(yù)分配內(nèi)存
W=zeros(SZ,1); %為W預(yù)分配內(nèi)存
R=zeros(SZ,1); %為R預(yù)分配內(nèi)存
RG2=zeros(SZ,1); %為RG2預(yù)分配內(nèi)存;RG2表示二水源地下徑流
RD2=zeros(SZ,1); %為RD2預(yù)分配內(nèi)存；RD2表示二水源地面徑流
RS3=zeros(SZ,1); %為RS3預(yù)分配內(nèi)存；RS3表示三水源地面徑流
RI3=zeros(SZ,1); %為RI3預(yù)分配內(nèi)存；RI3表示三水源壤中水徑流
RG3=zeros(SZ,1); %為RG3預(yù)分配內(nèi)存；RG3表示三水源地下徑流
PE=zeros(SZ,1); %為PE預(yù)分配內(nèi)存；PE表示凈雨量
QI3=zeros(SZ,1); %為QI3預(yù)分配內(nèi)存；QI3表示三分水源壤中流匯流
QG3=zeros(SZ,1); %為QG3預(yù)分配內(nèi)存；QG3表示三分水源地下徑流匯流
Q3=zeros(SZ,1); %為Q3預(yù)分配內(nèi)存；Q3表示單元面積河網(wǎng)匯流
Q=zeros(SZ,1); %為Q預(yù)分配內(nèi)存；Q表示呈村流域出口斷面流量過程
TTT=zeros(SZ,18); %為TTT預(yù)分配內(nèi)存；
WU(1)=0; %流域上層初始土壤含水量 范圍：0～20
WL(1)=50;WD(1)=WDM; %WL、WD都表示本時(shí)段初相應(yīng)土層土壤含水量
W(1)=WU(1)+WL(1)+WD(1); %W表示本時(shí)段初土壤含水量
%%%%%%%%%%%%%計(jì)算10個(gè)子流域%%%%%%%%%%%%%
for J=1:10
P=(1-IM)*data(:,J+3); %P表示某子流域降雨過程，data數(shù)組中第4列到第13列分別為
% 呈村、汪村、樟源口、棣甸、董坑塢、用功城、左龍、馮村、田里、大連的降雨資料
RB=IM*data(:,J+3); %RB表示降落在不透水面上的降雨量
PTT=PTT+sum(data(:,J+3))*Fr(J);
for I=1:SZ
%%%%%%%%%%%%%三層蒸散發(fā)模型%%%%%%%%%%%%%
if WU(I)+P(I)>=EP(I)
EU(I)=EP(I);EL(I)=0;ED(I)=0;
else if WL(I)>=C*WLM
EU(I)=WU(I)+P(I);EL(I)=(EP(I)-EU(I))*WL(I)/WLM;ED(I)=0;
else if WL(I)>=C*(EP(I)-(WU(I)+P(I)))
EU(I)=WU(I)+P(I);EL(I)=C*(EP(I)-EU(I));ED(I)=0;
else
EU(I)=WU(I)+P(I);EL(I)=WL(I);ED(I)=C*(EP(I)-EU(I))-EL(I);
end
end
end
E(I)=EU(I)+EL(I)+ED(I);PE(I)=P(I)-E(I);
%%%%%%%%%%%%%三層蒸散發(fā)模型%%%%%%%%%%%%%
%%%%%%%%%%%%%計(jì)算總徑流深%%%%%%%%%%%%%
a=WMM*(1-(1-W(I)/WM)^(1/(B+1))); %a表示流域初始土壤含水量最大值
if PE(I)<=0
R(I)=0;
else if a+PE(I)<=WMM
R(I)=PE(I)+W(I)-WM+WM*(1-(PE(I)+a)/WMM)^(B+1);
else
R(I)=PE(I)+W(I)-WM;
end
end
%%%%%%%%%%%%%計(jì)算總徑流深%%%%%%%%%%%%%
%%%%%%%%%%%%%計(jì)算下一時(shí)段土壤含水量%%%%%%%%%%%%%
WU(I+1)=WU(I)+P(I)-EU(I)-R(I);
WL(I+1)=WL(I)-EL(I);
WD(I+1)=WD(I)-ED(I);
if WU(I+1)>WUM
WL(I+1)=WL(I)-EL(I)+WU(I+1)-WUM;
WU(I+1)=WUM;
end
if WL(I+1)>WLM
WD(I+1)=WD(I)-ED(I)+WL(I+1)-WLM;
WL(I+1)=WLM;
end
if WD(I+1)>WDM
WD(I+1)=WDM;
end
W(I+1)=WU(I+1)+WL(I+1)+WD(I+1);
%%%%%%%%%%%%%計(jì)算下一時(shí)段土壤含水量%%%%%%%%%%%%%
%%%%%%%%%%%%%計(jì)算產(chǎn)流面積FR=R/PE%%%%%%%%%%%%%
if PE(I)==0
FR(I)=0;
else if R(I)/PE(I)>1
FR(I)=1;
else
FR(I)=R(I)/PE(I);
end
end
%%%%%%%%%%%%%計(jì)算產(chǎn)流面積FR=R/PE%%%%%%%%%%%%%
%%%%%%%%%%%%%二分水源計(jì)算%%%%%%%%%%%%%
if PE(I)<=FC
RG2(I)=R(I);RD2(I)=0;
else
RG2(I)=FC*FR(I);RD2(I)=(PE(I)-FC)*FR(I);
end
%%%%%%%%%%%%%二分水源計(jì)算%%%%%%%%%%%%%
%%%%%%%%%%%%%三分水源計(jì)算%%%%%%%%%%%%%
%%%%%%%%%%%%%第一時(shí)段數(shù)據(jù)計(jì)算%%%%%%%%%%%%%
if I==1
AU=SMM*(1-(1-(S1(I)*FR0/FR(I))/SM)^(1/(1+EX)));
if PE(I)<=0
RS3(I)=0;RI3(I)=0;RG3(I)=0;S1(I+1)=S1(I)*(1-KI-KG);
else if PE(I)+AU
數(shù)據(jù)文件
模型運(yùn)行所需文件請從新安江模型數(shù)據(jù)下載。
總結(jié)
以上是生活随笔為你收集整理的matlab液体湿润模拟,【水文模型】01 三水源新安江模型的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： VS 2010 OpenGL 配置与实例
- 下一篇： 5.1 入门整合案例(SpringBoo