海綿城市綠色屋頂降雨徑流模型構(gòu)建及建議研究**

李鵬程，郝 放，趙建利

（1.河北水利電力學(xué)院滄州市遙感與智慧水利技術(shù)創(chuàng)新中心，河北 滄州 061001；2.中水北方勘測設(shè)計(jì)研究有限責(zé)任公司，天津 300222；3.河海大學(xué)水利水電學(xué)院，江蘇 南京 210013；4.北京市水科學(xué)技術(shù)研究院，北京 100048）

城市水生態(tài)文明體系構(gòu)建要求對水資源開展有效保護(hù)及合理利用，綜合整治城市發(fā)展過程的水環(huán)境問題，是城市可持續(xù)發(fā)展的關(guān)鍵因素之一。海綿城市是通過構(gòu)建低影響開發(fā)措施，包括綠色屋頂、植草溝、生物滯留池等，凈化并截留降雨，對城市雨洪進(jìn)行調(diào)控，對城市水生態(tài)文明體系構(gòu)建具有積極促進(jìn)作用。

對海綿城市低影響措施的研究主要通過理論分析、水文模型及實(shí)驗(yàn)分析進(jìn)行，其中，水文模型是通過一定的數(shù)學(xué)物理模型及邊界條件的設(shè)定，模擬降雨產(chǎn)流過程[1-5]，將實(shí)際降雨過程概化處理，因此模型構(gòu)建過程中參數(shù)的選取直接影響模型模擬的精度。本文通過建立綠色屋頂降雨產(chǎn)流模型，并設(shè)定初始條件和邊界條件，進(jìn)行模型率定及參數(shù)敏感性分析，以對提高模擬降雨產(chǎn)流精度。

1 模型構(gòu)建及研究方法

1.1 降雨徑流模型

運(yùn)用HYDRUS-1D計(jì)算模擬綠色屋頂削減降雨徑流過程，用它可以解算在不同邊界條件制約下的數(shù)學(xué)模型，模型能夠模擬溶質(zhì)遷移以及多孔介質(zhì)飽和-非飽和滲流區(qū)水、熱過程，全面融入植物根系吸收、熱運(yùn)動、水分運(yùn)動、溶質(zhì)運(yùn)移以及土壤持蓄作用[6-7]，可不受限制地選擇大氣邊界、滲水邊界恒定水頭、排水溝、通量邊界、自由排水邊界等，對構(gòu)建的計(jì)算模型方程進(jìn)行運(yùn)算求解，并運(yùn)用Crank-Nicholson有限差分法來進(jìn)行時間離散計(jì)算，運(yùn)用Galerkin有限元法進(jìn)行空間離散，此模型可用于模擬水在土壤中的移動過程。

1.2 土壤水運(yùn)動方程

土壤水運(yùn)動方程：局部飽和的多孔剛性介質(zhì)中的一維土壤水運(yùn)動，若假設(shè)空氣對水流過程影響不大且忽略熱力梯度的影響，則可用Richard方程進(jìn)行描述。

其中，h為壓力水頭[L]，θ為體積含水率[L3L-3]，t為時間[T]，x為空間坐標(biāo)[L]，S為SinkTeam[L3L-3T-1]，α為水流方向與豎直方向的夾角。非飽和水力傳導(dǎo)度函數(shù)：

式中，KS為飽和水力傳導(dǎo)度，Kr為相對水力傳導(dǎo)度。

1.3 非飽和土壤水力特性

1.3.1 模型原理

非飽和土壤的水力特性，土壤含水量θ(h)與水力傳導(dǎo)度K(h)與壓力水頭高度呈非線性關(guān)系。HYDRUS允許使用3種不同的水力特性模型[Brooks and Corey 1964；van Genuchtem 1980；Vogel and Cislerova 1988]。

本文模型采用van Genuchtem模型。

1.3.2 van Genuchtem模型

由統(tǒng)計(jì)上的氣孔分布模型Mualem獲得根據(jù)持水參數(shù)所獲得的非飽和水力傳導(dǎo)度函數(shù)：

上述方程包括5個獨(dú)立參數(shù)θr，θS，α，n，K；其中，水力傳導(dǎo)度方程的氣孔連接系數(shù)l對多種土壤估計(jì)得0.5。

以上這些參數(shù)，以現(xiàn)場取土在室內(nèi)測得的數(shù)據(jù)及HYDRUS提供的基本土壤數(shù)據(jù)中的數(shù)值作為初始值，根據(jù)實(shí)際觀測的土壤水分運(yùn)動過程經(jīng)過率定獲得。

1.4 初始條件和邊界條件

在模擬土壤水運(yùn)動過程中需設(shè)定土壤剖面的初始條件，設(shè)定初始含水率值。含水率：θ(z,t)=θo(z)，t=0。

在土壤剖面頂部（z=L）或底部（z=0）必須指定以下邊界條件之一：

式中，h0[L]與q0[LT-1]為邊界上壓力水頭與土壤水流的規(guī)定值。

HYDRUS-1D中又將上述邊界條件進(jìn)行了細(xì)化，其中，上邊界條件包括定通量邊界、土壤表層大氣邊界、定水頭邊界、土壤表層徑流大氣邊界、變通量邊界和變水頭邊界；下邊界條件包括定通量邊界、變通量邊界、定水頭邊界、變水頭邊界、滲漏面邊界、自由排水邊界、水平排水邊界條件和深層排水邊界。本研究中模型采用自由排水邊界和滲漏面邊界。

2 模型率定及參數(shù)敏感性分析

分別選取不同綠色屋頂在不同降雨條件下所模擬得出的出流值與實(shí)測值進(jìn)行對比，完成對模型的率定和驗(yàn)證。

檢驗(yàn)率定和驗(yàn)證是否達(dá)到要求，通過以下3個統(tǒng)計(jì)量評定：1）均方誤差RMSE；2）相關(guān)系數(shù)R2；3）平均相對誤差MRE。

式中：N是觀測值的個數(shù)，Oi表示第i個觀測值，Pi表示第i個觀測值的模擬值。

模型在模擬土壤水運(yùn)動過程中的輸入數(shù)據(jù)主要包括時間信息、幾何信息、土壤水力特征參數(shù)、初始條件、邊界條件以及土壤剖面信息等[8-9]。綠色屋頂?shù)哪M時間依據(jù)試驗(yàn)的降雨和產(chǎn)流時間來確定[10-11]；模擬的土壤層厚度為填料層厚度和種植層厚度之和，且不同的綠色屋頂有不同的厚度，分為10 cm和20 cm。輸入各填料層對應(yīng)的土壤水力特征參數(shù)；模擬的初始條件為初始時刻土壤剖面的含水率值；上邊界條件均為土壤表層的大氣邊界條件，由于模擬時間較短（小于5 h），不考慮土壤蒸發(fā)和作物蒸騰作用，只輸入降雨數(shù)據(jù)，其中積水深度為蓄水層厚度（15 cm）；下邊界條件均選擇滲漏面邊界和自由出流邊界。在模型率定過程中，可以對土壤水力特征參數(shù)和土壤初始含水率進(jìn)行調(diào)整，以取得更好的率定效果。

由于此次降雨受前期降雨影響比較小，所以初始含水率設(shè)置在0.35左右。土壤水分特征曲線值以試驗(yàn)測量值為初始值進(jìn)行率定。如表1所示。

表1 初始土壤水分特征參數(shù)表

2.1 初始含水率進(jìn)行率定

由于基質(zhì)的上表層含水率和下底層含水率有一定差別，首先對含水率的差值進(jìn)行率定，此處選取六組模擬數(shù)據(jù)進(jìn)行研究觀察，模擬結(jié)果如圖1所示。

由圖1可知，當(dāng)含水率介于0.32~0.36之間時（其平均含水率為0.34），與含水率介于0.25~0.43之間時候（平均含水率為0.34）模擬徑流過程基本一致；同樣，當(dāng)上表面與下表面的平均含水率為0.345和0.35時，兩組徑流過程線也會重合。因此，可以看出基質(zhì)上層含水率與下層含水率設(shè)置不同數(shù)值時，主要為平均含水率對模擬結(jié)果產(chǎn)生影響。

圖1 初始含水率率定圖

2.2 土壤水分特征參數(shù)敏感性分析

對于飽和含水率，取飽和含水率為0.54、0.548（試驗(yàn)值）和0.556進(jìn)行分析，模擬結(jié)果如圖2所示。

圖2 飽和含水率敏感分析圖

三種模擬值由小到大累計(jì)徑流深對應(yīng)值分別為13.1 mm、12.4 mm、11.7 mm，而實(shí)測值為4 mm。飽和含水率分別取了0.54、0.548（試驗(yàn)值）和0.556，差值為0.008，該差值對于峰值而言為敏感性參數(shù)。三者產(chǎn)流的起始過程線基本一致，在峰值附近的產(chǎn)流過程差異較大，當(dāng)飽和含水率增大時產(chǎn)流峰值降低，產(chǎn)流速率減慢；當(dāng)飽和含水率減小時，產(chǎn)流峰值增加，產(chǎn)流速率加快。當(dāng)飽和含水率為0.54時模擬值與實(shí)際值較接近。

殘留含水率差值為0.004徑流過程線圖，如圖3所示。

三種模擬值由小到大累計(jì)徑流深分別為13.1 mm、13.1 mm、13 mm，無太大變化。由圖3可知，三種模擬值的過程基本重合，說明就差值0.004而言，對產(chǎn)流基本無影響。故增大差值進(jìn)一步比較，此次選擇差值0.01和0.03，模擬0.088和0.048兩個產(chǎn)流過程，其結(jié)果如圖4所示。

圖3 殘留含水率差值為0.004徑流過程線圖

圖4 殘留含水率差值為0.004、0.01、0.03產(chǎn)流過程對比圖

以上五種模擬值只有當(dāng)殘留含水率為0.048時累計(jì)徑流深變化最大，為13.4 mm。通過圖4殘留含水率差值分別為0.004、0.01、0.03產(chǎn)流過程對比圖可知，過程線基本重合，殘留含水率對降雨過程的產(chǎn)流并無有效影響，故以下的率定過程暫不予以考慮。

固定除飽和導(dǎo)水率之外的土壤水分參數(shù)值，選擇初始值0.059（試驗(yàn)值），差值初定為0.005，即模擬0.049、0.054、0.059和0.064四個產(chǎn)流過程線。其結(jié)果如圖5所示。

圖5 飽和導(dǎo)水率差值為0.005產(chǎn)流過程對比圖

四種模擬值由小到大累計(jì)徑流深分別為12.5 mm、12.8 mm、13.1 mm、13.4 mm。由圖5可知，飽和導(dǎo)水率對模型的產(chǎn)流時間以及初期、后期的產(chǎn)流過程有較小的影響，對峰值有一定的影響，飽和導(dǎo)水率變小，徑流率最大值也隨之變小，變化值為0.000 5 cm/min左右。飽和導(dǎo)水率越小，模擬值的總徑流深越接近實(shí)測值。

四組模擬數(shù)據(jù)在峰值附近有小幅變化，但是不明顯。由于此參數(shù)對該降雨模型有一定影響，故模擬差值初定為0.001產(chǎn)流過程進(jìn)行比照分析。其過程對比如圖6所示。

圖6 α差值為0.001和0.005產(chǎn)流過程對比

α值為0.003 6~0.009 6時，累計(jì)徑流深均為12.5 mm，可見α值對徑流總量影響較小。由圖6產(chǎn)流過程線可知，α值增大，峰現(xiàn)時間增長，體態(tài)變“胖”即產(chǎn)流稍微提前，峰值略有增加。

綜上，通過參數(shù)的敏感性分析，做匯總表，結(jié)果如表2所示。

表2 人工改良土特征參數(shù)敏感性分析表

3 小結(jié)

通過HYDRUS-1D對綠色屋頂不同降雨條件的模擬，觀察徑流深、徑流峰值及徑流過程的變化，對不同參數(shù)下的降雨產(chǎn)流過程進(jìn)行對比分析可以得出：

模擬參數(shù)中飽和含水率、N值對于綠色屋頂?shù)膹搅魈匦杂绊懽顬槊黠@，是決定產(chǎn)流時間和產(chǎn)流量的主要因素，是模擬精度提高的關(guān)鍵參數(shù)。殘留含水率對各種雨強(qiáng)的模擬基本沒有影響。

由于綠色屋頂?shù)慕Y(jié)構(gòu)組成比較復(fù)雜，影響因素眾多，模擬值與實(shí)測值有一定差值，本模擬只是產(chǎn)流過程線基本一致，這取決于模型的理論算法，有待于探索更合適的模型算法。

研究成果為海綿城市綠色屋頂降雨徑流模型構(gòu)建提供了針對性參考，為提高模型模擬精度提供技術(shù)支撐，后續(xù)其他海綿城市低影響開發(fā)措施降雨徑流的模型構(gòu)建、模擬及參數(shù)敏感性分析可參考本文的研究方法。