亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        格林-安普特降雨徑流模型改進(jìn)及初步應(yīng)用

        2020-09-24 10:39:36李致家霍文博
        關(guān)鍵詞:志丹普特產(chǎn)流

        李致家,霍文博,張 珂

        (河海大學(xué)水文水資源學(xué)院,江蘇 南京 210098)

        水文模型是對(duì)自然界中水文過(guò)程的模擬和概化。自19世紀(jì)50年代以來(lái),水文模型在國(guó)內(nèi)外發(fā)展迅速并廣泛應(yīng)用[1-6]。大部分水文模型在濕潤(rùn)地區(qū)應(yīng)用效果較好[7],而在干旱和半干旱地區(qū),由于降雨時(shí)空分布極不均勻、實(shí)測(cè)資料不足、產(chǎn)流模式比較復(fù)雜,導(dǎo)致水文模型模擬效果不佳[8-11]。在干旱和半干旱地區(qū),由于土層相對(duì)較厚,土壤含水量較低,在一場(chǎng)降雨中土壤很難蓄滿,因此徑流主要來(lái)源于當(dāng)降雨強(qiáng)度大于土壤下滲能力時(shí)形成的超滲地面徑流。李致家等[12]在新安江模型中增加了超滲產(chǎn)流模型,使其可以應(yīng)用于干旱和半干旱地區(qū)。超滲產(chǎn)流模型包括霍頓(Horton)模型、格林-安普特(Green-Ampt)模型、菲利普(Philip)模型和霍爾坦(Holtan)模型等,其中格林-安普特模型是由格林和安普特于1911年提出的基于物理基礎(chǔ)的下滲模型[13],國(guó)內(nèi)外對(duì)于格林-安普特模型的應(yīng)用和改進(jìn)研究較多。包為民[14]提出了一條具有流域分布特征的格林-安普特下滲曲線,并在黃河中游13個(gè)流域應(yīng)用,發(fā)現(xiàn)這條分布曲線結(jié)構(gòu)合理,能使格林-安普特降雨徑流模型模擬精度顯著提高;Stewart[15]通過(guò)研究土壤收縮和膨脹對(duì)飽和導(dǎo)水率的影響,提出了動(dòng)態(tài)格林-安普特下滲模型,與經(jīng)典格林-安普特模型相比,動(dòng)態(tài)格林-安普特模型在多尺度上計(jì)算入滲和徑流;Gowdish等[16]開(kāi)發(fā)了三維格林-安普特模型(3DMGAR),用于模擬土壤水的三維分布。

        大多數(shù)對(duì)格林-安普特模型的研究集中在下滲計(jì)算,而用其計(jì)算產(chǎn)流的研究相對(duì)較少。目前格林-安普特降雨徑流模型中的下滲曲線為一條拋物線形的經(jīng)驗(yàn)曲線,該曲線能反映流域內(nèi)降雨和下墊面特征分布不均的特點(diǎn),但與實(shí)際下墊面分布情況存在差距。本文通過(guò)流域土壤分布資料和降雨資料,計(jì)算流域內(nèi)各點(diǎn)在不同時(shí)刻的下滲能力,從而得到一條具有物理基礎(chǔ)的下滲能力分布曲線,將該曲線應(yīng)用于格林-安普特模型,探究改進(jìn)的格林-安普特模型在半干旱地區(qū)的應(yīng)用效果。

        1 格林-安普特降雨徑流模型

        格林-安普特降雨徑流模型是基于格林-安普特下滲公式計(jì)算產(chǎn)流的超滲產(chǎn)流模型。在下滲過(guò)程中,假定土壤濕潤(rùn)區(qū)與未濕潤(rùn)區(qū)之間存在一個(gè)水平的濕潤(rùn)鋒,濕潤(rùn)區(qū)土壤含水率為飽和含水率,未濕潤(rùn)區(qū)為初始含水率,隨著下滲過(guò)程的進(jìn)行,濕潤(rùn)鋒不斷地向下移動(dòng)。格林-安普特下滲方程為

        (1)

        式中:f(t)——土壤下滲能力;K——飽和水力傳導(dǎo)度;ψ——濕潤(rùn)鋒處土壤吸力;Δθ——土壤飽和含水率與初始含水率之差;F(t)——累計(jì)下滲量;t——時(shí)間。

        圖1 格林-安普特模型中的下滲能力分布曲線Fig.1 Infiltration distribution curve of Green-Ampt model

        在格林-安普特降雨徑流模型中,產(chǎn)流計(jì)算使用格林-安普特下滲公式,蒸散發(fā)計(jì)算使用新安江模型中的三層土壤水模型[1],河道匯流計(jì)算使用馬斯京根匯流演算法。

        1.1 下滲能力分布曲線

        這條下滲能力分布曲線的表達(dá)式為

        (2)

        式中:BX——分布曲線的指數(shù),控制曲線形狀。

        1.2 改進(jìn)的下滲能力分布曲線

        在新安江模型中,拋物線形的蓄水容量曲線很好地反映了流域蓄水容量分布不均的特點(diǎn)。為驗(yàn)證格林-安普特超滲產(chǎn)流模型中,流域內(nèi)的下滲能力分布是否也呈拋物線形,本文提出一條改進(jìn)的具有物理基礎(chǔ)的下滲能力分布曲線,與原經(jīng)驗(yàn)曲線進(jìn)行對(duì)比,研究哪種曲線更適用于半濕潤(rùn)-半干旱地區(qū)水文模擬。根據(jù)美國(guó)農(nóng)業(yè)部(USDA)土壤分類(lèi)標(biāo)準(zhǔn),Rawls等[17]給出了不同土壤類(lèi)型的格林-安普特下滲模型中的4個(gè)參數(shù)值:土壤總孔隙度、土壤有效孔隙度、濕潤(rùn)鋒處土壤吸力和飽和水力傳導(dǎo)度。利用這些土壤參數(shù)值,即可通過(guò)降雨資料和土壤分布資料計(jì)算出流域內(nèi)各點(diǎn)在任意時(shí)刻的下滲能力,得到一條具有物理基礎(chǔ)的下滲能力分布曲線。具體計(jì)算步驟如下:

        圖2 流域地理位置Fig.2 Location of the Zhidan and Walnut Gulch Basins

        a. 利用土壤類(lèi)型資料和Rawls等給出的下滲參數(shù)值,可由式(1)計(jì)算出流域內(nèi)各點(diǎn)的初始下滲能力,由此得到流域的初始下滲能力分布曲線。

        b. 由流域初始下滲能力分布曲線,結(jié)合降雨資料,計(jì)算得到初始時(shí)刻流域產(chǎn)流量及下滲量,將下滲量帶入式(1)計(jì)算得到下一時(shí)段的下滲能力分布曲線。

        c. 重復(fù)步驟a和b,依次計(jì)算各時(shí)段產(chǎn)流量及下滲能力分布曲線,改進(jìn)的下滲曲線形狀隨降雨和時(shí)間的變化而不斷地變化。

        將這條下滲曲線應(yīng)用于格林-安普特降雨徑流模型,并與原始基于經(jīng)驗(yàn)下滲曲線的格林-安普特降雨徑流模型進(jìn)行對(duì)比。

        2 研究流域與資料選取

        選擇中國(guó)陜西省志丹流域和美國(guó)亞利桑那州Walnut Gulch流域(圖2)作為研究流域,兩流域均為半干旱流域。志丹流域出口控制站為志丹水文站,該站是黃河流域北洛河水系周河控制站。志丹流域面積774 km2,多年平均降水量為509.8 mm,流域內(nèi)植被覆蓋差,水土流失較為嚴(yán)重。洪水多由暴雨形成,漲落較快,洪峰一般呈尖瘦形,洪水歷時(shí)較短。Walnut Gulch流域位于美國(guó)亞利桑那州San Pedro河,流域面積150 km2,流域內(nèi)共有16個(gè)雨量站,多年平均降水量為312 mm,氣候環(huán)境更接近干旱地區(qū),地表植被覆蓋較差,多為低矮的草本植物。志丹流域的流量、降雨量、蒸發(fā)量等資料由當(dāng)?shù)厮木痔峁?,Walnut Gulch流域資料由美國(guó)農(nóng)業(yè)部(USDA)網(wǎng)站提供(https://www.usda.gov/)。

        表1列出了志丹流域和Walnut Gulch流域的土壤類(lèi)型占比,以及各類(lèi)土壤所對(duì)應(yīng)的格林-安普特下滲參數(shù),圖3為兩流域土壤類(lèi)型分布。本研究使用的土壤數(shù)據(jù)為USDA提供的1 km×1 km數(shù)據(jù),在進(jìn)行模型計(jì)算時(shí),將每個(gè)流域劃分成若干子流域,分別計(jì)算各子流域內(nèi)的產(chǎn)匯流過(guò)程。由于子流域的土壤類(lèi)型分布不同,因此改進(jìn)的格林-安普特下滲分布曲線在不同的子流域中形狀也不相同,并且其形狀會(huì)隨著時(shí)間及降雨的變化而改變。

        表1 兩流域土壤類(lèi)型及格林-安普特下滲參數(shù)

        3 模 型 率 定

        將原格林-安普特模型和改進(jìn)的格林-安普特模型模擬結(jié)果分別與新安江模型模擬結(jié)果對(duì)比,研究其在半干旱地區(qū)的適用性。

        圖3 流域土壤類(lèi)型分布Fig.3 Distribution of soil types in the Zhidan and Walnut Gulch Basins

        在3個(gè)模型中,只有改進(jìn)格林-安普特模型中的3個(gè)參數(shù)飽和含水率θS、飽和水力傳導(dǎo)度K和濕潤(rùn)鋒處的土壤吸力ψ是根據(jù)土壤性質(zhì)得出的,其余參數(shù)均使用人工優(yōu)選法和SCE-UA自動(dòng)優(yōu)選法[18-19]結(jié)合率定得到[20-21]。自動(dòng)優(yōu)選法具有高精度和高效率的特點(diǎn),但由于模型中參數(shù)較多,存在異參同效現(xiàn)象,單純使用自動(dòng)優(yōu)選法率定參數(shù)容易陷入局部最優(yōu)解,得不到最佳參數(shù)。而只使用人工優(yōu)選法率定效率較低,精度不高。因此,對(duì)于模型中相對(duì)不敏感的參數(shù),使用人工優(yōu)選法根據(jù)參數(shù)物理意義及流域特點(diǎn)確定參數(shù)值;對(duì)于敏感參數(shù),使用SCE-UA自動(dòng)優(yōu)選法率定參數(shù)值,這樣可以減少自動(dòng)率定參數(shù)的數(shù)量,有效避免異參同效問(wèn)題,提高參數(shù)精度。例如,對(duì)新安江模型中的蓄水容量分布曲線指數(shù)B、深層蒸散發(fā)擴(kuò)散系數(shù)C、地下水出流系數(shù)KG、地下水消退系數(shù)CG等參數(shù)和原格林-安普特模型中的下滲能力分布曲線指數(shù)BX、飽和含水率θS等參數(shù),根據(jù)經(jīng)驗(yàn)及流域氣候特點(diǎn)、下墊面特征人工率定參數(shù)值;對(duì)于新安江模型中蒸散發(fā)折算系數(shù)E、自由水容量SM、河網(wǎng)水流消退系數(shù)CS等和原格林-安普特模型中飽和水力傳導(dǎo)度K、濕潤(rùn)鋒處的土壤吸力ψ等參數(shù)使用SCE-UA方法率定參數(shù)值。

        各模型使用日模型和時(shí)段模型分別模擬,在志丹流域時(shí)段模型計(jì)算步長(zhǎng)為1 h。Walnut Gulch流域面積較小,匯流時(shí)間短,時(shí)段模型計(jì)算步長(zhǎng)選擇10 min。兩流域日模型計(jì)算步長(zhǎng)均為24 h,日模型為時(shí)段模型提供初始土壤含水量等下墊面狀態(tài)初值。

        4 模型模擬結(jié)果比較

        選擇徑流深合格率、洪峰合格率、峰現(xiàn)時(shí)間誤差和確定性系數(shù)4 個(gè)指標(biāo)來(lái)評(píng)價(jià)模型模擬結(jié)果,其中徑流深和洪峰是否合格按照GB/T 22482—2008《水文情報(bào)預(yù)報(bào)規(guī)范》[22]進(jìn)行評(píng)判。對(duì)各場(chǎng)洪水峰現(xiàn)時(shí)間誤差和確定性系數(shù)取平均值,綜合評(píng)價(jià)各模型模擬結(jié)果。

        4.1 志丹流域模擬結(jié)果比較

        在志丹流域選擇2000—2010年共13場(chǎng)洪水,其中前8場(chǎng)洪水用來(lái)率定模型參數(shù),后5場(chǎng)洪水驗(yàn)證模型結(jié)果,3個(gè)模型模擬結(jié)果見(jiàn)表2,模型參數(shù)值見(jiàn)表3。

        表2 志丹流域模擬結(jié)果對(duì)比

        表3 志丹流域各模型參數(shù)值

        對(duì)于徑流深模擬,改進(jìn)格林-安普特模型具有很好的模擬結(jié)果,在率定期和驗(yàn)證期合格率都為100.0%。新安江模型率定期合格率為75.0%,驗(yàn)證期合格率為80.0%,在3個(gè)模型中合格率相對(duì)最低,但仍具有較高的合格率??傮w來(lái)看3個(gè)模型對(duì)徑流深模擬結(jié)果較好,主要由于志丹流域洪水徑流深普遍偏小。根據(jù)《水文情報(bào)預(yù)報(bào)規(guī)范》[22],當(dāng)一場(chǎng)洪水實(shí)測(cè)徑流深小于15 mm時(shí),模擬值與實(shí)測(cè)值誤差在3 mm以內(nèi)為合格。對(duì)于志丹流域的13場(chǎng)洪水,徑流深最大的洪水為2001081612號(hào)洪水,其實(shí)測(cè)徑流深為9.08 mm,另外有7場(chǎng)洪水實(shí)測(cè)徑流深都小于3 mm,因此3個(gè)水文模型對(duì)于徑流深的模擬值與實(shí)測(cè)值誤差大多在3 mm以內(nèi),徑流深合格率較高。對(duì)于志丹流域這樣的半干旱中小流域,其洪水特點(diǎn)為歷時(shí)短、漲落較快、總徑流深小但洪峰流量較大,較大的洪峰可能導(dǎo)致山洪等災(zāi)害發(fā)生,因此對(duì)于洪峰的準(zhǔn)確模擬更為重要。

        在洪峰合格率方面,改進(jìn)格林-安普特模型模擬結(jié)果最好,率定期合格率為50%,驗(yàn)證期合格率為40%,相比原格林-安普特模型有顯著提高。但總體來(lái)看,3個(gè)模型對(duì)洪峰模擬的精度普遍偏低,合格率不高。目前在世界范圍內(nèi),水文模型在濕潤(rùn)及半濕潤(rùn)地區(qū)具有較好的模擬效果,而在干旱和半干旱地區(qū)模擬精度較低,主要由于干旱和半干旱地區(qū)降雨時(shí)空差異大,觀測(cè)資料不足,產(chǎn)流模式較為復(fù)雜。對(duì)于志丹流域,洪峰模擬合格率不高主要是由于降雨觀測(cè)精度不高。大部分雨量站的觀測(cè)時(shí)間步長(zhǎng)為6 h,而實(shí)際產(chǎn)生洪峰過(guò)程的降雨常常發(fā)生在1~2 h內(nèi),將6 h的降雨資料插值為1 h數(shù)據(jù)時(shí)大幅度削弱了降雨強(qiáng)度,造成模型對(duì)洪峰模擬結(jié)果不佳。

        改進(jìn)前后的格林-安普特模型在率定期和驗(yàn)證期具有相同的峰現(xiàn)時(shí)間誤差,并且誤差稍小于新安江模型,這說(shuō)明改進(jìn)格林-安普特模型在模擬洪峰時(shí)間方面與原模型差異不大。改進(jìn)格林-安普特模型確定性系數(shù)比原模型有提高,說(shuō)明改進(jìn)模型對(duì)洪水形狀的模擬與觀測(cè)值更接近。

        總體來(lái)看,改進(jìn)前后的格林-安普特模型在志丹流域模擬結(jié)果均好于新安江模型,主要因?yàn)槟P彤a(chǎn)流機(jī)制不同。新安江模型是蓄滿產(chǎn)流模型,在濕潤(rùn)地區(qū)模擬效果良好,而格林-安普特模型為超滲產(chǎn)流模型,更適用于干旱和半干旱地區(qū)。模擬結(jié)果表明志丹流域產(chǎn)流過(guò)程以超滲產(chǎn)流為主,這與其半干旱流域的特性相符。原格林-安普特模型模擬結(jié)果略好于新安江模型,但優(yōu)勢(shì)并不明顯,主要原因還是降雨數(shù)據(jù)的時(shí)間分辨率不夠。由于格林-安普特模型為超滲產(chǎn)流模型,當(dāng)流域降雨強(qiáng)度大于下滲強(qiáng)度時(shí)產(chǎn)流,否則不產(chǎn)流,所以降雨強(qiáng)度觀測(cè)精度對(duì)格林-安普特模型模擬結(jié)果的影響比對(duì)新安江模型影響更大。另外,志丹流域多年平均降水量為509.8 mm,并不是完全的干旱流域,產(chǎn)流以超滲產(chǎn)流為主,同時(shí)還會(huì)伴隨少量蓄滿產(chǎn)流過(guò)程,因此用格林-安普特模型計(jì)算產(chǎn)流時(shí)會(huì)造成一定的誤差。改進(jìn)格林-安普特模型在徑流深、洪峰和確定性系數(shù)方面模擬結(jié)果均好于原模型,說(shuō)明加入具有物理基礎(chǔ)的下滲能力分布曲線發(fā)揮了作用。這條下滲曲線比原來(lái)的經(jīng)驗(yàn)下滲曲線更好地反映出流域下墊面狀態(tài)隨時(shí)間和空間的變化特點(diǎn)。原始經(jīng)驗(yàn)曲線形狀是固定的,不隨時(shí)間變化,而具有物理基礎(chǔ)的下滲能力分布曲線是基于土壤和降雨數(shù)據(jù)計(jì)算得出的,由于流域上每一點(diǎn)的下滲能力都會(huì)隨時(shí)間和降雨變化而變化,因此改進(jìn)的下滲曲線形狀也是時(shí)刻變化的,這樣更符合實(shí)際情況,也能更精確地計(jì)算出流域內(nèi)的下滲和產(chǎn)流量,使改進(jìn)格林-安普特模型在志丹流域模擬準(zhǔn)確率得到提高。

        圖4展示了格林-安普特模型改進(jìn)前后在志丹流域2010081103號(hào)洪水中野雞岔子流域內(nèi)的下滲能力分布曲線。由于在該子流域內(nèi)主要有兩種不同的土壤類(lèi)型,計(jì)算得到的下滲能力分布曲線呈兩階梯形(圖4(a)(b)),隨著降雨和時(shí)間的進(jìn)行,該基于物理基礎(chǔ)的下滲曲線形狀會(huì)不斷變化。原始的經(jīng)驗(yàn)下滲曲線呈拋物線形(圖4(c)(d)),隨著時(shí)間變化,曲線上各點(diǎn)的下滲能力不斷變化,但曲線形狀是固定不變的。此外,由于各子流域內(nèi)土壤分布情況不同,改進(jìn)的基于物理基礎(chǔ)的下滲曲線在不同子流域內(nèi)具有不同的形狀,而原始經(jīng)驗(yàn)曲線在不同子流域內(nèi)形狀相同。

        圖4 2010081103號(hào)洪水改進(jìn)前后格林-安普特模型下滲能力分布曲線Fig.4 Infiltration distribution curves of Green-Ampt models before and after the improvement for the flood 2010081103

        4.2 Walnut Gulch流域模擬結(jié)果比較

        表4 Walnut Gulch流域模擬結(jié)果對(duì)比

        在Walnut Gulch流域選擇2000—2019年共20場(chǎng)洪水,其中前12場(chǎng)洪水用來(lái)率定模型參數(shù),后8場(chǎng)洪水作為驗(yàn)證。由于實(shí)測(cè)流量資料只有徑流深和洪峰流量值,因此選擇徑流深合格率與洪峰合格率作為評(píng)價(jià)指標(biāo),模擬結(jié)果見(jiàn)表4,模型參數(shù)值見(jiàn)表5。

        在Walnut Gulch流域,新安江模型的模擬結(jié)果不盡人意,尤其對(duì)洪峰模擬誤差較大。由于該流域年均降雨量為312 mm,氣候更接近干旱地區(qū),新安江模型的蓄滿產(chǎn)流理論在該地區(qū)并不適用。

        表5 Walnut Gulch流域各模型參數(shù)值

        原始格林-安普特模型與改進(jìn)格林-安普特模型在率定期和驗(yàn)證期具有相同的徑流深合格率,總體合格率略高于新安江模型。對(duì)于洪峰模擬,原始格林-安普特模型率定期合格率為33.3%(略高于新安江模型),驗(yàn)證期合格率為12.5%(與新安江模型相同)。改進(jìn)格林-安普特模型率定期和驗(yàn)證期的洪峰合格率分別為41.7%和37.5%,均高于原始格林-安普特模型。

        總體來(lái)看,改進(jìn)格林-安普特模型在Walnut Gulch流域模擬效果最好,但洪峰合格率依然較低,主要原因是實(shí)測(cè)資料精度不夠。該流域面積較小,匯流時(shí)間短,洪峰漲落大多在1 h以內(nèi),因此模型模擬時(shí)間步長(zhǎng)選為10 min。而降雨數(shù)據(jù)觀測(cè)時(shí)間間隔較長(zhǎng),甚至超過(guò)1 h,當(dāng)把降雨資料插值為10 min時(shí),會(huì)造成較大的插值誤差,削弱降雨強(qiáng)度。格林-安普特模型為超滲產(chǎn)流模型,降雨強(qiáng)度對(duì)產(chǎn)流計(jì)算影響很大,降雨強(qiáng)度誤差會(huì)帶來(lái)較大的洪水模擬誤差。

        5 結(jié) 語(yǔ)

        將一條具有物理基礎(chǔ)的下滲能力分布曲線引入傳統(tǒng)格林-安普特模型對(duì)其進(jìn)行改進(jìn),并選擇2個(gè)半干旱流域,將改進(jìn)格林-安普特模型與原始格林-安普特模型和新安江模型進(jìn)行對(duì)比。

        研究結(jié)果表明,在志丹流域,改進(jìn)前后的格林-安普特模型模擬效果均比新安江模型更好,并且改進(jìn)格林-安普特模型比原始格林-安普特模型有顯著提高,說(shuō)明引入基于物理基礎(chǔ)的下滲能力分布曲線能更好地反映流域內(nèi)降雨及下墊面變化特征,能夠更加精確地計(jì)算各時(shí)刻流域內(nèi)各個(gè)位置的下滲和產(chǎn)流量。并且由于不同子流域內(nèi)土壤特性不同,改進(jìn)的下滲能力分布曲線在各子流域的形狀也不相同,這比原始的形狀固定的經(jīng)驗(yàn)下滲曲線更符合實(shí)際。在更接近于干旱地區(qū)的Walnut Gulch流域,改進(jìn)格林-安普特模型比原始格林-安普特模型模擬精度提高,但仍然不夠理想。半干旱、干旱地區(qū)洪水多由短時(shí)暴雨引起,陡漲陡落,因此需要更高精度的觀測(cè)資料來(lái)進(jìn)一步提高洪水模擬及預(yù)報(bào)精度。

        在半干旱地區(qū)洪水預(yù)報(bào)中,改進(jìn)格林-安普特模型可以發(fā)揮作用,以后還需要將此模型應(yīng)用于干旱地區(qū),研究其在干旱流域的洪水模擬及預(yù)報(bào)效果。

        猜你喜歡
        志丹普特產(chǎn)流
        產(chǎn)流及其研究進(jìn)展
        生鮮電商平臺(tái)“日日鮮”構(gòu)建研究
        不同坡面單元人工降雨產(chǎn)流試驗(yàn)與分析
        左志丹老師點(diǎn)評(píng)“全國(guó)第十八屆《少兒美術(shù)》杯年度藝術(shù)展評(píng)”獲獎(jiǎng)作品
        “嚴(yán)肅”的“戲仿”——從《普特邁瑟故事集》看辛西婭·奧茲克的“禮拜式”敘事
        探索小日歷
        北京山區(qū)側(cè)柏林地坡面初始產(chǎn)流時(shí)間影響因素
        地表粗糙度對(duì)黃土坡面產(chǎn)流機(jī)制的影響
        Flow regime and head loss in a drip emitter equipped with a labyrinth channel*
        馬來(lái)西亞舉行普特拉賈亞花卉節(jié)
        少妇无码一区二区三区免费| 日韩午夜三级在线视频| 亚洲精品一区二在线观看| 少妇人妻中文久久综合| 日韩日韩日韩日韩日韩| 欧美日韩国产成人高清视频| 久精品国产欧美亚洲色aⅴ大片| 亚洲高清一区二区三区在线观看 | 日本熟妇另类一区二区三区| 国产精品99无码一区二区| 色婷婷综合久久久久中文| 国产精品九九久久一区hh| 亚洲中文字幕诱惑第一页| 国产精品成人av一区二区三区| 中文字幕亚洲欧美在线不卡| 成人免费毛片aaaaaa片| 激情人妻在线视频| 国产av熟女一区二区三区老牛| 日本免费一二三区在线| 日本一区二区在线播放| 国产精品jizz视频| 成人免费无码a毛片| 国产成人一区二区三区| 蜜桃av在线免费网站| 中国女人内谢69xxxx免费视频| 免费无码成人av在线播| 国产人妖一区二区在线| 中文字幕女同人妖熟女| 亚洲av高清在线观看一区二区| 免费观看黄网站在线播放| 99精品国产闺蜜国产在线闺蜜| 亚洲精品综合久久国产二区| 麻豆精品国产av在线网址| 欧美人妻日韩精品| 国产激情一区二区三区| 熟女人妇交换俱乐部| 无码人妻精品一区二区三区下载| 日本一区二区日韩在线| 精品综合一区二区三区| 最近在线更新8中文字幕免费| 中文字幕一区二区人妻出轨|