• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    基于氣象水文水動(dòng)力耦合模型的流域尺度洪水預(yù)報(bào)方法: 以北江流域?yàn)槔?/h1>
    2024-05-20 07:33:24韓赟希張昕王自法周良辰朱龍祥
    科學(xué)技術(shù)與工程 2024年11期
    關(guān)鍵詞:北江水文徑流

    韓赟希, 張昕, 王自法, 周良辰, 朱龍祥

    (1.河南大學(xué)土木建筑學(xué)院, 開封 475004; 2.深圳防災(zāi)減災(zāi)技術(shù)研究院, 深圳 518003; 3.中震科建(廣東)防災(zāi)減災(zāi)研究院, 韶關(guān) 512026; 4.中國地震局工程力學(xué)研究所, 哈爾濱 150080)

    近年來,全球氣候變化引發(fā)了極端暴雨事件的頻繁發(fā)生和強(qiáng)度增加[1],這導(dǎo)致洪水災(zāi)害的風(fēng)險(xiǎn)日益突出,在這種背景下,有效的洪水預(yù)報(bào)[2]成為保障人民生命財(cái)產(chǎn)安全的重要舉措。因此,發(fā)展和應(yīng)用先進(jìn)的洪水預(yù)報(bào)技術(shù)和方法,能夠更準(zhǔn)確地預(yù)測(cè)洪水的發(fā)生和發(fā)展趨勢(shì),最大限度地減少洪水災(zāi)害對(duì)人們生活和社會(huì)經(jīng)濟(jì)的影響。

    洪水是一種受多種自然物理和人為因素影響的復(fù)雜現(xiàn)象[3],其發(fā)生過程通常涉及3個(gè)主要方面:降雨、徑流和淹沒。在實(shí)際的洪水研究中,降雨通常使用大氣數(shù)值模式[4]進(jìn)行模擬,如WRF (weather research and forecasting model) 模式[5-7];徑流通常由水文模型模擬,如新安江模型[8]、HEC-HMS (hydrologic engineering center-hydrologic modeling system) 模型[9]等;洪水淹沒這一過程通常采用水動(dòng)力模型模擬,如MIKE 21[10]、HEC-RAS (river analysis system)[11]等。這種將復(fù)雜的流域系統(tǒng)拆分并分別處理的方法,雖然有助于深入研究每個(gè)環(huán)節(jié)的影響機(jī)制,但難以完整地考慮氣象、水文、水動(dòng)力之間的緊密聯(lián)系。

    大氣數(shù)值模式的輸出產(chǎn)品通常以較大尺度分辨率為主,而流域水文過程具有更為復(fù)雜的空間變化特征,直接將其輸出應(yīng)用于精細(xì)的流域水文預(yù)報(bào)存在一定困難和限制。近些年來,出現(xiàn)了兩者耦合建模方法(即陸氣耦合水文模型),例如,Cerbelaud等[12]采用WRF-hydro模式對(duì)新喀里多尼亞的6個(gè)流域進(jìn)行了為期兩年的模擬,取得了不錯(cuò)的長期徑流模擬效果,證明WRF-Hydro是一個(gè)合適的框架,可以更好地了解水文過程,并高效分析新喀里多尼亞的土壤和氣候在未來的變化。劉洪波等[13]將WRF-Hydro應(yīng)用于半濕潤的陳河流域和濕潤的屯溪流域,并結(jié)合新安江模型以探索其在中小流域中的徑流模擬能力,結(jié)果表明WRF-Hydro模式在濕潤地區(qū)有較好的應(yīng)用潛力;而應(yīng)用在半濕潤地區(qū),可能需要更高質(zhì)量的輸入數(shù)據(jù)。另一方面,流域和河流之間往往存在著天然水力聯(lián)系[14],為了保持這種聯(lián)系,水文水動(dòng)力耦合模型的方法將成為流域洪水的精細(xì)化預(yù)測(cè)及災(zāi)情的快速評(píng)估[15]的重要基礎(chǔ),近年來國內(nèi)外多位學(xué)者展開了該方面的研究。例如,Fleischmann等[16]以亞馬遜流域?yàn)檠芯繉?duì)象,通過耦合MGB (modelo de grandes bacias) 水文模型和CaMa-Flood水動(dòng)力模型,模擬了該流域在1980—2014年期間的逐日淹沒情況,并結(jié)合實(shí)測(cè)影像資料,首次分析了40年來洪泛區(qū)的淹沒范圍與全流域降雨變化的關(guān)系。沈澤宇等[17]以溫德河流域?yàn)檠芯繉?duì)象,首先建立HEC-HMS水文模型并進(jìn)行參數(shù)優(yōu)化,之后建立HEC-RAS二維水動(dòng)力學(xué)模型對(duì)重要子流域進(jìn)行淹沒模擬,表明所建水文水動(dòng)力模耦合型,對(duì)具有復(fù)雜水文、水力條件流域的洪水預(yù)報(bào)準(zhǔn)確率較高。付曉花等[18]以淮河流域中下游結(jié)合部的洪澤湖周邊滯洪區(qū)為研究對(duì)象,基于松耦合技術(shù),將氣候模型與水文模型、水文模型與水動(dòng)力模型進(jìn)行耦合,實(shí)現(xiàn)了流域、河網(wǎng)、蓄滯洪區(qū)和湖泊復(fù)雜系統(tǒng)氣候水文水動(dòng)力單向耦合,并對(duì)氣候變化影響下的復(fù)雜河網(wǎng)地區(qū)水流演進(jìn)進(jìn)行了模擬。盡管許多研究都得到了不錯(cuò)的效果,但其應(yīng)用范圍仍然受到限制,僅集中于模擬單個(gè)子流域或特定區(qū)域的淹沒情況[19],難以預(yù)測(cè)整個(gè)流域的洪水淹沒情況。在全球氣候變化,極端洪水事件頻發(fā)的背景下,有必要開展陸氣耦合水文模型和大尺度水動(dòng)力模型的耦合預(yù)報(bào)研究,進(jìn)而模擬整個(gè)流域的徑流和洪水淹沒過程。

    在前人研究的基礎(chǔ)之上,針對(duì)其中的不足之處,將北江流域作為研究對(duì)象,采用高分辨率的氣象模式輸出產(chǎn)品[20-21]、陸氣耦合水文模式WRF-Hydro[22-23]、大尺度水動(dòng)力模型CaMa-Flood[24-26],構(gòu)建基于降雨徑流和洪水淹沒的綜合洪水預(yù)報(bào)模型,模擬洪水過程、洪峰流量、洪峰抵達(dá)時(shí)間、洪泛區(qū)的淹沒深度和淹沒范圍等要素,以期為保障流域防洪安全和人員及財(cái)產(chǎn)安全提供重要的技術(shù)支撐。

    1 研究區(qū)域概況與數(shù)據(jù)收集

    1.1 研究區(qū)域概況

    北江為珠江流域的第二大水系,發(fā)源于江西省信豐縣石碣大茅坑,經(jīng)大余縣進(jìn)入廣東,自東北往西南穿山越嶺,流經(jīng)南雄、始興、曲江等市(縣),至韶關(guān)市沙洲尾與支流武江匯合,始稱北江;再自北向南流經(jīng)英德、清新、清遠(yuǎn)至三水河口,在思賢滘與西江相通,注入珠江三角洲網(wǎng)河區(qū)。北江干流至三水區(qū)思賢滘全長4 600 km(廣東省境內(nèi)4 500 km),平均坡降0.26%,流域集雨面積46 710 km2,其中廣東省境內(nèi)集雨面積達(dá)42 930 km2,其余位于湖南、江西等省境內(nèi)。北江干流總比降平緩,洪水漲快退慢,持續(xù)時(shí)間長。上游高山峽谷眾多,集水面積超過1 000 km2的支流有墨江、錦江、武江、南花溪、南水、滃江、煙嶺河、連江、青蓮水、潖江、濱江、綏江、鳳崗河等13條,其中一級(jí)支流9條,按葉脈狀排列,從東西兩側(cè)匯入干流。

    北江流域背靠南嶺山脈,正處在山脈的向風(fēng)坡,加之河流水系呈闊葉脈狀分布,洪水匯流集中迅猛,洪水暴澆暴落,峰高量不大,范圍廣但歷時(shí)長,洪水過程大都呈單峰型或雙峰型,復(fù)峰型的洪水過程較少,每年汛期發(fā)生洪水3~4次。一般4—7月份出現(xiàn)洪水的機(jī)會(huì)多,為北江的主汛期。北江流域的高程和水系分布如圖1所示。

    圖1 北江流域高程和水系分布Fig.1 Elevation and water system distribution in Beijiang Basin

    1.2 數(shù)據(jù)收集

    收集的數(shù)據(jù)時(shí)間跨度覆蓋了整個(gè)模擬期,包括全球陸面數(shù)據(jù)同化系統(tǒng) (global land data assimilation system, GLDAS) 中的必要?dú)庀髷?shù)據(jù),其水平分辨率為0.25°,時(shí)間分辨率為3 h;中國區(qū)域地面氣象要素驅(qū)動(dòng)數(shù)據(jù)集(China meteorological forcing dataset, CMFD)中的降水?dāng)?shù)據(jù),其水平分辨率為0.1°,時(shí)間分辨率為3 h;北江流域石角水文站的實(shí)測(cè)徑流數(shù)據(jù),由于收集到的水文站徑流數(shù)據(jù)在時(shí)間上具有不連續(xù)性,因此在時(shí)間上利用線性插值方法處理,最終得到石角水文站1 h的徑流數(shù)據(jù)。同時(shí),為了使驅(qū)動(dòng)數(shù)據(jù)與陸面模式的分辨率保持一致,采用雙線性插值方法使兩套驅(qū)動(dòng)數(shù)據(jù)插值到3 km分辨率。模型的靜態(tài)地理數(shù)據(jù)來自WRF模式前處理系統(tǒng)(weather research and forecasting preprocessing system, WPS)數(shù)據(jù)庫,各嵌套層區(qū)域的靜態(tài)地表數(shù)據(jù)分辨率均為30″,驅(qū)動(dòng)水文模型的初始場(chǎng)資料采用NCEP FNL(Final)再分析資料,高程數(shù)據(jù)來自HydroSHEDS數(shù)據(jù)集,CaMa-Flood模型的下墊面資料由模型內(nèi)置。

    2 研究方法

    2.1 研究技術(shù)路線

    WRF-Hydro(weather research and forecasting model hydrological modeling system)模式和CaMa-Flood模型的耦合屬于松散耦合[27],即WRF-Hydro的輸出需要經(jīng)過一系列處理之后再傳遞給CaMa-Flood,耦合模型的構(gòu)建流程如圖2所示,其具體包括:①收集建模所需的地形地貌、土地利用、再分析降雨和水文站實(shí)測(cè)徑流等數(shù)據(jù)并進(jìn)行預(yù)處理工作;②采用ArcGIS軟件,根據(jù)數(shù)字高程模型(digital elevation model, DEM)建立河網(wǎng),再根據(jù)流域出口位置對(duì)流域進(jìn)行劃分;③建立WRF-Hydro水文模型,選取場(chǎng)次洪水,完成水文模型的參數(shù)率定,得到優(yōu)化后的研究區(qū)域的水文模擬結(jié)果,對(duì)洪水的量級(jí)做出相對(duì)準(zhǔn)確的評(píng)價(jià);④針對(duì)耦合的需求,調(diào)整水文模型的一些選項(xiàng),在陸面模式的計(jì)算結(jié)果中選取地表徑流和地下徑流之和作為水動(dòng)力模型的邊界條件,其余條件均在CaMa-Flood模型中構(gòu)建,經(jīng)過該耦合過程可以快速地把洪水淹沒過程進(jìn)行系統(tǒng)性的模擬。

    圖2 研究技術(shù)路線圖Fig.2 Research technology roadmap

    2.2 WRF模式

    WRF模式是一種由美國國家大氣研究中心(national center for atmospheric research, NCAR)主導(dǎo)開發(fā)的中尺度氣象模式。該模式由兩種動(dòng)力核組成,分別是ARW(advanced research WRF)和NMM(non-hydrostatic mesoscale model),分別用于科學(xué)研究和業(yè)務(wù)應(yīng)用。

    選擇使用ARW動(dòng)力核,并以FNL數(shù)據(jù)作為輸入來驅(qū)動(dòng)WRF模式,以獲得WRF-Hydro模式所需的初始邊界場(chǎng)數(shù)據(jù)。此外,為了得到更高的模擬精度和空間分辨率的初始條件,將WRF模式設(shè)置為三層嵌套結(jié)構(gòu)。這三層分別對(duì)應(yīng)的水平分辨率為27、9、3 km,第三層的區(qū)域范圍即為所關(guān)注的水文模型區(qū)域。

    2.3 WRF-Hydro模式

    WRF-Hydro模式是新一代的分布式氣象水文模式,由美國國家大氣研究中心于2015年提出。該模式是在WRF模式的基礎(chǔ)上進(jìn)行擴(kuò)展,通過整合大氣動(dòng)力學(xué)、陸面過程和水文過程,實(shí)現(xiàn)了對(duì)流域水文過程的綜合建模和模擬。設(shè)計(jì)目標(biāo)是提供高分辨率的水文預(yù)報(bào)和水資源管理能力,以更準(zhǔn)確地模擬地表和地下水的水文過程。該模式主要由陸面模型、壤中匯流、地表匯流、地下水模型及河道匯流模塊構(gòu)成。

    本文中采用了WRF-Hydro V5.1.1版本,陸面模式選擇Noah-MP模式,采用默認(rèn)的2 m土壤柱,并分為4個(gè)土壤層,各土壤層的厚度分別為10、30、60、100 cm,其粗網(wǎng)格分辨率與WRF模式中的最內(nèi)層分辨率相同,為3 km。而匯流模塊則采用更細(xì)的網(wǎng)格分辨率,為500 m。為了處理不同網(wǎng)格分辨率之間的信息傳遞,使用了子網(wǎng)格聚合因子為6。模擬的時(shí)間步長設(shè)定為10 s,以更精細(xì)地捕捉流域的水文過程。地表匯流模塊,采用D8流向算法,用于確定水流的流向和流量分布。河道匯流模塊選擇了擴(kuò)散波法,用于模擬河道中水流的匯流和洪水傳播過程,地下徑流模塊采用基于概念的水桶模型。

    Noah-MP陸面模式主要模擬大氣下層、植被冠層和淺層土層之間的水量和能量的交互過程[28]。冠層截留、超滲地面徑流、土壤含水量和蒸散發(fā)均在陸面模式中計(jì)算。淺層土層被稱為次表面層,主要計(jì)算水量、能量的分布情況。當(dāng)雨強(qiáng)大于最大下滲率時(shí),地表出現(xiàn)超滲地面徑流。最大下滲率的計(jì)算公式為

    (1)

    式(1)中:Imax為最大下滲率,m/s;Pd為經(jīng)冠層截留后的雨強(qiáng),m/s;Dx為次表面層總?cè)彼?m/s;Ks為飽和水力傳導(dǎo)度,m/s;Kref= 2.0×10-6m/s,k為可調(diào)參數(shù);δi為與時(shí)段轉(zhuǎn)換有關(guān)的系數(shù)。次表面層的溫度和濕度分別由熱擴(kuò)散方程和Richard方程計(jì)算。Noah-MP模式的總蒸發(fā)分為地表直接蒸發(fā)、冠層截留量蒸發(fā)和植物根系散發(fā)3部分,由Penman公式計(jì)算。

    匯流網(wǎng)格主要計(jì)算壤中匯流、坡面匯流、地下徑流和河道匯流等,利用達(dá)西定律計(jì)算側(cè)向飽和壤中流,計(jì)算公式為

    (2)

    式(2)中:Qi為飽和壤中流的流量,m3/s;αx和αy分別為x和y方向的水面坡度;d為土層厚度,m;w為網(wǎng)格邊長,m;h為水面深度,m;n為可調(diào)參數(shù),默認(rèn)為1。在計(jì)算步長Δt內(nèi),根據(jù)水量平衡計(jì)算當(dāng)前網(wǎng)格的水深變化Δh。模型認(rèn)為每個(gè)網(wǎng)格的地表積水由3部分構(gòu)成:超滲地面產(chǎn)流、飽和地面流以及與其他網(wǎng)格的水量交換。當(dāng)且僅當(dāng)網(wǎng)格水深hs超過參數(shù)地表積水深hret時(shí),出現(xiàn)可自由流動(dòng)的地表徑流。地表徑流可由擴(kuò)散波方程計(jì)算,其中單寬流量與糙率的關(guān)系由曼寧公式計(jì)算,即

    (3)

    式(3)中:q為單寬流量,m2/s;nov為坡面的曼寧糙率;Sf為摩阻坡度。使用概念水桶模型描述流域地下徑流的調(diào)蓄作用,該水桶的水深并不代表真實(shí)含水層中的水深。水桶的流入量為次表面層的最下層向下的滲透量,流出量采用經(jīng)驗(yàn)蓄泄關(guān)系計(jì)算。河道形狀均采用等腰梯形,其底寬、邊坡系數(shù)、糙率為河道等級(jí)的函數(shù),在相關(guān)參數(shù)文件中有預(yù)定值,也可在隨后過程中率定。河道水演算采用一維擴(kuò)散波方程,即

    (4)

    式(4)中:Q為河道流量,m3/s;nch為河道曼寧糙率;A為過水?dāng)嗝婷娣e,m2;H為河道水面高程,m;R為水力半徑,m;sgn為符號(hào)函數(shù)。

    2.4 CaMa-Flood模型

    CaMa-Flood模型[29]是一個(gè)分布式的全球水動(dòng)力模型,CaMa-Flood模型的特點(diǎn)之一是其全球范圍的適用性,可以在不同的地理環(huán)境和流域條件下進(jìn)行應(yīng)用,模型利用高分辨率的數(shù)字高程模型和全球覆蓋的水系河網(wǎng)數(shù)據(jù),能夠精確地模擬河流的水動(dòng)力過程。CaMa-Flood模型的另一個(gè)特點(diǎn)是其在全球河流模擬中具有高計(jì)算效率,通過引入子網(wǎng)格地形參數(shù),將洪泛平原復(fù)雜的淹沒過程近似視為單元集水區(qū)的計(jì)算過程。CaMa-Flood模型的高計(jì)算效率對(duì)于需要進(jìn)行集合或長期高計(jì)算要求的實(shí)驗(yàn)特別有利,同時(shí)也方便了與其他水文方案的動(dòng)態(tài)耦合。

    現(xiàn)采用CaMa-Flood V4.10版本,首先將北江流域的河網(wǎng)數(shù)據(jù)離散成多個(gè)單元集水區(qū)(unit-catchments),基于1 min的全球河網(wǎng)數(shù)據(jù)和90 m的高程數(shù)據(jù)生成流域內(nèi)河道和洪泛區(qū)的高分辨率網(wǎng)格地形參數(shù),結(jié)合來自水文模型的3 km分辨率徑流量,利用局部慣性方程和自適應(yīng)時(shí)間步長方案,能夠快速模擬洪泛區(qū)的淹沒情況,其中采用自適應(yīng)時(shí)間步長方案以滿足CFL(Courant-Friedrichs-Lewy)條件,局部慣性方程通過離散化處理[30]并修正為

    (5)

    式(5)中:Qt+Δt為t次與t+Δt次之間的流量;Qt為前一時(shí)間步長的流量;A為水流截面積,m2;h為水流深度,m;g為重力加速度,m/s2;n為曼寧系數(shù),m1/3/s;S為上下游單元集水區(qū)之間的水面坡度。

    各時(shí)間步的計(jì)算流程為:① 根據(jù)各單元流域的蓄水量估算水深和淹沒面積;② 利用局部慣性方程計(jì)算河網(wǎng)圖中各單元流域到其下游單元流域的流量;③ 考慮上游單元流域的流量輸入、下游單元流域的流量輸出和地表徑流的輸入,通過質(zhì)量守恒方程更新每個(gè)單元流域的蓄水量。

    3 案例研究與應(yīng)用

    3.1 水文模型建立與優(yōu)化

    3.1.1 場(chǎng)次洪水選取

    WRF-Hydro模式需要預(yù)熱以獲得與實(shí)際較為相符的流域初始狀態(tài)[31],每一年的預(yù)熱期為汛期中有記錄的第一場(chǎng)洪水的前兩個(gè)月左右,模型預(yù)熱可以獲得相對(duì)準(zhǔn)確的初始條件,提高模擬精度。在選取洪水事件時(shí),考慮到洪峰較大的洪水可以更好地反映洪水事件的極端情況[32],因此挑選了洪峰流量超過10 000 m3/s的特大洪水作為研究事件。并將2001—2010年期間的7場(chǎng)洪水作為率定期洪水事件,將2011—2014年期間的5場(chǎng)洪水作為驗(yàn)證期洪水事件,如表1所示。率定期和驗(yàn)證期的驅(qū)動(dòng)數(shù)據(jù)采用GLDAS中除降水外的必要?dú)庀髷?shù)據(jù)和CMFD中的降水?dāng)?shù)據(jù),所提及的時(shí)間均為世界協(xié)調(diào)時(shí)間(universal time coordinated, UTC )。

    表1 模擬洪水編號(hào)和具體過程信息Table 1 Simulated flood numbers and specific process information

    3.1.2 率定參數(shù)選取

    WRF-Hydro中有許多物理參數(shù),但每個(gè)參數(shù)對(duì)洪水模擬的影響程度并不相同,因此參考Ryu等[33]、李致家等[34]對(duì)WRF-Hydro模式參數(shù)敏感性的分析,并結(jié)合預(yù)熱時(shí)北江流域的實(shí)際情況,選取對(duì)降雨-徑流模擬結(jié)果影響較大的下滲系數(shù)(REFKDT)、曼寧糙率系數(shù)(MANNFAC)、地表持水深度(RETDEPRTFAC)、地表糙率(OVROUGHRT)這4個(gè)參數(shù)進(jìn)行率定,具體率定方案如表2所示。

    表2 WRF-Hydro模型4個(gè)主要參數(shù)介紹以及率定的設(shè)置Table 2 The introduction of four main parameters of WRF-Hydro model and the setting of the rate

    3.1.3 率定過程分析

    考慮到WRF-hydro的計(jì)算成本較高,采用手動(dòng)逐步率定法[23]進(jìn)行參數(shù)率定,首先在合理的取值范圍內(nèi)[34]對(duì)某個(gè)參數(shù)按一定的步長進(jìn)行調(diào)參,其次根據(jù)所選擇的評(píng)價(jià)指標(biāo)來評(píng)價(jià)結(jié)果,選取出最優(yōu)的參數(shù)值,最終在此基礎(chǔ)上,對(duì)下一個(gè)參數(shù)進(jìn)行率定。采用皮爾遜相關(guān)系數(shù)(R)、納什效率系數(shù)(NSE)、峰現(xiàn)時(shí)間誤差(ΔT)、總流量相對(duì)誤差(RE)和洪峰流量相對(duì)誤差(PRE)作為評(píng)價(jià)指標(biāo)評(píng)估每次調(diào)參后的模擬結(jié)果,其函數(shù)表達(dá)式分別為

    (6)

    (7)

    (8)

    (9)

    (10)

    參考丁光旭[35]、高玉芳等[36]對(duì)WRF-Hydro模型的率定過程和目標(biāo)函數(shù)的研究,對(duì)于場(chǎng)次洪水的模擬,綜合考慮7場(chǎng)洪水過程平均的皮爾遜相關(guān)系數(shù)(fR)、納什效率系數(shù)(fNSE)、總流量相對(duì)誤差絕對(duì)值(f|RE|)、洪峰流量相對(duì)誤差絕對(duì)值(f|PRE|)、峰現(xiàn)時(shí)間誤差絕對(duì)值(f|ΔT|)這5個(gè)客觀函數(shù),對(duì)模擬結(jié)果進(jìn)行評(píng)估并選擇最優(yōu)值,函數(shù)表達(dá)式分別為

    (11)

    (12)

    (13)

    (14)

    (15)

    4個(gè)主要參數(shù)的率定過程如圖3所示,其中虛線標(biāo)記為選出的最優(yōu)值。下滲系數(shù)主要決定產(chǎn)流量的大小,因此考慮流量相對(duì)誤差絕對(duì)值和納什效率系數(shù)進(jìn)行優(yōu)選。如圖3(a)所示,在下滲系數(shù)為6時(shí),總流量相對(duì)誤差絕對(duì)值和洪峰誤差絕對(duì)值最小,納什效率系數(shù)最大,故取下滲系數(shù)為6。曼寧糙率主要影響水在河道的匯流速度,因此考慮峰現(xiàn)時(shí)間誤差絕對(duì)值和納什效率系數(shù)。從圖3(b)中可以觀察到,首先隨著曼寧糙率的增大,納什效率系數(shù)也增大,同時(shí)峰現(xiàn)時(shí)間誤差的絕對(duì)值逐漸減小,但在曼寧糙率達(dá)到1.3之后,納什效率開始降低,故取最優(yōu)值1.5。地表持水深度表征地表的蓄水能力,參考納什效率系數(shù)和相關(guān)系數(shù),由圖3(c)可以看出納什效率系數(shù)隨著地表持水深度的增大先增大后減小,相關(guān)系數(shù)先減小后增大,故取最優(yōu)值為3。地表糙率主要影響坡面的匯流速度,參考納什效率系數(shù)和峰現(xiàn)時(shí)間誤差絕對(duì)值,從圖3(d)中可以觀察到,隨著地表糙率的增大,納什效率系數(shù)先增大后減小,峰現(xiàn)時(shí)間誤差絕對(duì)值先減小后增大,故取最優(yōu)值2,最終得到4個(gè)主要參數(shù)的取值結(jié)果,如表3所示。

    圖3 4個(gè)主要參數(shù)的率定過程分析Fig.3 Analysis of the calibration process of four main parameters

    表3 4個(gè)主要參數(shù)的取值結(jié)果Table 3 Values of the four main parameters

    3.2 水文模型結(jié)果分析

    3.2.1 驗(yàn)證期結(jié)果分析

    2011—2014年5場(chǎng)洪水過程的徑流模擬結(jié)果如圖4所示,其中Q為站點(diǎn)流量、P為面平均降雨量,模擬徑流過程與降雨量的時(shí)間分布較為吻合。編號(hào)20110512洪水和編號(hào)20130513洪水為典型的單洪峰過程,降雨-徑流過程相對(duì)簡單,洪水的漲水過程模擬較優(yōu),對(duì)洪峰流量也把握較準(zhǔn),但落水過程相對(duì)較慢,可能是因?yàn)閷?duì)應(yīng)時(shí)段內(nèi)降雨量仍然較大,而模型對(duì)降雨條件比較敏感,因此落水過程較緩慢。編號(hào)20120616洪水和編號(hào)20130809洪水這兩次過程的模擬與實(shí)測(cè)過程線基本吻合,但在漲水期模擬流量偏小,從圖4(b)和圖4(e)中可以觀察到,這兩場(chǎng)洪水在降雨期間都是出現(xiàn)了多個(gè)時(shí)段的降雨峰值,對(duì)應(yīng)的降雨徑流過程相對(duì)復(fù)雜,導(dǎo)致洪水的漲水過程較難模擬。編號(hào)20140513該場(chǎng)洪水為多峰洪水,第一次模擬洪峰出現(xiàn)的較晚且洪峰要大,但第二次洪峰的漲退水過程模擬良好,可能是由于模型的初始狀態(tài)與實(shí)際狀態(tài)存在一定誤差,從而延遲了對(duì)第一次洪峰的到達(dá)。整體來說,該模型可以較好地捕捉到洪峰,且峰現(xiàn)時(shí)間誤差時(shí)間較小,但整體的漲退水過程還有一定的偏差。

    圖4 驗(yàn)證期5場(chǎng)洪水的降雨徑流模擬結(jié)果Fig.4 Simulation results of rainfall runoff from 5 floods in verification period

    驗(yàn)證期5場(chǎng)洪水過程的逐小時(shí)評(píng)價(jià)結(jié)果如表4所示,模擬徑流和實(shí)測(cè)徑流過程的相關(guān)系數(shù)均在0.8以上,納什效率系數(shù)介于0.52~0.71,峰現(xiàn)時(shí)間誤差絕對(duì)值都在10 h以內(nèi),洪峰流量誤差穩(wěn)定在12.6%以內(nèi),對(duì)于洪水期間總流量誤差,有3場(chǎng)均在7.93%以下,兩場(chǎng)分別為38.54%、-21.82%。

    表4 驗(yàn)證期5場(chǎng)洪水模擬的指標(biāo)評(píng)價(jià)Table 4 Evaluation of 5 flood simulation in verification phase

    3.2.2 預(yù)報(bào)研究分析

    為了進(jìn)一步研究優(yōu)化后的WRF-Hydro模型的預(yù)報(bào)性能,將經(jīng)過率定的模型參數(shù)應(yīng)用到后續(xù)研究中。由于CMFD數(shù)據(jù)集僅覆蓋到2018年,因此將CMFD降水?dāng)?shù)據(jù)替換為GLDAS降水?dāng)?shù)據(jù),并以2022年6月北江流域“龍舟水”期間的3次編號(hào)洪水作為案例,旨在探究更換降水產(chǎn)品后模型的預(yù)報(bào)效果是否仍然良好,由于缺乏完整的實(shí)測(cè)徑流過程,所以考慮用實(shí)測(cè)洪峰流量來對(duì)比。對(duì)于北江流域“龍舟水”期間的徑流模擬結(jié)果如圖5所示,可以看到出現(xiàn)了3次洪峰[37-38],其中提到了北江流域“龍舟水”期間出現(xiàn)了3次洪水,分別在6月15日、6月22日和7月7日,預(yù)報(bào)結(jié)果與實(shí)際情況相符,并且文獻(xiàn)[37]提供了石角水文站實(shí)測(cè)3次編號(hào)洪水的天然洪峰流量分別為15 400、20 700、15 100 m3/s。

    圖5 2022年北江流域“龍舟水”期間降雨徑流關(guān)系Fig.5 Relationship of rainfall runoff during Dragon Boat Water in Beijiang River Basin in 2022

    該次模擬并沒有考慮人為影響和水庫的調(diào)度,所以屬于天然洪水過程。模擬洪峰流量和實(shí)測(cè)洪峰流量的對(duì)比如圖6所示,這3場(chǎng)洪峰的誤差分別為-31.9%、15.2%、-19.5%,平均誤差絕對(duì)值為22.2%,其中1號(hào)洪峰誤差較大,從圖6中可以看到在1號(hào)洪水期間,降水量較后兩場(chǎng)偏小,這可能和GLDAS降水?dāng)?shù)據(jù)與實(shí)際降水之間的差異性有關(guān),導(dǎo)致對(duì)1號(hào)洪峰模擬的誤差較大。綜合來說,該水文模型在更換較粗分辨率的降水產(chǎn)品之后,預(yù)報(bào)性能也達(dá)到了合理的效果,為北江流域洪水預(yù)報(bào)提供有價(jià)值的參考。

    圖6 2022年北江流域“龍舟水”期間3次洪峰流量對(duì)比Fig.6 Comparison of three flood peaks discharge during “Dragon Boat Water” in Beijian Basin in 2022

    3.3 水文水動(dòng)力模型耦合

    3.3.1 耦合模型建立

    CaMa-Flood模型的邊界條件通?;谒哪P偷膹搅鬟^程[39]。在WRF-Hydro模式中,Noah-MP模式考慮了地表徑流和地下徑流的形成過程[40]。因此,在耦合的過程中實(shí)際上是通過采用WRF-Hydro模式中的Noah-MP陸面模式來與CaMa-Flood模型進(jìn)行耦合。這種單向耦合方法是基于Noah-MP模塊輸出的3 km水平分辨率和3 h時(shí)間分辨率的地表和地下總徑流過程,作為CaMa-Flood模型的邊界條件。借助這樣的耦合方式,CaMa-Flood模型能夠以0.1°分辨率模擬整個(gè)流域內(nèi)的洪水淹沒深度和范圍。

    3.3.2 洪水淹沒分析

    2022年6月13日—2022年6月28日期間北江流域的洪水淹沒深度變化如圖7所示,洪水淹沒范圍變化如圖8所示,該淹沒制圖能大致反映出洪泛區(qū)的洪水淹沒變化過程。在這段時(shí)間內(nèi),洪水淹沒過程的演變與降雨-徑流的變化情況密切相關(guān),如圖5所示,在6月15日一號(hào)洪峰到達(dá)后,仍持續(xù)數(shù)天強(qiáng)降雨,在強(qiáng)降雨的作用下,洪水流量也逐漸增大,從而導(dǎo)致洪水淹沒區(qū)域的擴(kuò)大。如圖7(h)、圖8(h)所示,在6月20日洪水淹沒達(dá)到了高峰,這說明強(qiáng)降雨的影響,導(dǎo)致洪水淹沒深度和范圍逐漸增大,然而隨著時(shí)間的推移,6月22日之后降雨量逐漸減小一直到停止,與此同時(shí)洪水流量也開始逐漸減小。需要注意的是,洪水淹沒情況的反映往往存在一定的滯后性且由于地勢(shì)等因素的影響,洪水淹沒區(qū)域并不會(huì)立即減少,而是需要一段時(shí)間才能逐漸恢復(fù)正常。由圖8可以觀察到,在6月22日之后,洪水淹沒范圍開始逐漸減少,這表明洪水強(qiáng)度與洪水淹沒之間存在合理的相關(guān)性。綜上所述,洪水淹沒過程與降雨強(qiáng)度和洪峰流量的變化密切相關(guān),這種耦合模型的合理性得到了驗(yàn)證,提供了預(yù)測(cè)洪水后洪泛區(qū)淹沒情況的參考依據(jù)。

    圖7 2022年北江流域“龍舟水”期間洪水淹沒深度逐日變化情況Fig.7 Daily variation of flood inundation depth during “Dragon Boat Water” in Beijiang Basin in 2022

    圖8 2022年北江流域“龍舟水”期間洪水淹沒范圍逐日變化情況Fig.8 Daily variation of flood inundation area during Dragon Boat Water in Beijiang River Basin in 2022

    4 結(jié)論

    以珠江流域的第二大水系北江為研究對(duì)象,首先利用WRF-Hydro模式構(gòu)建北江流域水文模型,經(jīng)過7場(chǎng)洪水事件的率定,得到優(yōu)化后的模型參數(shù),對(duì)5場(chǎng)洪水事件進(jìn)行驗(yàn)證,并探究其水文預(yù)報(bào)性能。接下來利用優(yōu)化后的WRF-Hydro模型耦合CaMa-Flood模型,快速計(jì)算出洪水演進(jìn)過程中洪泛區(qū)的淹沒深度和淹沒范圍,得到如下主要結(jié)論。

    (1) 水文模型的結(jié)果表明徑流模擬的效果較好,驗(yàn)證期的5場(chǎng)洪水過程納什效率系數(shù)在0.52~0.71,相關(guān)系數(shù)均在0.8以上,洪峰流量誤差均在10%以內(nèi),峰現(xiàn)時(shí)間誤差也都在10 h以內(nèi)。對(duì)于預(yù)報(bào)期洪水,3次洪峰的平均誤差為22.2%。由于WRF-Hydro模式對(duì)于降水?dāng)?shù)據(jù)較為敏感,未來可以采用更高的質(zhì)量的降水?dāng)?shù)據(jù)進(jìn)行預(yù)報(bào),以更好地發(fā)揮模式性能。

    (2) 水文水動(dòng)力耦合模型中,前者為后者提供數(shù)據(jù),前者的參數(shù)優(yōu)化和模擬的準(zhǔn)確性尤為重要,一旦確保水動(dòng)力模型的驅(qū)動(dòng)數(shù)據(jù)相對(duì)準(zhǔn)確,可以進(jìn)行迅速大范圍的洪泛區(qū)淹沒模擬,對(duì)受災(zāi)區(qū)提供有價(jià)值的信息,如洪水淹沒深度、淹沒范圍和淹沒持續(xù)時(shí)間。

    (3) 洪水淹沒通常被視為洪水到來之后的突發(fā)事件,在北江流域這種地形復(fù)雜、降水分布不均勻的流域,構(gòu)建的水文水動(dòng)力耦合模型對(duì)洪泛區(qū)的淹沒過程也做出了合理的模擬。但由于缺乏實(shí)測(cè)淹沒數(shù)據(jù),未能詳細(xì)驗(yàn)證洪泛區(qū)的淹沒情況,未來的研究將考慮與其他水動(dòng)力模型的計(jì)算結(jié)果、衛(wèi)星影像等進(jìn)行對(duì)比。

    (4) 當(dāng)前水文水動(dòng)力耦合預(yù)報(bào)領(lǐng)域的研究相對(duì)不足,在前人對(duì)WRF-Hydro模式和CaMa-Flood模型的研究基礎(chǔ)上,首次將WRF-Hydro模式與CaMa-Flood模型進(jìn)行耦合,以經(jīng)過優(yōu)化的WRF-Hydro水文模型中的Noah-MP模式的輸出徑流過程作為CaMa-Flood模型的邊界條件,構(gòu)建了降雨徑流-洪水淹沒綜合預(yù)報(bào)模型。在未來的研究中,將進(jìn)一步探索多種耦合方案,以提升預(yù)報(bào)的準(zhǔn)確性和精度。

    (5) 值得進(jìn)一步指出的是,WRF-Hydro與WRF“在線”耦合模式下,可以再結(jié)合CaMa-Flood模型,實(shí)現(xiàn)氣象模式與水文水動(dòng)力模型的耦合預(yù)報(bào),延長洪水預(yù)報(bào)的預(yù)見期。

    構(gòu)建的WRF-Hydro水文模型在洪峰捕捉方面較優(yōu)秀,雖然漲退水過程有略有偏差,但總體能夠大致反映實(shí)際洪水過程,對(duì)洪水的量級(jí)的模擬較為準(zhǔn)確。水文水動(dòng)力模型耦合之后,隨著隆雨-徑流過程的變化,能夠快速、合理地模擬洪泛區(qū)的淹沒情況。然而,在以精細(xì)分辨率模擬大尺度洪水淹沒方面一直存在困難,CaMa-Flood模型本身是一種簡化的水動(dòng)力模型,基于數(shù)字高程模型假設(shè)河道和洪泛平原的斷面形狀,這種簡化意味著包括永久性湖泊和水庫在內(nèi)的一些地勢(shì)低洼區(qū)域會(huì)被模型視為河流或洪泛平原。此外,當(dāng)前版本的CaMa-Flood 模型僅能模擬自然洪水,未考慮防洪設(shè)施如堤壩、水庫等的影響,盡管WRF-Hydro模式已納入這些因素,但鑒于在缺乏實(shí)測(cè)數(shù)據(jù)和人為影響難以準(zhǔn)確把握的情況下,暫未引入這些影響因素,在未來的研究中應(yīng)考慮將這些因素納入,以更準(zhǔn)確地評(píng)估洪水的影響??傮w而言,在集水面積4萬多平方公里、地形復(fù)雜、降水時(shí)空分布不均勻的北江流域,模擬結(jié)果較為合理,為北江流域洪水預(yù)報(bào)、水資源管理、風(fēng)險(xiǎn)評(píng)估等研究提供一定參考價(jià)值。

    猜你喜歡
    北江水文徑流
    2022年《中國水文年報(bào)》發(fā)布
    水文
    水文水資源管理
    北江,向前
    奔騰北江
    泥娃娃
    暢談(2018年11期)2018-08-26 02:10:28
    水文
    Topmodel在布哈河流域徑流模擬中的應(yīng)用
    基于MODIS北江流域土地覆蓋變化研究及CA-Markov預(yù)測(cè)
    河北遙感(2015年2期)2015-07-18 11:11:14
    探秘“大徑流”

    午夜福利在线在线| 亚洲 欧美 日韩 在线 免费| 国产精品免费视频内射| 少妇 在线观看| 欧美精品亚洲一区二区| 午夜激情av网站| 久久精品91蜜桃| 日韩中文字幕欧美一区二区| 精品国产亚洲在线| 日本五十路高清| 一a级毛片在线观看| 波多野结衣高清作品| 亚洲狠狠婷婷综合久久图片| 中文字幕另类日韩欧美亚洲嫩草| 人成视频在线观看免费观看| 国产精品久久视频播放| 好男人在线观看高清免费视频 | 欧美乱码精品一区二区三区| 国产亚洲精品综合一区在线观看 | 一级a爱视频在线免费观看| 精品国产超薄肉色丝袜足j| 亚洲 国产 在线| 成人三级黄色视频| 午夜免费观看网址| 亚洲精品国产精品久久久不卡| 亚洲无线在线观看| 少妇裸体淫交视频免费看高清 | 99国产综合亚洲精品| 国产三级黄色录像| 嫩草影院精品99| 免费看美女性在线毛片视频| av视频在线观看入口| 国产精品 国内视频| 大型av网站在线播放| 久久人妻福利社区极品人妻图片| 欧美zozozo另类| 免费搜索国产男女视频| 国产成人影院久久av| a在线观看视频网站| 禁无遮挡网站| 三级毛片av免费| 日韩一卡2卡3卡4卡2021年| 在线av久久热| 国产精品久久电影中文字幕| 亚洲av电影不卡..在线观看| 黄色a级毛片大全视频| 可以在线观看毛片的网站| 精品久久久久久久人妻蜜臀av| 精品久久久久久久人妻蜜臀av| 免费一级毛片在线播放高清视频| 久久久久久久久中文| 国产精品国产高清国产av| 欧美激情久久久久久爽电影| 手机成人av网站| www日本在线高清视频| 欧美不卡视频在线免费观看 | 亚洲一区二区三区色噜噜| 午夜影院日韩av| 久久国产精品男人的天堂亚洲| 久久精品国产清高在天天线| 久久久水蜜桃国产精品网| 90打野战视频偷拍视频| 精品免费久久久久久久清纯| 亚洲国产精品成人综合色| 在线观看午夜福利视频| tocl精华| 99久久国产精品久久久| 女人被狂操c到高潮| www.熟女人妻精品国产| 精品一区二区三区av网在线观看| 国产1区2区3区精品| 欧美性长视频在线观看| 一区二区三区激情视频| 男女午夜视频在线观看| av在线播放免费不卡| 97人妻精品一区二区三区麻豆 | 麻豆av在线久日| 欧美色视频一区免费| 韩国精品一区二区三区| 人人妻人人澡欧美一区二区| 亚洲国产精品999在线| 国产精品久久久av美女十八| 亚洲成人国产一区在线观看| videosex国产| 大型av网站在线播放| 久9热在线精品视频| 国产亚洲欧美98| 日韩欧美国产在线观看| √禁漫天堂资源中文www| 人人澡人人妻人| av电影中文网址| 精品一区二区三区视频在线观看免费| 久久精品人妻少妇| av欧美777| 亚洲成av片中文字幕在线观看| 哪里可以看免费的av片| 欧美日韩中文字幕国产精品一区二区三区| 日日干狠狠操夜夜爽| 日韩高清综合在线| 美女 人体艺术 gogo| 亚洲天堂国产精品一区在线| 精品国产国语对白av| 黄色成人免费大全| 99久久久亚洲精品蜜臀av| 亚洲国产精品sss在线观看| 精品国产美女av久久久久小说| 亚洲av电影在线进入| 亚洲国产欧美一区二区综合| 麻豆成人午夜福利视频| 国产区一区二久久| 久久热在线av| 欧美性长视频在线观看| 久久久久久九九精品二区国产 | 香蕉av资源在线| 91成年电影在线观看| 91国产中文字幕| 亚洲男人的天堂狠狠| 色综合婷婷激情| 精品一区二区三区视频在线观看免费| 一级a爱片免费观看的视频| 欧美激情极品国产一区二区三区| 亚洲精品久久国产高清桃花| 日韩一卡2卡3卡4卡2021年| 看免费av毛片| 欧美激情久久久久久爽电影| 老汉色∧v一级毛片| 国产成人系列免费观看| 啦啦啦 在线观看视频| 免费人成视频x8x8入口观看| 久99久视频精品免费| 亚洲精品国产一区二区精华液| 男男h啪啪无遮挡| 国产视频一区二区在线看| 国产在线观看jvid| 色综合站精品国产| 国产色视频综合| 成人三级做爰电影| 9191精品国产免费久久| 亚洲中文字幕日韩| 精华霜和精华液先用哪个| 国产一区二区激情短视频| 日韩大码丰满熟妇| 亚洲第一欧美日韩一区二区三区| 日韩大尺度精品在线看网址| 757午夜福利合集在线观看| 国产精品久久久久久亚洲av鲁大| 宅男免费午夜| 欧美日韩亚洲综合一区二区三区_| 国产色视频综合| 精品久久久久久久久久久久久 | 草草在线视频免费看| а√天堂www在线а√下载| 麻豆成人午夜福利视频| 精品第一国产精品| 亚洲国产精品合色在线| 人人妻人人澡欧美一区二区| 久久久久久大精品| 成年人黄色毛片网站| www.精华液| 久久精品国产综合久久久| 一进一出抽搐gif免费好疼| 婷婷六月久久综合丁香| 国产高清有码在线观看视频 | 亚洲国产毛片av蜜桃av| 一边摸一边做爽爽视频免费| 亚洲美女黄片视频| 国产亚洲欧美在线一区二区| 最近在线观看免费完整版| 一本精品99久久精品77| 日本熟妇午夜| 久久国产精品男人的天堂亚洲| 黄频高清免费视频| 欧美丝袜亚洲另类 | 午夜激情福利司机影院| 亚洲精品中文字幕一二三四区| 国产成人一区二区三区免费视频网站| 级片在线观看| 欧美精品啪啪一区二区三区| 亚洲国产精品999在线| 露出奶头的视频| 午夜精品久久久久久毛片777| 久久久久久久久免费视频了| 日本成人三级电影网站| 亚洲第一av免费看| 亚洲七黄色美女视频| 性欧美人与动物交配| 亚洲精品av麻豆狂野| 最新在线观看一区二区三区| √禁漫天堂资源中文www| 色老头精品视频在线观看| av天堂在线播放| www.999成人在线观看| 99久久无色码亚洲精品果冻| 亚洲第一欧美日韩一区二区三区| 亚洲av成人一区二区三| 国产精品久久久久久亚洲av鲁大| 亚洲熟妇熟女久久| 在线观看午夜福利视频| 国产精品永久免费网站| 国产精品一区二区精品视频观看| 精华霜和精华液先用哪个| 国产伦人伦偷精品视频| 女警被强在线播放| 精品国内亚洲2022精品成人| 亚洲 欧美 日韩 在线 免费| 国产精品久久久久久精品电影 | 窝窝影院91人妻| 熟女电影av网| 99热只有精品国产| 这个男人来自地球电影免费观看| 午夜成年电影在线免费观看| 激情在线观看视频在线高清| 午夜a级毛片| 亚洲黑人精品在线| 久久精品国产99精品国产亚洲性色| 久久精品亚洲精品国产色婷小说| 日本五十路高清| 欧美日韩乱码在线| 亚洲欧美精品综合久久99| 韩国av一区二区三区四区| 最好的美女福利视频网| 哪里可以看免费的av片| 亚洲色图 男人天堂 中文字幕| 黄色a级毛片大全视频| 国产成人欧美| 亚洲av成人av| 免费在线观看黄色视频的| 免费在线观看日本一区| 18禁观看日本| 亚洲人成77777在线视频| 无人区码免费观看不卡| 中文字幕高清在线视频| 男人舔奶头视频| 极品教师在线免费播放| av片东京热男人的天堂| 成人永久免费在线观看视频| 日韩成人在线观看一区二区三区| 成人午夜高清在线视频 | 村上凉子中文字幕在线| 黄频高清免费视频| 国产精品久久电影中文字幕| 国产亚洲精品久久久久久毛片| 国产91精品成人一区二区三区| 久久久久九九精品影院| 国产aⅴ精品一区二区三区波| 在线观看www视频免费| 国产av不卡久久| 国产成人一区二区三区免费视频网站| 国产亚洲精品av在线| www日本黄色视频网| 色精品久久人妻99蜜桃| 男女那种视频在线观看| 国产成人欧美| 国产精品爽爽va在线观看网站 | 亚洲午夜精品一区,二区,三区| 在线永久观看黄色视频| www.www免费av| 国产成人av激情在线播放| 波多野结衣高清无吗| 久久久久久久精品吃奶| 国产日本99.免费观看| 在线十欧美十亚洲十日本专区| 国产蜜桃级精品一区二区三区| 嫁个100分男人电影在线观看| 国产aⅴ精品一区二区三区波| 亚洲成人久久性| 99国产综合亚洲精品| 国产成人啪精品午夜网站| 变态另类成人亚洲欧美熟女| 88av欧美| 国产麻豆成人av免费视频| 中文资源天堂在线| 国产精品久久久av美女十八| 欧美黑人巨大hd| 亚洲在线自拍视频| 亚洲专区中文字幕在线| 90打野战视频偷拍视频| 人妻久久中文字幕网| 欧美成人午夜精品| 亚洲精品国产精品久久久不卡| 日本一区二区免费在线视频| 韩国av一区二区三区四区| 精华霜和精华液先用哪个| 久久精品成人免费网站| 日日爽夜夜爽网站| 日韩欧美国产一区二区入口| 亚洲精品中文字幕一二三四区| 成人18禁在线播放| 丝袜在线中文字幕| 搡老妇女老女人老熟妇| 久久午夜综合久久蜜桃| 亚洲成人久久性| 国产黄色小视频在线观看| 精品久久久久久久久久久久久 | 亚洲九九香蕉| 国产久久久一区二区三区| 每晚都被弄得嗷嗷叫到高潮| 99精品欧美一区二区三区四区| 大香蕉久久成人网| www.熟女人妻精品国产| 麻豆久久精品国产亚洲av| 黑丝袜美女国产一区| 国产亚洲欧美精品永久| 亚洲 欧美 日韩 在线 免费| 亚洲av日韩精品久久久久久密| xxxwww97欧美| 侵犯人妻中文字幕一二三四区| 亚洲无线在线观看| 亚洲国产中文字幕在线视频| 日韩国内少妇激情av| 精品少妇一区二区三区视频日本电影| 亚洲精品色激情综合| 一级片免费观看大全| 黄色丝袜av网址大全| 大型黄色视频在线免费观看| 极品教师在线免费播放| 黄色丝袜av网址大全| 国产精品久久久久久亚洲av鲁大| 99国产精品一区二区蜜桃av| a级毛片在线看网站| 欧美zozozo另类| 成人18禁高潮啪啪吃奶动态图| 亚洲一区二区三区色噜噜| 一个人免费在线观看的高清视频| 欧美成人免费av一区二区三区| 久久精品91蜜桃| 国产精品av久久久久免费| 老司机福利观看| 国产成人欧美| 性色av乱码一区二区三区2| 日韩欧美在线二视频| 少妇粗大呻吟视频| 亚洲aⅴ乱码一区二区在线播放 | 精品国产乱子伦一区二区三区| 老司机靠b影院| 国产精品永久免费网站| 超碰成人久久| 国产精品二区激情视频| 国产成人欧美在线观看| 国产亚洲av高清不卡| 免费高清在线观看日韩| 欧美国产日韩亚洲一区| 91字幕亚洲| 久久精品国产亚洲av高清一级| 黄色丝袜av网址大全| 亚洲真实伦在线观看| 欧美乱妇无乱码| 无限看片的www在线观看| 亚洲精品色激情综合| 夜夜夜夜夜久久久久| 麻豆成人av在线观看| 叶爱在线成人免费视频播放| 久久久久精品国产欧美久久久| 亚洲 国产 在线| 老鸭窝网址在线观看| 最近最新中文字幕大全电影3 | 淫秽高清视频在线观看| 久久久久精品国产欧美久久久| 男人舔奶头视频| 嫁个100分男人电影在线观看| 不卡av一区二区三区| 国产av不卡久久| 久久久国产成人精品二区| 一二三四在线观看免费中文在| 琪琪午夜伦伦电影理论片6080| 久久人妻av系列| 精品久久久久久,| 男人操女人黄网站| 欧美精品啪啪一区二区三区| 长腿黑丝高跟| www.精华液| 成人国产一区最新在线观看| 老鸭窝网址在线观看| 91在线观看av| 亚洲成人国产一区在线观看| 首页视频小说图片口味搜索| 国产亚洲精品久久久久久毛片| 日韩欧美免费精品| 夜夜躁狠狠躁天天躁| 美女午夜性视频免费| bbb黄色大片| 亚洲av五月六月丁香网| 久久香蕉国产精品| 亚洲男人的天堂狠狠| 亚洲专区中文字幕在线| 老鸭窝网址在线观看| 波多野结衣巨乳人妻| 午夜福利一区二区在线看| 国产日本99.免费观看| 动漫黄色视频在线观看| 国产高清有码在线观看视频 | 国产日本99.免费观看| 在线观看www视频免费| 俄罗斯特黄特色一大片| 久久国产精品人妻蜜桃| 久久国产精品影院| 99热只有精品国产| 非洲黑人性xxxx精品又粗又长| 黄片播放在线免费| 国产精品久久久av美女十八| 亚洲色图 男人天堂 中文字幕| 老熟妇乱子伦视频在线观看| videosex国产| 久久伊人香网站| 人人妻,人人澡人人爽秒播| 久久久久九九精品影院| 在线十欧美十亚洲十日本专区| 国产黄a三级三级三级人| 欧美另类亚洲清纯唯美| 此物有八面人人有两片| 久久性视频一级片| 亚洲精品av麻豆狂野| 18禁国产床啪视频网站| xxxwww97欧美| 国语自产精品视频在线第100页| 欧美性猛交黑人性爽| 高清在线国产一区| 精品久久久久久,| 久99久视频精品免费| 国产视频一区二区在线看| 国产精品av久久久久免费| 欧美一级a爱片免费观看看 | 国产精品久久久久久精品电影 | 亚洲国产精品sss在线观看| 国产色视频综合| 国产精品国产高清国产av| 久久久久久久午夜电影| 色尼玛亚洲综合影院| 日韩欧美国产一区二区入口| 成人手机av| 国产精品99久久99久久久不卡| 午夜福利免费观看在线| 日本在线视频免费播放| 国产精品久久久久久亚洲av鲁大| 一区二区三区精品91| 国产精品国产高清国产av| 免费高清在线观看日韩| 国产成人av激情在线播放| 草草在线视频免费看| 久久久国产成人精品二区| 精品电影一区二区在线| 午夜福利视频1000在线观看| 欧美日韩瑟瑟在线播放| 黑丝袜美女国产一区| 国产1区2区3区精品| 国产精品二区激情视频| 一卡2卡三卡四卡精品乱码亚洲| 久久人人精品亚洲av| 禁无遮挡网站| 精品国产国语对白av| 老汉色av国产亚洲站长工具| 日韩欧美在线二视频| 亚洲国产毛片av蜜桃av| 热99re8久久精品国产| 在线观看免费视频日本深夜| 亚洲欧美一区二区三区黑人| 国产成人精品久久二区二区免费| 99热只有精品国产| 久久久久精品国产欧美久久久| 丁香欧美五月| 日本撒尿小便嘘嘘汇集6| 国产高清视频在线播放一区| 亚洲五月婷婷丁香| 精品福利观看| 亚洲精品中文字幕一二三四区| netflix在线观看网站| 19禁男女啪啪无遮挡网站| 制服人妻中文乱码| 久久精品夜夜夜夜夜久久蜜豆 | 麻豆av在线久日| 亚洲va日本ⅴa欧美va伊人久久| 国产精品自产拍在线观看55亚洲| 可以在线观看毛片的网站| 宅男免费午夜| 成人18禁在线播放| 99re在线观看精品视频| 美女 人体艺术 gogo| 一本大道久久a久久精品| 一级a爱片免费观看的视频| 香蕉丝袜av| www国产在线视频色| 免费观看精品视频网站| 变态另类成人亚洲欧美熟女| www日本黄色视频网| 免费高清在线观看日韩| 日韩有码中文字幕| 最近最新中文字幕大全电影3 | 一区二区日韩欧美中文字幕| 99久久精品国产亚洲精品| 亚洲 国产 在线| 可以在线观看的亚洲视频| 午夜视频精品福利| 后天国语完整版免费观看| 身体一侧抽搐| 两性夫妻黄色片| 久久精品91蜜桃| 国产成人影院久久av| 亚洲精品一卡2卡三卡4卡5卡| a级毛片在线看网站| 国产精品亚洲一级av第二区| 日韩高清综合在线| 麻豆成人av在线观看| 午夜福利在线观看吧| 国产精品98久久久久久宅男小说| 一个人观看的视频www高清免费观看 | 丰满人妻熟妇乱又伦精品不卡| 亚洲片人在线观看| 欧美一级毛片孕妇| 97人妻精品一区二区三区麻豆 | 757午夜福利合集在线观看| 亚洲七黄色美女视频| 色婷婷久久久亚洲欧美| 满18在线观看网站| 亚洲一区中文字幕在线| 欧美zozozo另类| 亚洲国产精品合色在线| а√天堂www在线а√下载| 两性午夜刺激爽爽歪歪视频在线观看 | 国产精品久久久av美女十八| 欧美丝袜亚洲另类 | av在线天堂中文字幕| 国产亚洲精品综合一区在线观看 | 精品免费久久久久久久清纯| 亚洲一区二区三区不卡视频| 中文字幕另类日韩欧美亚洲嫩草| 不卡av一区二区三区| 99国产精品一区二区蜜桃av| 婷婷亚洲欧美| 啪啪无遮挡十八禁网站| 欧美日韩精品网址| 满18在线观看网站| 手机成人av网站| 一边摸一边抽搐一进一小说| 变态另类成人亚洲欧美熟女| 在线观看午夜福利视频| 欧美激情久久久久久爽电影| 亚洲精品国产精品久久久不卡| 精品第一国产精品| 岛国视频午夜一区免费看| 狂野欧美激情性xxxx| 亚洲精品国产区一区二| 久久久国产成人免费| 18禁观看日本| 久久久精品国产亚洲av高清涩受| 午夜免费成人在线视频| 在线av久久热| 男女午夜视频在线观看| 欧美中文日本在线观看视频| 欧美在线黄色| www.熟女人妻精品国产| 1024香蕉在线观看| 欧美一级a爱片免费观看看 | 黑丝袜美女国产一区| 国产亚洲精品一区二区www| 亚洲人成电影免费在线| 变态另类成人亚洲欧美熟女| 亚洲成av人片免费观看| 一区二区三区高清视频在线| 香蕉国产在线看| 亚洲五月色婷婷综合| 国产成人影院久久av| 中文亚洲av片在线观看爽| 亚洲色图 男人天堂 中文字幕| 亚洲国产欧美日韩在线播放| 在线观看免费午夜福利视频| 婷婷精品国产亚洲av在线| 夜夜夜夜夜久久久久| 99re在线观看精品视频| 日韩欧美在线二视频| 91麻豆精品激情在线观看国产| xxx96com| 美女高潮到喷水免费观看| 国产伦一二天堂av在线观看| 国产高清有码在线观看视频 | 一级毛片精品| 亚洲av电影在线进入| 亚洲自偷自拍图片 自拍| 啪啪无遮挡十八禁网站| 夜夜看夜夜爽夜夜摸| 在线观看日韩欧美| 自线自在国产av| 亚洲一区高清亚洲精品| 日韩精品青青久久久久久| 国产精品亚洲av一区麻豆| 在线观看一区二区三区| 国产高清videossex| 宅男免费午夜| 成人一区二区视频在线观看| av在线天堂中文字幕| 免费在线观看视频国产中文字幕亚洲| 满18在线观看网站| 又紧又爽又黄一区二区| 久久香蕉激情| 午夜影院日韩av| 欧美日韩亚洲国产一区二区在线观看| 国产精品久久久人人做人人爽| 午夜免费鲁丝| 久久久国产精品麻豆| www日本在线高清视频| 视频在线观看一区二区三区| 国产一区二区三区视频了| 99热6这里只有精品| 精品国产一区二区三区四区第35| 欧美一级毛片孕妇| 久久天堂一区二区三区四区| 啦啦啦观看免费观看视频高清| 久久久国产精品麻豆| 女性被躁到高潮视频| 大香蕉久久成人网| 久久欧美精品欧美久久欧美| 精品午夜福利视频在线观看一区| 精品国产美女av久久久久小说| 免费在线观看日本一区| 久久人妻福利社区极品人妻图片| 亚洲专区国产一区二区|