999精品在线视频,手机成人午夜在线视频,久久不卡国产精品无码,中日无码在线观看,成人av手机在线观看,日韩精品亚洲一区中文字幕,亚洲av无码人妻,四虎国产在线观看 ?

黑河流域高寒草甸生態系統水分收支及蒸散發拆分研究

2018-11-15 05:14:30童雅琴佩1李小雁1張賜成
生態學報 2018年20期
關鍵詞:模型研究

童雅琴,王 佩1,,李小雁1,,*,張賜成,白 巖

1 北京師范大學地表過程與資源生態國家重點實驗室,北京 100875 2 北京師范大學地理科學學部自然資源學院,北京 100875

根據IPCC第五次評估報告,自工業革命以來,全球氣候變化是當前人類面臨的重要議題,其主要表現為氣溫升高和降水格局的變化[1],是改變水循環的重要驅動因素。水分收支是對水循環要素降水、蒸發蒸騰、徑流以及土壤儲水量變化等的定量刻畫[2]。在水分收支過程中,降水是陸地生態系統水分的主要來源,直接影響植被生產力;蒸散發(ET)包括土壤蒸發(E)和植被蒸騰(T)[3],是陸地生態系統水分的主要耗散形式,土壤蒸發和植被蒸騰受不同的生物物理過程控制[4],對于其在生態系統中的相對貢獻的定量研究對準確預測生態系統功能對氣候變化的響應至關重要,當前關于生態系統蒸散發拆分的方法包括實地觀測法[5- 6],同位素法[7- 10]以及模型估算法[11- 13];土壤水分是土壤-植被-大氣連續體的關鍵因子[14- 15],土壤貯水量對維持“土壤水庫”功能,調節降水分布不均等具有重要作用[16]。高寒草甸/草原占中國西部高寒區面積的64%[17],被認為是對氣候變化最為敏感的區域[18- 19]。在氣候變化背景下,準確評估其水分收支過程對理解氣候、植被、水文交互作用至關重要。

黑河流域是典型的干旱區內陸河流域,是以水循環過程為紐帶的冰川/凍土、綠洲和荒漠為主要特征的多元生態系統[20],上游高寒區被認為是干旱區水塔,水循環過程復雜。已有的研究多是從流域尺度運用模型進行水分收支估算,或者側重水分收支的個別要素[21- 24],如Jia等[25]用分布式水文模型分析了黑河中上游水分收支和能量收支平衡,表明在上游地區年均蒸散發和年均徑流分別占年均降水的63%和37%;Qin等[26]表明在黑河上游52.46%的出山徑流量產生于裸地景觀,上游森林景觀幾乎不產流;Yang等[27]運用分布式生態水文模型研究表明,黑河上游地下水徑流和壤中流占主導,地表徑流極少,且水分收支特征受植被類型的影響顯著;以上模型雖然能很好的分析流域尺度水分收支,但是,一般流域尺度上空間異質性較強,缺少在生態系統尺度上對水分收支過程的驗證,且缺少直接觀測的生態系統實際蒸散發,可能導致對生態系統水分以及植被需水量的估算存在偏差。渦動相關技術實地觀測的水汽通量數據,為準確的分析生態系統水分收支提供了可能,且當前運用雙源模型對蒸散發進行拆分來評價生態系統水分收支狀況的研究較少。因此,本文以黑河流域高寒草甸生態系統為研究對象,運用蒸散發拆分模型和渦動相關技術直接觀測的水汽通量數據,估算2014和2015年黑河上游高寒草甸生態系統水分收支及其組分的季節變化和年際變化旨在定量評估黑河流域高寒草甸生態系統的水分收支狀況,為流域水資源承載力評估及水資源管理提供科學依據。

1 研究區概況與數據來源

1.1 研究區概況

黑河是我國第二大內陸河,發源于青海省祁連山區的冰川積雪帶,流域總面積14.3萬km2,是集合高山冰川、森林草原、平原綠洲及戈壁荒漠的復合生態系統[28]。本研究站點位于黑河上游青海省祁連縣中東地區的阿柔鄉,黑河上游支流八寶河南側的河谷高地上,經緯度為100.46°E,38.05°N,海拔3044 m,周圍地勢相對平坦開闊,自東南向西北略有傾斜下降,南北兩側約3km外是連綿的山丘和高山。研究區多年平均降水為400.3 mm[21],降水量年內分配不均,主要集中在生長季(5—9月)。研究區屬于高寒草甸區,植被類型以小嵩草草甸為主。

1.2 數據來源

研究數據包括2014和2015年自動氣象站數據,渦動相關數據及相關的遙感數據,本研究的氣象數據和渦度相關數據來源于阿柔凍融觀測站,黑河綜合遙感聯合試驗:水文氣象觀測網數據集(http://hiwater.westgis.ac.cn/)[29]。其中,渦動相關數據用來計算實際蒸散發和檢驗模型的效果,該數據是經過Edire軟件后處理的30 min通量數據,處理過程包括野點剔除、坐標旋轉、頻率響應修正、WPL修正以及初步質量控制等,本文對缺失的渦動數據運用查表法和平均晝夜法進行插補[30]。氣象數據包括氣溫、降水、相對濕度、太陽輻射、風速、土壤溫濕度(0—320 cm)等[31- 32]。其中,自動氣象站的降水采用翻斗式雨量計觀測,翻斗式雨量計安裝在阿柔超級站40 m觀測塔上。土壤溫濕度探頭埋設在地表0 cm和地下2、4、6、10、15、20、30、40、60、80、120、160、200、240、280、320 cm處,用于土壤水分和溫度的觀測。葉面積指數(LAI)數據來自2014和2015黑河生態水文遙感試驗:黑河流域1km分辨率5d合成葉面積指數(LAI)數據集(http://www.heihedata.org/data)[33- 35],模型中需要的日葉面積指數(LAI)采用三次樣條函數插值所得,植被高度數據參考黑河流域生態水文樣帶調查:2013年上游植被數據[36]。

2 研究方法

2.1 水分收支計算

陸地生態系統水分收支的基本方程[37]為:

P+Sd+I+R上入+R下入=ET+R上出+R下出+D+ΔW

(1)

式中,P為降水量,Sd為地下深層補給量,I為灌溉量,R上入、R下出分別為計算時段內經地面和地下流入生態系統內的徑流量;ET(ET=E+T)為蒸散發,包含有土壤蒸發E和植物蒸騰T,R上出、R下入分別為計算時段內經地面、地下流出的徑流量,ΔW為土壤貯水量的變化量,D為土壤水分的下滲。在祁連山高寒草甸沒有灌溉且地下水埋藏較深,可將I和Sd忽略,本文用于計算土壤貯水量的是0—320 cm處的土壤水分,所以可將D忽略,R上入、R下入和R上出、R下出可分別表示為R入和R出,所以水分收支公式可以簡化為:

P+R入=ET+R出+ΔW

(2)

簡化后公式的物理意義可以理解為到達生態系統的水分消耗于土壤和植被的蒸散發、徑流以及土壤貯水量的變化,若P-ET+ΔW>R出-R入>0表示生態系統水分盈余;反之,表明生態系統水分虧缺。

2.2 土壤貯水量及其變化的估算

土壤貯水量(W)的計算采用各層土壤體積含水量的平均值與其代表的深度乘積累加求得[38],具體計算公式如下:

(3)

2.3 蒸散發組分估算

本文使用Wang和Yamanaka雙源模型[39],該模型以植被冠層和土壤表面能量平衡為基本原理,包括太陽輻射傳輸過程,植被氣孔導度的模擬,而且其在不同的生態系統對蒸散發的模擬與分割都取得了很好的結果[39- 41]。模型的主要公式表征如下:

(4)

(5)

式中,RnV是植被冠層的凈輻射(W/m2),HV是植被冠層的顯熱通量(W/m2),T是蒸騰量 (kg m-2s-1),fV為植被冠層的電介常數,Sd是向下的短波輻射(W/m2),Ld是向下的長波輻射(W/m2),σ為 Stefan-Boltzmann常數(=5.67×10-8W m-2K-4),TG和TL分別為地表溫度和葉片表面溫度(℃),RnG為地表凈輻射(W/m2),G是土壤熱通量(W/m2),HG是地表顯熱通量(W/m2),E是蒸發通量(kg/ m2s),αG為地表反照率,αV是植冠層反照率,αG,αV可近似為常數0.1和0.2。

能量平衡中的顯熱通量可以通過以下等式表示:

HV=ρa(TL-Ta)/raV

(6)

HG=ρa(TG-Ta)/raG

(7)

式中,cp為干燥空氣的熱容量,raV為植被冠層的空氣動力阻抗,raG是地表的空氣動力阻抗(s/m),它們主要與風速,植被高度和風速測量儀的高度有關,詳細描述可見王佩等[39]研究。

模型結果中的土壤蒸發(E)和植被蒸騰計算公式如下:

T=ρa([q(TL)-qa])/(raV+rc)

(8)

E=ρa([q(TG)-qa])/(raG+rSS)

(9)

式中,qsat(TL) 為葉片溫度的飽和比濕 (kg/kg),qsat(TG)為地表溫度的飽和比濕,qa為空氣比濕,raV是植被冠層的空氣動力阻抗,rc為氣孔導度(s/m),以及rSS是土壤表面阻抗(s/m)。

rc=rst/LAI

(10)

式中,rst為葉片氣孔阻抗。

(11)

式中,rst_min和rst_max分別為土壤完全濕潤時葉片氣孔阻力的最大值和最小值,公式中的常數CSW可表示為:

(12)

式中,θ為土壤體積含水量(m3/m3),θmax是土壤最大體積含水量。土壤的表面阻力與土壤含水量滿足指數關系式:

rsoil=a(θs/θ)b+c

(13)

式中,a,b,c為經驗系數,a=3.5,b=2.3 和c=433.5。

2.4 模型的敏感性分析

為了評價模型指定參數和輸入變量可能存在的誤差,對模型進行敏感性分析,敏感性系數si[39,42- 43]定義如下:

(14)

式中,Pi是指影響模型輸出結果(如ET或者T/ET)的i個驅動參數或變量,其偏導?O/(?Pi)計算公式如下:

(15)

3 結果分析

3.1 模型的有效性和敏感性分析

圖1為高寒草甸渦動相關系統觀測的和模型模擬的蒸散發以及地表溫度值的1∶1線性擬合關系,結果發現模型模擬值與實測值之間擬合度高,相關系數為分別0.90和0.89,表明該雙源模型對高寒草甸生態系統蒸散發的模擬是有效的。

圖1 實測值與模型模擬蒸散發(潛熱)值與溫度對比2014—2015Fig.1 Comparison of evapotranspiration (latent heat flux) and temperature between measured and the predicted by model during 2014—2015

表1總結了高寒草甸潛熱通量(LE)和蒸散比(T/ET)對模型輸入變量的敏感性。結果表明潛熱通量(LE)對相對濕度(RH)最為敏感,為-0.74 ± 0.32,即相對濕度5%的誤差會導致潛熱通量平均3.7%的誤差,對于蒸散比T/ET,葉面積指數(LAI)和土壤水分(θ)對其較為敏感,其5%的誤差分別可導致蒸散比T/ET2.2%和1.7%的誤差。因此可以看出模型估算的蒸散發(LE)和蒸散比(T/ET)存在的誤差都在可以忽略的范圍之內,表明模型模擬的結果較為可靠。

表1模型輸入參數對蒸散發(ET)和蒸散比(T/ET)敏感性系數(Si)的平均值±標準偏差(SD)

Table1Meanandstandarddeviation(SD)ofthesensitivitycoefficients(Si)ofevapotranspiration(ET)andtranspirationfraction(T/ET)fortheinputvariablesinmode

變量名稱Variables 平均值±標準差Mean ± Standard deviation潛熱通量 Latent heat Fluex LE蒸散發比 Transpiration fraction T/ET氣溫 Temperature Ta/(℃)0.52±0.380.19±0.27相對濕度 Relative humid RH/%-0.74±0.320.00±0.22葉面積指數 Leaf area index LAI/(m3/m3)0.16±0.090.44±0.28根系層土壤水分 Soil moisture θ/%0.22±0.060.34±0.28

3.2 高寒草甸生態系統水分收支要素特征

3.2.1 降水特征分析

高寒草甸生態系統2014和2015年降水量為分別為520 mm和401 mm,降水量年內分配不均,降水主要集中在5—9月(圖2),2014年和2015年生長季降水分別占全年總降水量的88.75%和90.21%,2014和2015年均溫度分別-0.43℃和0.28℃,氣溫、飽和氣壓差(VPD)以及潛在蒸散發(PET)的季節變化趨勢一致,呈現單峰變化規律,都是從5月開始上升,7月達到最大值,隨后開始緩慢下降。

圖2 2014—2015年研究區氣象要素的季節變化Fig.2 Monthly variation of meteorological variables in study area from 2014 to 2015

圖3 2014和2105年模擬蒸散發與觀測蒸散發對比Fig.3 Comparison of Predicted evapotranspiration and simulated evapotranspiration in 2014 and 2015

3.2.2 蒸散發及其組分變化特征分析

通過渦動相關觀測的高寒草甸生態系統2014年和2015年的實際蒸散發量分別為520 mm和551 mm,生長季(5—9)蒸散發分別為408 mm和422 mm,分別占全年蒸散發總量的78%和76%,圖3顯示為高寒草甸蒸散發的時間序列,從中可以看出觀測蒸散發、模型模擬蒸散發在時間變化上具有較好的一致性。日蒸散發量在全年尺度上呈現明顯的單峰變化規律,4月份蒸散發量開始增加,7月達到全年最大值,2014年和2015年月最大蒸散發量ET分別為85 mm和89 mm,日最大蒸散發量分別為5.30 mm/d和5.24 mm/d,隨后開始下降,這與生長季植被物候變化特征相一致。

已有研究表明植被蒸騰(T)是陸地生態系統水分收支的重要組分,蒸散比(T/ET)可以很好的指示生態系統的生產力。通過模型對蒸散發的拆分發現,在高寒草甸生態系統2014年和2015年蒸散比(T/ET)的日變化也呈現明顯的單峰變化規律(圖4),在非生長季,缺少植被覆蓋,蒸散發(ET)主要是土壤蒸發,蒸散比(T/ET)很小0.3左右;在生長季(5—9月),隨著植被的生長,植被蓋度增加,植被蒸騰在蒸散發中占主導,T/ET迅速增加,2014和2015年生長季平均蒸散比(T/ET)分別為0.74和0.79,日最大值可達0.91。

3.2.3 土壤水分變化特征

2014年和2015年高寒草甸生態系統0—320 cm土壤貯水量季節及年際變化特征如圖5所示,從圖中可以看出該高寒草甸生態系統土壤貯水量存在年內波動,呈現雙峰變化規律。每年3月開始,土層解凍再加上積雪融化,土壤貯水量快速增加,到5—6月份達到每年的最大值,7—8月草甸生態系統具有較強的植被蒸騰,植被的生長消耗了土壤中儲存的水分,使得土壤貯水量略有下降,9—10月蒸騰減弱,土壤貯水量又有所上升,每年11月至次年3月份,受土層凍結的影響,實測土壤含水量僅為土壤中的液態含水量,固態水無法探測,該時段內的土壤貯水量很小,即土壤貯水量受凍融作用影響顯著。圖6為不同深度凍融前后土壤體積含水量的變化,表層至10 cm處土壤體積含水量凍融前后變化顯著,消融后土壤體積含水量明顯高于凍結前,20 cm以下至120 cm處則相反,凍融前后土壤體積含水量變化不大,至160 cm處土壤無凍融,表明冬季凍結對高寒草甸生態系統表層土壤水分保持具有重要作用。

圖4 2014和2015年高寒草甸蒸散比(T/ET)的日變化 Fig.4 The diurnal variation of transpiration fraction T/ET of Alpine meadow ecosystem in 2014 and 2015

圖5 2014和2015年高寒草甸0—320 cm土壤貯水量的日變化 Fig.5 Diurnal variation of soil water storage (0—320 cm) of Alpine meadow in 2014 and 2015

圖6 高寒草甸生態系統不同土層深度凍結前消融后土壤體積含水量Fig.6 Compare to the soil water content between the before freezing and after thawing at different soil depth

在垂直變化上,高寒草甸生態系統土壤體積含水量在0—40 cm處波動較大(圖7),越往表層,土壤含水量越高,40 cm以下土壤體積含水量小(10%以下)且變化不大,說明0—40 cm為高寒草甸生態系統主要的貯水層。從年際變化上看,2014年該高寒草甸生態系0—320 cm處的土壤貯水量大于2015年,表明土壤貯水量的變化深受當年降水量的影響。

圖7 不同深度處的土壤體積含水量的變化Fig.7 Variation of volumetric soil water content at different depths

3.3 高寒草甸生態系統水分收支特征

圖8為2014和2015年高寒草甸生態系統水分收支要素降水(P),蒸散發(ET)和土壤貯水量變化(ΔW)的月變化及水分收支狀況,每年10月至次年2月,觀測到的降水量與蒸散發量很小,土壤貯水量的變化由于冬季土層凍結,使(P-ET-ΔW)為正值;在每年3—5月,春季氣溫回升,再加上降水少,該時段凍結的土層及冬季降雪相繼融化,補充土壤水分,土壤儲水量急劇增加,導致(P-ET-ΔW)為負值,可見,水分收支的季節變化受凍融過程影響顯著;在完全不受凍融影響的月份6—9月(表2),2014年降水較多,生態系統有地表徑流形成,6—9月地表徑流量分別為6、15、0.45和22 mm;2015年降水量少,6—8月份,蒸散發不僅消耗全部的降水,還消耗儲存在土壤中的水分,生態系統呈現水分虧缺,6—8月虧缺量分別為22、12、13 mm,9月份到植被生長的末期,植被蒸騰減弱,該時段生態系統產流量為7 mm。

圖8 2014和2015高寒草甸生態系統水分收支月變化Fig.8 Monthly variation of water budget in Alpine meadow ecosystem in 2014 and 2015

從表3可以看到,在全年尺度上,2014年和2015年高寒草甸呈現不同程度的水分虧缺,2014年降水多,虧缺量為18 mm,2015年虧缺量為134 mm;在完全不受凍融影響的月份,2014年該高寒草甸生態系統產生地表徑流量為42 mm;2015年降水量少,生態系統呈現為水分虧缺,總的虧缺量為26 mm。

表2 阿柔高寒草甸2014至2015年凍融起止時間

根據不同深度土壤最高溫度(Tmax)和最低溫度(Tmin),秋季始凍期(Tmax>0℃且Tmin<0℃)、完全凍結期(Tmax<0 ℃且Tmin<0℃)和春季解凍期(Tmax>0℃且Tmin<0℃),完全解凍(Tmax>0℃且Tmin>0℃)。

表3 2014和2015高寒草甸生態系統水分收支要素年變化

4 討論

4.1 蒸散比T/ET及影響因素分析

圖9 葉面積指數(LAI)與蒸散比(T/ET)的關系 Fig.9 Relationship of LAI with daily mean transpiration fraction (T/ET)

通過對蒸散發ET的拆分發現在高寒草甸生態系統植被蒸騰,尤其生長季植被蒸騰是生態系統水分消耗的主要形式,2014和2015年高寒草甸生態系統年均蒸散比(T/ET)分別為0.54和0.53。Wei等[44]研究表明在全球尺度上T/ET為0.57 ± 0.06。Hu等[4]基于S-W雙源模型研究的2004年高寒草甸的T/ET為0.44,王海波等[45]基于Penman-Monteith模型研究發現2008—2009年高寒草甸T/ET約為0.6。本研究估算T/ET與其他相關研究基本接近,但也存在一定的差異,可能與LAI數據集有關,已有研究表明T/ET在季節尺度上主要受植被葉面積指數(LAI)的影響[4,39,44,46- 48],這與我們的研究結果相一致(圖9),葉面積指數LAI與T/ET呈現明顯的對數關系。Wei等[44]表明,不同的來源的LAI數據集,可能會導致T/ET存在一定的誤差,因此,高精度的LAI數據對蒸散發ET拆分影響重大。

4.2 水分收支估算

本研究利用水量平衡方程對高寒草甸生態系統2014和2015年水分收支進行估算。在每年的非生長季,各月蒸散發量略大于降水量,生態系統水分虧缺,這與Yao等[49]運用2004—2011實測數據對海北高寒草甸生態系統水分收支研究結果相一致。在2014年不受凍融影響的生長季,生態系統水分盈余。在正常年份的生長季,生態系統水分呈現虧缺,蒸發和蒸騰不僅消耗全部的降水,還消耗儲存在土壤中的水分,儲存在土壤中的水分對生態系統水分收支起重要調節作用。Yan等[50]通過對森林生態系統的研究發現,土壤水分是旱季蒸散發的重要水源,這與本研究結果較為一致。在年尺度上,阿柔站高寒草甸生態系統并無地表徑流產生,只有在降水量大的濕潤年份的生長季,有地表徑流形成。Yang等[51]運用分布式水文模型對黑河上游水量平衡研究表明2001—2012年,地下水徑流占總徑流量的80%左右,壤中流占20%左右,地表徑流量極小,占徑流量的2%—3%。陳仁升等[17]在黑上游高寒地區通過長時間小流域觀測實驗表明,高山寒漠帶為山區流域的主要產流區,而高寒草甸/草原區徑流貢獻很少,其水源涵養功能大于水文功能;通常平緩地形下的高寒草甸/草原具有較多的水源涵養功能,對徑流的貢獻少[17],這與本研究結果相當。

在高寒草甸生態系統,土壤水分的劇烈變化主要發生在0—40 cm處,與郭淑海[45]在天山高寒草甸的試驗結果一致,0—40 cm土壤貯水,是蒸散發消耗的重要水源,對維持植被生長、發育具有重要作用。此外,在生態系統水分收支中,外界徑流補給是生態系統水分消耗的重要來源,而且已有研究[52- 53]表明高寒草甸生態系統地下水埋藏較深,且壤中流較為發育,已有相關研究[26,54-55]表明黑河上游,受凍土的影響,流域產流以壤中流為主,因此,高寒草甸生態系統虧缺的水分可能是壤中流補給。

5 結論

本文運用雙源模型對黑河上游高寒草甸生態系統蒸散發組分進行拆分及以生態系統水量平衡方程為原理,估算了黑河流域高寒草甸生態系統2014年和2015年水分收支狀況,得到以下結論:(1)生長季植被蒸騰是生態系統水分主要的消耗形式,在季節尺度上,蒸散比(T/ET)主要受葉面積指數(LAI)的控制;(2)0—40 cm土壤儲水對維持生長季旺期植被蒸騰具有重要的意義;(3)在濕潤年(2014年)高寒草甸生態系統水分收支基本平衡,且生長季有產流;在正常年份(2015年)生態系統呈現水分虧缺,虧缺的水分主要來自外界壤中流的補給。總體而言,實測數據可以準確的刻畫生態系統水分收支過程,但也存在一定的不確定性,如儀器的觀測精度和觀測年限等,因此,在未來的研究中應加強對生態系統的連續觀測,為遙感數據和模型模擬研究提供更精確地驗證數據。

猜你喜歡
模型研究
一半模型
FMS與YBT相關性的實證研究
2020年國內翻譯研究述評
遼代千人邑研究述論
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
視錯覺在平面設計中的應用與研究
科技傳播(2019年22期)2020-01-14 03:06:54
EMA伺服控制系統研究
新版C-NCAP側面碰撞假人損傷研究
3D打印中的模型分割與打包
主站蜘蛛池模板: 超清无码熟妇人妻AV在线绿巨人| 国产精品免费p区| 国产一区成人| 亚洲伊人久久精品影院| 女人爽到高潮免费视频大全| 伊人久久婷婷五月综合97色| 国产无码制服丝袜| 91成人精品视频| 亚洲国产天堂久久九九九| 尤物视频一区| 992tv国产人成在线观看| 高清无码不卡视频| 国产中文在线亚洲精品官网| 精品国产一二三区| 蜜桃臀无码内射一区二区三区 | 亚洲天堂精品视频| 99久久精品美女高潮喷水| 99久久精品久久久久久婷婷| 亚洲精品手机在线| 国产视频一区二区在线观看| 国产成人啪视频一区二区三区| 国产玖玖视频| 欧美日韩一区二区在线免费观看| 亚洲一区免费看| 亚洲欧洲天堂色AV| 国产中文一区二区苍井空| 四虎免费视频网站| aⅴ免费在线观看| 国产欧美专区在线观看| 一级毛片基地| 激情乱人伦| 日本黄色a视频| 动漫精品中文字幕无码| 一级毛片免费不卡在线| 亚洲欧美日韩另类| 97国产在线观看| 老熟妇喷水一区二区三区| 成人蜜桃网| 国产精品免费p区| 中文一级毛片| 91久久偷偷做嫩草影院精品| 国产精品手机在线观看你懂的| 婷婷亚洲视频| 91精品啪在线观看国产91| 国产伦精品一区二区三区视频优播 | 國產尤物AV尤物在線觀看| 一级黄色网站在线免费看| 欧美曰批视频免费播放免费| 一区二区三区在线不卡免费| 91色爱欧美精品www| 色成人综合| 熟妇人妻无乱码中文字幕真矢织江| 在线看免费无码av天堂的| 婷婷六月天激情| 久久semm亚洲国产| 国产性爱网站| 精品久久久无码专区中文字幕| 男女性色大片免费网站| 中国一级特黄视频| a级毛片网| 99ri国产在线| 免费va国产在线观看| 99精品久久精品| 欧美亚洲一二三区| 伊大人香蕉久久网欧美| 国产精品成人不卡在线观看 | 国产91色| 中文字幕首页系列人妻| 久久99这里精品8国产| 亚洲AV无码久久精品色欲 | 人妻中文字幕无码久久一区| 2021天堂在线亚洲精品专区| 欧美成人午夜在线全部免费| 久久午夜夜伦鲁鲁片无码免费| 91久久精品日日躁夜夜躁欧美| 91久久国产综合精品女同我| 亚洲无码高清一区| 大陆国产精品视频| 欧美一级特黄aaaaaa在线看片| 在线欧美日韩| 国产一区二区三区夜色| 成年免费在线观看|