水文預(yù)報課程設(shè)計報告-江西滁州流域不同參數(shù)組合下產(chǎn)流量模  擬效果對比分析

水文預(yù)報課程設(shè)計報告心得體會1緒論1.1設(shè)計任務(wù)根據(jù)已給的降水、蒸發(fā)資料，選擇不同的產(chǎn)流參數(shù)范圍，編寫程序，將滁州流城作為整體進行產(chǎn)流量計算；分析不同參數(shù)組合對產(chǎn)流量計算結(jié)果的影響程度。1.2設(shè)計目的1、通過課程設(shè)計熟悉降水徑流預(yù)報的內(nèi)容，熟悉新安江模型的原理和結(jié)構(gòu)，掌握降水徑流預(yù)報方案編制的基本步驟；2、掌握參數(shù)的物理意義，了解產(chǎn)流參數(shù)大小設(shè)置對結(jié)果的影響。3、加深對水文預(yù)報方法的掌握，學(xué)會各種方法的綜合運用；4、鍛煉從事專業(yè)工作的基本能力，培養(yǎng)分析問題和解決問題的能力。5、掌握計算機編程的基本方法和基本操作；熟練OFFICE文檔操作和研究報告的規(guī)范編制。1.3技術(shù)要求1、采用泰森單元法進行面平均降水量計算；2、產(chǎn)流采用三江源新安江模型；3、每人計算3種組合情景下的產(chǎn)流量，并進行對比。分析參數(shù)改變對計算徑流量的影響。1.4工作步驟1、產(chǎn)流量計算；2、劃分水源，得出三種水源徑流量(徑流深)；3、根據(jù)參數(shù)設(shè)置，分析參數(shù)改變對徑流深的影響。1.5設(shè)計要求根據(jù)已給資料、參數(shù)，自己獨立編寫程序，將流域作為整體進行產(chǎn)流量計算、水源劃分;繪制徑流過程；給出參數(shù)組合情景下的產(chǎn)流量對比結(jié)果，分析參數(shù)設(shè)置對產(chǎn)流量的影響。設(shè)計時間為2周，要求獨立完成，遇到問題可以相互討論，人工計算與計算機計算相結(jié)合。所有數(shù)據(jù)資料共享，成果表達方式要求多樣化(包含文字、圖、表等)。報告撰寫要認真，不得抄襲?？梢韵热斯び嬎阋欢〞r間長度的降水徑流過程，以便于計算機編程并作程序調(diào)試。1.6有關(guān)說明1、程序代碼編制中建議變量名：PJ表平均雨量，Ep表示蒸散發(fā)能力；其他參數(shù)名稱設(shè)置參考教材。2材料與方法2.1流域基本情況2.1.1自然地理概況遂川江流域主要位于江西省遂川縣境內(nèi)，是贛江一級支流，是吉安市五大河流之一，發(fā)源于湖南省桂東縣北部的龍?zhí)赌X，河源位于東經(jīng)113o56′，北緯26o11′。遂川江流域以山地和丘陵為主，流域內(nèi)植被較好，森林覆蓋率達70%。流域上游為山區(qū)，河床多為礫石、巖石，河面較窄，比降較大，水淺流速大，中下游為低丘區(qū)，河床多為礫沙、卵石、細沙。遂川流域?qū)僦衼啛釒駶櫦撅L氣候區(qū)。遂川江由左、右溪組成，從西南向東北縱貫全線，總長249.5km，流域面積2824km2，流域平均高程511m，流域平均坡度0.31m/km2，流域長度100km，流域形狀系數(shù)0.29，主河道長度176.0km，河道縱比降2.36‰。2.1.2水文特征流域多年平均降水量為1637mm。降水的年際變化較大，降水主要集中在汛期4～9月，特別其中在4、5、6月，即主汛期，約占全年降水的50%。年降水天數(shù)為138～183天，降雨量總的分布是東部和西部大于中部，北部大于南部。多年平均蒸發(fā)量為700～1100mm，蒸發(fā)量的年變化與氣溫的年變化趨勢大體一致，一般是12～1月份的蒸發(fā)量為最小，7～9月份的蒸發(fā)量為最大。多年平均徑流量1416×108m3，多年平均徑流深847mm。徑流年內(nèi)分配是：從1月的2～3%開始逐月上升，5月或6月達到18～25%，從6月或7月開始逐月下降，12月約占2～3%，徑流年際變化較大。2.1.3暴雨洪水特征暴雨是江西省遂川縣主要災(zāi)害之一，以4～8月出現(xiàn)較多。暴雨過程一般持續(xù)1～2d，個別年份暴雨帶在流域內(nèi)南北擺動，持續(xù)10d以上，例如1962年6月，1998年7月。暴雨強度一般是日降水量為50～100mm。遂川縣各河流主汛期大致在4～7月，各站最大洪峰流量和最大峰現(xiàn)時間多發(fā)生在5月和6月，干流各站一次洪水過程多表現(xiàn)為峰高量大，歷時一般為15～20d，單峰上游出現(xiàn)較多，歷時一般為1～3d，由于暴雨帶的南北擺動，暴雨范圍大、強度大、歷時長，往往遂川江的洪水為幾條河流的洪水同時組成。2.1.4水文資料遂川江流域面積288.24km2，包括小夏、桐古、營盤墟、沙湖里、阡陌、淋洋、七嶺、滁州8個水文站的集水面積。本次課程設(shè)計選取滁州水文站（東經(jīng)114o08′，北緯26o21′）以上集水區(qū)間。站點分布情況參見圖1和圖2。站點權(quán)重見表1。模型參數(shù)取值范圍見表2。圖1遂川江流域圖圖2遂川江流域滁州區(qū)間泰森多邊形表1遂川江流域滁州區(qū)間泰森多邊形面積權(quán)重站號站名面積（km2）面積權(quán)重1小夏44.9115.58%2桐古14.715.10%3營盤墟42.2814.67%4沙湖里37.7213.09%5阡陌43.2114.99%6淋洋53.7418.65%7七嶺27.519.54%8滁州24.168.38%合計288.24100%表2模型參數(shù)取值范圍KCWMWUMWLMWDMBCIMPSMEX0.901302540650.20.120.0013510.951402550650.250.140.002401.21.001202560350.30.160.003451.3初始張力水含量初始SWWUWLWDS110104060352.2新安江流域水文模型2.2.1模型介紹1973年，河海大學(xué)趙人俊教授領(lǐng)導(dǎo)的研究組在編制新安江洪水預(yù)報方案時，匯集了當時在產(chǎn)匯流理論方面的研究成果，并結(jié)合大流域洪水預(yù)報的特點，設(shè)計了國內(nèi)第一個完整的流域水文模型——新安江流域水文模型，以下簡稱新安江模型[]。最初研制的是二水源新安江模型，20世紀80年代中期，借鑒山坡水文學(xué)的概念和國內(nèi)外產(chǎn)匯流理論的研究成果，提出了三水源新安江模型。三水源新安江模型蒸散發(fā)計算采用三層模型；產(chǎn)流計算采用蓄滿產(chǎn)流模型；用自由水蓄水庫結(jié)構(gòu)將總徑流劃分為地表徑流、壤中流和地下徑流3種；流域匯流計算采用線性水庫；河道匯流采用馬斯京根分段連續(xù)演算或滯后演算法[]。2.2.2模型結(jié)構(gòu)為了考慮降水和流域下墊面分布不均勻的影響，新安江模型的結(jié)構(gòu)設(shè)計為分散性的，分為蒸散發(fā)計算、產(chǎn)流計算、分水源計算和匯流計算四個層次結(jié)構(gòu)。每塊單元流域的計算流程如圖3所示。圖3三水源新安江模型流程圖圖3中方框外為參數(shù)，方框內(nèi)為狀態(tài)變量。模型各層次結(jié)構(gòu)的功能、計算采用的方法和相應(yīng)參數(shù)見表3。表3新安江模型各層次結(jié)構(gòu)功能、計算采用的方法和相應(yīng)參數(shù)表層次（第一層次）（第二層次）（第三層次）（第四層次）功能蒸散發(fā)計算產(chǎn)流計算水源劃分匯流計算二水源三水源坡面匯流河道匯流方法三層模型蓄滿產(chǎn)流穩(wěn)定下滲率自由水蓄水庫單位線或線性水庫或滯后演算法馬斯京根或滯后演算法參數(shù)KC、UM、LM、CWM、B、IMFCSM、EX、KG、KIUH或CS、CI、CGKE、XE或L2.2.3計算方法2.2.3.1降雨量計算為了考慮降雨分布不均勻和下墊面分布的不均勻性，采用泰森多邊形法將遂川江流域滁州區(qū)間劃分為八個單元流域，在每塊單元流域內(nèi)各有一個雨量站，以每個多邊形內(nèi)雨量站的雨量代表該多邊形面積上的降雨量，最后按面積加權(quán)推求流域平均降雨量。計算公式如下：P=式中：P為流域平均降雨量；Pi為流域內(nèi)第i個雨量站同時段降雨量；f為第i個雨量站所在多邊形的面積；F2.2.3.2蒸發(fā)量計算在新安江模型中，流域蒸散發(fā)計算沒有考慮流域內(nèi)土壤含水量在面上分布的不均勻性，而是按土壤垂向分布的不均勻性將土層分為三層，用三層蒸散發(fā)模型計算蒸散發(fā)量。參數(shù)有流域平均張力水容量WM（mm），上層張力水容量UM（mm），下土層張力水容量LM（mm），深層張力水容量DM（mm），蒸散發(fā)折算系數(shù)KC和深層蒸散發(fā)擴散系數(shù)C，計算公式如下：WM=W=WUE=EP式中：W為總的張力水蓄量，mm；WU為上層張力水蓄量，mm；WL為下層張力水蓄量，mm；WD為深層張力水蓄量，mm；E為總的蒸散發(fā)量，mm；EU為上層蒸散發(fā)量，mm；EL為下層蒸散發(fā)量，mm；ED為深層蒸散發(fā)量，mm；EP為蒸散發(fā)能力，mm；蒸散發(fā)折算系數(shù)三層蒸發(fā)模式按照先上層后下層的次序，具有分以下四種情況計算：（1）當WU+EU（2）當WU+EU（3）當WU+EU（4）當WU+EU2.2.3.3產(chǎn)流計算產(chǎn)流計算中采用蓄滿產(chǎn)流機制。蓄滿是指包氣帶的土壤含水量達到田間持水量。蓄滿產(chǎn)流是指：降水在滿足田間持水量以前不產(chǎn)流，所有的降水都被土壤所吸收；降水在滿足田間持水量以后，所有的降水（扣除同期蒸發(fā)量）都產(chǎn)流。其概念就是設(shè)想流域具有一定的蓄水能力，當這種蓄水能力滿足以后，全部降水變?yōu)閺搅?，產(chǎn)流表現(xiàn)為蓄量控制的特點。濕潤地區(qū)產(chǎn)流的蓄量控制特點，解決了產(chǎn)流計算在這些地區(qū)處理雨強和入滲動態(tài)過程的問題;而降雨徑流理論關(guān)系的建立，解決了考慮流域降雨不均勻的分布式產(chǎn)流計算問題。按照蓄滿產(chǎn)流的概念，采用蓄水容量-面積分配曲線來考慮土壤缺水量分布不均勻的問題。所謂蓄水容量—面積分配曲線是指：部分產(chǎn)流面積隨蓄水容量而變化的累計頻率曲線。應(yīng)用蓄容量-面積分配曲線可以確定降雨空間分布不均勻情況下蓄滿產(chǎn)流的總徑流量。對于閉合流域，流域蓄水容量-面積分配曲線采用拋物線形為宜。流域蓄水容量-面積分配曲線與降雨徑流相互轉(zhuǎn)換關(guān)系如圖4所示。(a)(b)圖4流域蓄水容量-面積分配曲線與降雨徑流間關(guān)系圖(a)流域蓄水容量-面積分配曲線(b)流域蓄水容量-面積分配曲線與降雨徑流關(guān)系主要計算公式如下：WM=a=若PE+R=若PE+R=2.2.3.4分水源計算按蓄滿產(chǎn)流模型計算出的總徑流量R中包括了各種徑流成分，由于各種水源的匯流規(guī)律和匯流速度不同，相應(yīng)采用的計算方法也不同。因此，必須進行水源劃分。本次使用的新安江模型采用三水源劃分。三水源的水源劃分結(jié)構(gòu)借鑒了山坡水文學(xué)的概念，去掉了FC，用自由水蓄水庫結(jié)構(gòu)解決水源劃分問題。自由水蓄水庫結(jié)構(gòu)如圖5所示。自由水蓄水庫結(jié)構(gòu)考慮了包氣帶的垂向調(diào)蓄作用。按蓄滿產(chǎn)流模型計算出總徑流量R，先進入自由水蓄水庫調(diào)蓄，再劃分水源。從圖可見，產(chǎn)流面積上自由水蓄水庫設(shè)置了兩個出口，一個為旁側(cè)出口，形成壤中流RI；另一個為向下出口，形成地下徑流RG。根據(jù)蓄滿產(chǎn)流的概念，只有在產(chǎn)流面積上才可能產(chǎn)生徑流，因為產(chǎn)流面積是變化圖5自由水蓄水庫結(jié)構(gòu)圖的，所以，自由水蓄水庫的底寬也是變化的。在圖中還設(shè)置了一個壤中流水庫，該水庫用于壤中流受調(diào)蓄作用大的流域，即將劃分出來的壤中流再進行一次調(diào)蓄計算。該水庫一般是不需要的，故在圖中用虛線表示。由于飽和坡面流的產(chǎn)流面積是不斷變化的，所以在產(chǎn)流面積FR上自由水蓄水容量分布是不均勻的。三水源水源劃分結(jié)構(gòu)是采用類似于流域蓄水容量—面積分配曲線的流域自由水蓄水容量—面積分配曲線來考慮流域內(nèi)自由水蓄水容量分布不均勻的問題。主要計算公式如下：流域平均自由水蓄水容量為：Sm時段初自由水蓄量S1S1相應(yīng)的縱坐標AU為：AU=當PE+AU<RS=當PE+AU≥RS=時段末的自由水蓄量S2S2相應(yīng)的壤中流RI為：RI=地下徑流RG為：RG=3結(jié)果與討論3.1產(chǎn)流計算根據(jù)VBA編程計算得出的數(shù)據(jù)，應(yīng)用oringe軟件繪制降雨及徑流的相關(guān)圖形如下。(a)81-83年徑流深過程圖(b)徑流深與降雨對比圖圖6組合一(a)81-83年徑流深過程圖(b)徑流深與降雨對比圖圖7組合二(a)81-83年徑流深過程圖(b)徑流深與降雨對比圖圖8組合三圖981-83年三種組合徑流深過程對比圖由圖6-9可知，三種參數(shù)組合情況下的徑流過程基本一致，變化較小，且徑流過程與降雨過程有明顯的對應(yīng)關(guān)系，徑流受降雨影響較大。3.2水源劃分利用VBA編程將徑流劃分為地表徑流、壤中流與地下徑流三種成分，應(yīng)用oringe軟件繪制的三種水源過程圖如圖10-12所示。圖1081-83年三水源過程圖(組合一)圖1181-83年三水源過程圖(組合二)圖1281-83年三水源過程圖(組合三)由圖10-12可知，三種參數(shù)組合情況下的地表徑流、壤中流及地下徑流過程基本一致。每種組合的三種水源成分占徑流的比例也基本相同。其中，地表徑流在集中降雨大時出現(xiàn)較高值，原因是流域已經(jīng)蓄滿從而產(chǎn)生大量地表徑流。壤中流與地下徑流受蓄水容量限制變化較為平緩，幾乎相等。3.3模型參數(shù)新安江模型參數(shù)可分蒸散發(fā)計算、產(chǎn)流計算、分水源計算和匯流計算四個層次，本次課設(shè)只進行到第三層次，相應(yīng)的參數(shù)如表4所示。表4新安江模型各層次參數(shù)表層次參數(shù)符號參數(shù)意義敏感程度參數(shù)值第一層次蒸散發(fā)計算KC流域蒸散發(fā)折算系數(shù)敏感0.90/0.95/1.00UM上層張力水容量（mm）敏感20/25/30LM下層張力水容量（mm）敏感40/50/60C深層蒸散發(fā)折算系數(shù)不敏感0.12/0.14/0.16第二層次產(chǎn)流計算WM流域平均張力水容量（mm）不敏感120/130/140B張力水蓄水容量曲線方次不敏感0.20/0.25/0.30IM不透水面積占全流域面積的比例不敏感第三層次水源劃分SM表層自由水蓄水容量（mm）敏感35/40/45EX表層自由水蓄水容量曲線方次不敏感1.0/1.2/1.3KG表層自由蓄水水庫對地下水的日出流系數(shù)敏感0.35KI表層自由蓄水水庫對壤中流的日出流系數(shù)敏感0.35分別進行三個層次的參數(shù)敏感性分析：1、蒸散發(fā)計算表5參數(shù)KC的敏感性分析KC0.900.951.00年徑流R/mm4432.74340.64250.1WM=130，WUM=25，WLM=40，WDM=65，WMM=156，B=0.2，C=0.12KC是影響產(chǎn)流量計算最為重要和敏感的參數(shù)，產(chǎn)流計算中KC控制著水量平衡，因此，對水量計算是最重要的。它主要反映流域平均高程與蒸發(fā)站高程之間差別的影響和蒸發(fā)皿蒸散發(fā)與陸地蒸散發(fā)間差別的影響。該參數(shù)值的改變，影響初始土壤含水量和雨期蒸發(fā)，導(dǎo)致計算產(chǎn)流量的改變。如表5所示，KC值增大，蒸發(fā)量增大，產(chǎn)流量減少；KC減小，蒸發(fā)量降低，產(chǎn)流量增大。表6參數(shù)WUM的敏感性分析WUM202530年徑流R/mm4341.84334.74321.6KC=0.95，WM=140，WLM=50，B=0.25，C=0.14WUM與作物根系土層厚度有關(guān)，其取值范圍為10-50mm。由表6可知，產(chǎn)流量R隨著上層張力水容量WUM值的增加而減小，且變化幅度明顯，較為敏感。表7參數(shù)WLM的敏感性分析WLM405060年徑流R/mm4432.74433.04443.2KC=0.90，WM=130，WUM=25，B=0.2，C=0.12WLM與包氣帶的土層結(jié)構(gòu)和物理特性有關(guān)，其值約在60-90mm（濕潤地區(qū)）之間。產(chǎn)流量R隨著下層張力水容量WLM值的增加而變大，其變化幅度較WUM更明顯，敏感性更強。表8參數(shù)C的敏感性分析C0.120.140.16年徑流R/mm4334.74334.74334.7KC=0.95，WM=140，WUM=25，WLM=50，B=0.25，C主要取決于流域內(nèi)深根植物的覆蓋面積，是反映深根植物作用的參數(shù)，決定深層散發(fā)。由表8可知，該值的變化對年徑流的影響很小，參數(shù)不敏感。2、產(chǎn)流計算表9參數(shù)WM的敏感性分析WM120130140年徑流R/mm4367.14341.64340.9KC=0.95，WUM=30，WLM=50，B=0.25，C=0.14流域平均張力水容量WM表示流域干旱程度和影響土壤水分變化的土層深度。WM在模型中不是很靈敏，如WM分別取值130和140時，年徑流變化不大，但WM值不能太小，可能會出現(xiàn)負值。表10參數(shù)B的敏感性分析B0.200.250.30年徑流R/mm4432.74435.44437.9KC=0.90，WM=130，WUM=25，WLM=40，C=0.12B值反映劃分單元流域張力蓄水分布的不均勻程度。B值越大表示流域越不均勻。其值的改變對年徑流的影響較小，參數(shù)敏感性不強。3、水源劃分表11參數(shù)SM的敏感性分析SM354045RS1111.3997.0902.8RI1520.51573.61617.4RG1520.51573.61617.4WM=130，EX=1.0，B=0.2，KI=0.35，KG=0.35SM反映表土蓄水能力，其值受降雨資料時段均化的影響明顯。當以日為計算時段長時，一般流域的SM值約為10～50mm，當所選取的計算時段長較小時，SM要增大，這個參數(shù)對地面徑流的多少起著決定性作用。如表11所示，SM增大，RS減小，RI與RG均增大，三者變化顯著。表12參數(shù)EX的敏感性分析EX1.01.21.3RS1111.31151.61168.7RI1520.51502.01494.2RG1520.51502.01494.2WM=130，SM=35，B=0.2，KI=0.35，KG=0.35EX反映流域自由水蓄水分布的不均勻程度，該值范圍一般為1.0～1.5，由表12可知，隨著EX的增大，RS逐漸增大，RI與RG均減小。3.4誤差分析影響流域降雨徑流過程的因素很多，三水源新安江模型的結(jié)構(gòu)與參數(shù)能夠反映濕潤地區(qū)降雨徑流過程的主要規(guī)律與特點，因而能獲得較好的精度。但是，模型本身以及模型計算中有許多概化，會造成誤差。造成遂川江流域年徑流量模擬誤差來源主要有以下幾方面：1、雨量站代表性的影響遂川江流域面積288.24km2，包括小夏、桐古、營盤墟、沙湖里、阡陌、淋洋、七嶺、滁州8個水文站的集水面積。對于多年平均降雨來講，基本上能控制降雨的空間分布。但對于不同年份不同時期，不同類型洪水，其差別較大。對于次洪水來說，在不同的氣象條件下計算出的面平均雨量也會有較大的誤差，面平均雨量的計量誤差是產(chǎn)流量誤差的主要影響因素。因此，在計算流域徑流過程中，要充分考慮流域的氣象條件和地質(zhì)的時空不均勻性。2、人類活動的影響隨著社會經(jīng)濟的快速發(fā)展，人類活動的影響加劇，流域內(nèi)先后修建了一些中小型水庫，這些中小型水庫大多數(shù)沒有固定調(diào)度原則和調(diào)洪方式，干旱季節(jié)或初汛的時候，流域內(nèi)的大小水利工程均需蓄水，會使實測的徑流量偏小。人類活動的影響也會造成流域下墊面的變化，下墊面的改變與不均勻分布將會產(chǎn)生一定的誤差。3、模型參數(shù)的影響模型參數(shù)是根據(jù)輸入，通過模型計算輸出，再將輸出過程與實測過程進行比較，用系統(tǒng)識別的方法作優(yōu)化調(diào)試的，因為需要率定的參數(shù)較多，有一部分參數(shù)是直接給定的，因此上述所率定的模型參數(shù)可能不是最優(yōu)。此外，用來率定參數(shù)的資料序列太短，資料的代表性不強，并且概念性水文模型本身精度也有限。3.5問題討論3.5.1流域劃分根據(jù)遂川江滁州水文站（東經(jīng)114o08′，北緯26o21′）以上控制流域水系及水文站雨量站網(wǎng)的分布特點，天然流域劃分為小夏、營盤墟、桐古、沙湖里、淋洋、阡陌、七嶺、滁州8塊單元面積，然后再在每一塊面積內(nèi)用泰森多邊形細分。這樣可以更符合降雨空間分布實際情況，解決面雨量計算誤差問題。3.5.2實時校正模型計算值與實測值（流量或水位）之間總是存在一定的誤差。造成兩者間誤差的因素很多，若針對每一個單一因素，它們是難于描述和預(yù)見的，一般采用實時校正模型來解決。實時校正模型的種類有很多，主要與所選擇的預(yù)報模型有關(guān)?？梢猿浞掷昧饔騼?nèi)水文（或水位）站的信息，對模型計算值進行實時校正。4參考文獻[1]劉新仁.流域水文模型的研究途徑[J].全國水文學(xué)術(shù)討論會論文集，2004，189-192.[2]趙人俊.流域水文模型-新安江模型與陜北模型[M].北京：水利電力出版社，1984.[3]趙人俊.新安江模型參數(shù)的分析[J].水文，1998（6）：2-9.[4]包為民.水文預(yù)報（第4版）[M].北京.中國水利水電出版社，2016.[5]《水文預(yù)報課程設(shè)計任務(wù)書指示書》.揚州大學(xué)水利與能源動力工程學(xué)院.5附錄：程序代碼5.1蓄滿產(chǎn)流計算Subtest1()DimPJAsDouble'流域日平均降雨量，mmDimE0AsDouble'流域水面蒸發(fā)，mmDimEPAsDouble'流域蒸發(fā)能力，mmDimEUAsDouble'土壤上層蒸發(fā)量，mmDimELAsDouble'土壤的下層蒸發(fā)量，mmDimEDAsDouble'土壤的深層蒸發(fā)量，mmDimEAsDouble'土壤總蒸發(fā)量，mmDimPEAsDouble'扣除雨期蒸發(fā)后的降雨量，mmDimWUAsDouble'上層土壤含水量，mmDimWLAsDouble'下層土壤含水量，mmDimWDAsDouble'深層土壤含水量，mmDimWAsDouble'總含水量，mmDimaAsDouble'流域初始含水量最大值，mmDimRAsDouble'流域蓄滿產(chǎn)流量，mmDimKCAsDouble'蒸發(fā)折算系數(shù)DimWMAsDouble'流域平均蓄水容量，mmDimWUMAsDouble'上層土壤含水容量，mmDimWLMAsDouble'下層土壤含水容量，mmDimWDMAsDouble'深層土壤含水容量，mmDimWMMAsDouble'流域蓄水容量最大值，mmDimBAsDouble'蓄水容量曲線指數(shù)，mmDimCAsDouble'蒸發(fā)擴散系數(shù)DimxAsIntegerForx=3To1097Step1PJ=Range("b"&x)E0=Range("c"&x)EP=Range("d"&x)EU=Range("e"&x)EL=Range("f"&x)ED=Range("g"&x)E=Range("h"&x)PE=Range("i"&x)WU=Range("j"&x)WL=Range("k"&x)WD=Range("l"&x)W=Range("m"&x)a=Range("n"&x)R=Range("o"&x)KC=0.9WM=130WUM=25WLM=40WDM=65WMM=156B=0.2C=0.12'三層蒸發(fā)模式計算IfRange("j"&x)+Range("b"&x)>=Range("d"&x)Then'第一種情況Range("e"&x)=Range("d"&x)Range("f"&x)=0Range("g"&x)=0ElseIfRange("j"&x)+Range("b"&x)<Range("d"&x)AndRange("k"&x)>=C*WLMThen'第二種情況Range("e"&x)=Range("j"&x)+Range("b"&x)Range("f"&x)=(Range("d"&x)-Range("e"&x))*Range("k"&x)/WLMRange("g"&x)=0ElseIfRange("j"&x)+Range("b"&x)<Range("d"&x)AndC*(Range("d"&x)-Range("e"&x))<=Range("k"&x)AndRange("k"&x)<C*WLMThen'第三種情況Range("e"&x)=Range("j"&x)+Range("b"&x)Range("f"&x)=C*(Range("d"&x)-Range("e"&x))Range("g"&x)=0ElseIfRange("j"&x)+Range("b"&x)<Range("d"&x)AndRange("k"&x)<C*(Range("d"&x)-Range("e"&x))Then'第四種情況Range("e"&x)=Range("j"&x)+Range("b"&x)Range("f"&x)=Range("k"&x)Range("g"&x)=C*(Range("d"&x)-Range("e"&x))-Range("f"&x)EndIfRange("h"&x)=Range("e"&x)+Range("f"&x)+Range("g"&x)'土壤總蒸發(fā)量計算Range("i"&x)=Range("b"&x)-Range("h"&x)'扣除雨期蒸發(fā)后的降雨量計算Range("m"&x)=Range("j"&x)+Range("k"&x)+Range("l"&x)'總含水量計算Range("n"&x)=WMM*(1-(1-Range("m"&x)/WM)^(1/(1+B)))'流域初始含水量最大值計算'流域蓄滿產(chǎn)流計算IfRange("i"&x)>0Then'產(chǎn)流SelectCaseRange("n"&x)+Range("i"&x).ValueCaseIs<=WMM'全流域蓄滿前產(chǎn)流量計算Range("o"&x)=Range("i"&x)+Range("m"&x)-WM+WM*(1-(Range("i"&x)+Range("n"&x))/WMM)^(B+1)CaseElse'全流域蓄滿后產(chǎn)流量計算Range("o"&x)=Range("i"&x)-(WM-Range("m"&x))EndSelectIfRange("j"&x)+Range("i"&x)-Range("o"&x)<=WUMThen'上層土壤未蓄滿時Range("j"&x+1)=Range("j"&x)+Range("i"&x)-Range("o"&x)Range("k"&x+1)=Range("k"&x)Range("l"&x+1)=Range("l"&x)ElseIfRange("j"&x)+Range("i"&x)-Range("o"&x)>WUMAndRange("k"&x)+Range("j"&x)+Range("i"&x)-Range("o"&x)-WUM<=WLMThen'上層土壤蓄滿、下層土壤未蓄滿時Range("j"&x+1)=WUMRange("k"&x+1)=Range("k"&x)+Range("j"&x)+Range("i"&x)-Range("o"&x)-WUMRange("l"&x+1)=Range("l"&x)ElseIfRange("k"&x)+Range("j"&x)+Range("i"&x)-Range("o"&x)-WUM>=WLMAndRange("l"&x)+Range("k"&x)+Range("j"&x)+Range("i"&x)-Range("o"&x)-WDM<WDMThen'下層土壤蓄滿、深層土壤未蓄滿時Range("j"&x+1)=WUMRange("k"&x+1)=WLMRange("l"&x+1)=Range("l"&x)+Range("k"&x)+Range("j"&x)+Range("i"&x)-Range("o"&x)-WDMElseIfRange("l"&x)+Range("k"&x)+Range("j"&x)+Range("i"&x)-Range("o"&x)-WDM>=WDMThen'深層土壤蓄滿時Range("j"&x+1)=WUMRange("k"&x+1)=WLMRange("l"&x+1)=WDMEndIfElseIfRange("i"&x)<=0Then'不產(chǎn)流IfRange("j"&x)+Range("i"&x)>0Then'上層土壤含水量未蒸發(fā)完時Range("j"&x+1)=Range("j"&x)+Range("i"&x)Range("k"&x+1)=Range("k"&x)Range("l"&x+1)=Range("l"&x)ElseIfRange("k"&x)+Range("j"&x)+Range("i"&x)>0Then'上層土壤含水量蒸發(fā)完、下層未蒸發(fā)完時Range("j"&x+1)=0Range("k"&x+1)=Range("k"&x)+Range("j"&x)+Range("i"&x)Range("l"&x+1)=Range("l"&x)ElseIfRange("l"&x)+Range("k"&x)+Range("j"&x)+Range("i"&x)>0Then'下層土壤含水量蒸發(fā)完、深層未蒸發(fā)完時Range("j"&x+1)=0Range("k"&x+1)=0Range("l"&x+1)=Range("l"&x)+Range("k"&x)+Range("j"&x)+Range("i"&x)ElseIfRange("l"&x)<=0Then'深層土壤含水量蒸發(fā)完時Range("j"&x+1)=0Range("k"&x+1)=0Range("l"&x+1)=0EndIfEndIfNextxEndSub5.2水源劃分Subtest2()DimPEAsDouble'扣除雨期蒸發(fā)后的降雨量，mmDimRAsDouble'流域蓄滿產(chǎn)流量，mmDimFRAsDouble'產(chǎn)流面積DimAUAsDouble'時段初始自由水蓄量相應(yīng)的縱坐標DimSAsDouble'自由水蓄量，mmDimRSAsDouble'地面徑流，mmDimRIAsDouble'壤中流，mmDimRGAsDouble'地下徑流，mmDimWAsLong'總含水量，mmDimWMAsDouble'流域平均蓄水容量，mmDimSMAsDouble'深度，mmDimEXAsDouble'自由水蓄量分布曲線指數(shù)DimBAsDouble'產(chǎn)流面積計算指數(shù)DimSmmAsDouble'流域最大自由蓄水量，mmDimKIAsDouble'壤中流的出流系數(shù)DimKGAsDouble'地下徑流的出流系數(shù)DimxAsIntegerForx=4To1098PE=Range("b"&x)R=Range("c"&x)FR=Range("d"&x)AU=Range("e"&x)S=Range("f"&x)RS=Range("g"&x)RI=Range("h"&x)RG=Range("i"&x)W=Range("j"&x)WM=130SM=35EX=1B=0.2Smm=70KI=0.35KG=0.35'三水源劃分SelectCaseRange("b"&x).ValueCaseIs<=0Range("d"&x)=1-(1-Range("j"&x)/WM)^(B/(1+B))Range("g"&x)=0Range("h"&x)=Range("f"&x)*KI*Range("d"&x)Range("i"&x)=Range("f"&x)*KG*Range("d"&x)Range("f"&x+1)=Range("f"&x)*(1-KI-KG)CaseElseRange("d"&x)=Range("c"&x)/Range("b"&x)Range("e"&x)=Smm*(1-(1-Range("f"&x)*Range("d"&x-1)/Range("d"&x)/SM)^(1/(1+EX)))SelectCaseRange("b"&x).Value+Range("e"&x).ValueCaseIs<SmmRange("g"&x)=Range("d"&x)*(Range("b"&x)+(Range("f"&x)*Range("d"&x-1)/Range("d"&x))-SM+SM*(1-((Range("b"&x)+Range("e"&x))/Smm)^(1+EX)))Range("h"&x)=(SM-SM*(1-(Range("b"&x)+Range("e"&x))/Smm)^(1+EX))*KI*Range("d"&x)Range("i"&x)=(SM-SM*(1-(Range("b"&x)+Range("e"&x))/Smm)^(1+EX))*KG*Range("d"&x)Range("f"&x+1)=(1-KI-KG)*(SM-SM*