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

鳳陽山針闊混交林通量觀測源區分布及特征

2020-11-13 02:01:38紀小芳魯建兵莊家堯葉立新劉勝龍方萬力何雪凱
生態學報 2020年20期
關鍵詞:大氣模型研究

紀小芳,龔 元,鄭 翔,魯建兵,馮 明,莊家堯,葉立新,劉勝龍,方萬力,王 丹,何雪凱,姜 姜,*

1 南方現代林業協同創新中心,南京林業大學林學院,南京 210037 2 南京林業大學生物與環境學院,南京 210037 3 鳳陽山-百山祖國家級自然保護區鳳陽山管理處,龍泉 323700

在全球氣候變化的背景下,2002—2011年間CO2濃度的增加速率達到了(2.0±0.1 )ppm/a,其對溫室效應的貢獻大約78%[1]。森林生態系統作為主要的碳匯之一,在陸地生態系統中扮演著重要的角色。過去幾十年大量的研究圍繞森林生態系統的碳循環展開,試圖了解和掌握其碳排放動態,為分析和解決溫室效應等全球環境問題尋找途徑。

渦度協方差(Eddy covariance,EC)方法常用來測量生態系統尺度上地表與大氣的CO2氣體交換。渦度協方差也稱渦度相關法是一種多用于測量植被下墊面與大氣間熱量、物質和動量交換的工具,目前為國際上主流的基于微氣象理論的通量觀測技術,20世紀90年代以來,隨著渦度相關法的發展和應用為直接觀測不同生態系統CO2排放和吸收提供了技術方法[2]。早期的渦度相關技術多應用于農田、濕地、草地和森林等幾個關鍵陸地生態系統與大氣之間的CO2交換研究[3],但受觀測高度,大氣邊界層高度,大氣穩定度等環境因素的影響,通量塔傳感器所測得的通量值僅代表下墊面的一定區域[4],因此在應用渦度協方差法時,需要確定其空間代表性和通量源區[5]。目前通量源區的模型主要有解析模型、大渦模擬模型、拉格朗日模型和閉合模型等幾類。解析模型的假設基礎較多,理論上僅適用于下墊面平緩的區域,主要采用梯度擴散理論、二維平流擴散方程以及相似理論得出通量貢獻區,較著名的有Kormann and Meixner(KM)模型、FSAM模型和Horst-Weil模型[6- 9];大渦模型最初用于大氣和環境科學的研究,但其物理機理、計算比較復雜繁瑣,耗時,消耗大量的存儲空間,不適宜于長期通量觀測數據下的計算[10- 11];拉格朗日模型是基于拉格朗日粒子擴散的數值模擬,理論上嚴格考慮了擴散的均勻分布約束,可以正確的反映出非均勻湍流中的擴散,較為著名的有Hsieh模型以及Kljun模型[12- 13]。2015年由Kljun[14]等提出的三維通量足跡模型更加適用于長期連續且多時間序列的通量塔觀測源區計算,運算速度更快。目前應用FSAM模型[5,9,15- 23]和KM模型[24- 27]的研究較多,但是對Kljun模型[28- 29]實際進行運用的研究較少。對通量源區的主要研究結果發現生長季的源區分布在任何狀態下均小于非生長季[5,24- 25,30- 31];但是也有研究表明在大氣不穩定狀態下時,非生長季的源區面積大于生長季[32];在穩定大氣條件下的通量源區長度要顯著大于不穩定條件下[33];不同時間其通量源區的大小有所不同[16],隨著研究的深入與實際需要針對下墊面不均一的森林生態系統碳通量的監測方法也逐漸成為研究的熱點[34]。Schmid H P 等[35]對位于尼日爾薩赫勒地區的虎灌木(粗砂和灌木交錯分布)進行研究發現從單一位置進行測量不能代表生態系統的空間湍流通量。Sogachev A等[36]利用SCADIS模型對山脊不同位置的通量值進行研究發現,相比較于同等高度,背山脊處的通量所受到的干擾最大,而通量塔的最佳建設位置在山脊頂部。而目前使用的最多的為ART (Agroscope Reckenholz Tanikon) Footprint Tool模型(基于KM模型)來實現對空間異質性不同土地類型對所測通量所占比例的確定[12,27]。

本研究利用Kljun通量足跡模型(后簡稱Kljun模型)和ART Footprint Tool(KM模型),對浙江鳳陽山通量塔2017年全年的觀測數據進行分析,探討該森林生態系統在不同時間跨度、不同大氣條件下的通量源區變化情況以及不同土地類型的對通量觀測值的貢獻,對研究區通量觀測值的空間代表性作出解釋,為其他森林生態系統CO2通量研究和類似下墊面的通量源區評估提供服務和參考。

1 材料與方法

1.1 研究區概況

研究區位于中國東南部的浙江麗水龍泉市鳳陽山自然保護區內,該區建于1975年,現有管理面積15171hm2[37- 38]。位于東經119°06′—119°15′、北緯27°46′—27°58′。保護區森林覆蓋率達到了90.8%,屬于中亞熱帶溫暖濕潤氣候區,同時受海洋性氣候和季風影響較大,研究區內雨量充沛,多年平均降水量達到2400mm,降水集中于4—6月,占全年的80%,濕度大,霧多;年均氣溫12.8℃,極端高低溫分別為30.2℃、-12.5℃,月均氣溫最高28℃左右,最低6℃—13℃,年蒸發量達到1170 mm以上,無霜期275 d,有效積溫約6500℃[39- 40]。四季分明,一般3月底4月初入春;7月初入夏;8月中旬入秋;11月中旬入冬,光照資源豐富[41]。土壤類型主要為黃壤土,分布在海拔800 m之上的高山坡地。地形以山地為主,地勢較高,溝壑交錯,保護區的黃茅尖,海拔高達1928 m,為江浙第一高峰。保護區內的通量塔于2016年10月搭建完成,位于研究區中心,其所在下墊面為陰坡24°,森林遍山野嶺,風浪區無限大,樹冠高度大部分分布均勻,符合Kljun模型的使用條件[42],周圍林地以木荷、杉木和木楠混交為主,平均冠層高度約15 m,平均林齡約40 a。植物資源豐富,形成多優勢種結構特征。人工杉木林和人工柳杉林分布在海拔1400—1500 m左右,木荷在海拔300 m,900 m和1500 m處均有分布,黃山松在海拔600 m和1300 m處均有分布,多脈青岡和黃山木蘭在海拔1500 m左右有分布[43]。

圖1 鳳陽山研究區土地利用類型地形圖Fig.1 Topographic map and land use types in Fengyangshan research area

1.2 觀測方法和數據處理

研究區中心建有高40m的觀測鐵塔,塔上裝有通量和梯度系統。通量觀測系統的傳感器有數據采集器(CSI CR300,Campbell Scientific Inc.,USA),三維超聲及CO2/H2O分析儀(IRGASON,Campbell Scientific Inc.,USA),四分量輻射表(CNR4,Campbell Scientific Inc.,USA),空氣溫濕度(HMP155A,Campbell Scientific Inc.,USA),土壤熱通量板(HFP01,Campbell Scientific Inc.,USA)。配套的氣象觀測系統為6層梯度觀測(空間高度2、8、16、24、32、40m,土壤深度10、20、30、40、60、90cm),梯度系統測量大氣中不同高度的溫度、濕度、風速、風向以及不同深度土壤含水量和土壤溫度,適用于不同的下墊面和大氣條件,是邊界層氣象、農林氣象、大氣環境監測最普遍運用和最基本的觀測手段。主要傳感器有土壤溫度傳感器(TCAV)、土壤水分傳感器(Campbell CS616,Campbell Scientific Inc.,USA)、風向傳感器(MetOne 020C,Campbell Scientific Inc.,USA)、數據采集器(CSI CR300,Campbell Scientific Inc.,USA)、6層風速傳感器(MetOne 010C,Campbell Scientific Inc.,USA)、HOBO雨量傳感器(Onset RG3-M,Campbell Scientific Inc.,USA)和空氣溫濕度傳感器(HMP155A,Campbell Scientific Inc.,USA)。整套觀測系統于2016年10月修建完成。研究人員每隔三個月去現場進行數據下載和維護儀器等工作。本研究選取鳳陽山通量塔2017全年的連續通量觀測數據。采集到的數據經過在線程序(EasyFluxTM-PC, Campbell Scientific Instruments, USA, https://www.campbellsci.com/)處理,EasyFlux程序中所用代碼均與EddyPro (Li-COR,NE, USA, https://www.licor.com/)軟件所用源代碼一致。主要處理過程包括:野點去除、坐標旋轉、WPL校正、建立數據質量等級指標(1—9)等最后得到30 min時間間隔的通量和微氣象數據序列[44- 45]。之后對數據進行篩選,主要是去除降雨及降雨前后1h數據;剔除紅外分析儀信號強度低于0.8的數據;去除質量控制等級標記為“9”的數據[46]。基于本文剔除數據的原則,2017年的通量數據有效率為51%,來進行通量貢獻區分析,一般通量數據有效率超過50%則具有代表性[26,28]。

依據鳳陽山當地的實際氣候狀況,季節劃分以3、4、5月為春季,6、7、8為夏季,9、10、11為秋季,1、2、12月為冬季[47]。

1.3 通量足跡模型

一般通量觀測塔上紅外分析儀所測得的通量值為某一時刻迎風方向上對觀測值產生影響的下墊面空間代表區域的碳源或碳匯強度,對觀測點通量值有貢獻的下墊面區域即為通量貢獻區[48- 49]。以觀測點通量塔為原點(0,0)建立坐標,x軸的正值方向代表上風距離:

(1)

式中Zm為有效觀測點高度,f為碳源或碳匯的轉換函數即footprint函數,代表了表面上某點(x,y)對Zm處觀測值的貢獻率密度,Qc代表的是表面碳源或碳匯的點源強度,Fc為在Zm處所測得的通量值,R代表對通量值有貢獻的下墊面區域[7]。

Zm=Zreceptor-Zd

(2)

式中Zreceptor為通量塔相對于地面的觀測高度,Zd為零平面位移高度,采用森林冠層高度的2/3[50- 51]。

而側風積分足跡函數fiy則主要受到上風距離(x),有效動力學高度(Zm),大氣邊界層高度(h),摩擦風速(u*)和垂直風速脈動的標準差(σw)等幾個參數決定。以上參數除了h需要單獨計算,主要與L(Obukhov,奧布霍夫長度)有關,詳細計算方法Kljun文章中給出[14],其他的都可以通過EasyFluxTM在線程序處理獲取[46]。由量綱分析(Π定理),將以上變量有可能構成的無量綱參數組如下所示:

Π1=fiyZm

(3)

(4)

(5)

(6)

側風積分足跡的無量綱函數F*可以由為上風距離X*表示,即F*=φ(X*),由X*=Π2Π3-1Π4,F*=Π1Π3-1Π4可得到:

(7)

式中的a,b,c,d為通過后向拉格朗日隨機粒子擴散模型(LPDM-B)不斷的實驗計算出的各擬合參數,均與粗糙度(z0)有關。用Zm/L來判別大氣穩定度,其中L為Obukhov長度(可由EasyFlux在線程序計算獲取),當Zm/L>0時,大氣為穩定狀態,而當Zm/L<0時,大氣為不穩定狀態[21]。一般在分析中常用貢獻率密度的P水平等值線所包圍的區域(一般取P=0.8或0.9)來表示EC系統的觀測范圍[21],這里參考Kim等[42]基于Kljun模型的森林生態系統通量源區的研究方法,本研究選取P=0.8[42],Kljun等[42]還提供了該模型的在線數據處理網站(http://www.footprint.kljun.net/index.php),該研究下載了該模型的開源Matlab (MathWorks,USA,https://www.mathworks.com/)函數代碼應用于浙江鳳陽山針闊混交林森林生態系統的通量足跡研究。

2 結果與分析

2.1 主風向分析

從圖2和表1可以清晰的看到鳳陽山森林全年和四季的風速和風向分布頻率,除秋季的主風向為東北風之外,其余季節與全年的主風向均為西南風。結合圖1中的地形可以看出,風向受到地形的很大影響,由于通量塔位于山坡上,所以其西北方向的風受到山坡高度的影響,很難傳輸到通量塔所在位置,而西南方向海拔從通量塔位置開始逐漸增高,因此不會完全阻擋風,給風的傳輸提供了可能,東北和西南方向也是如此。同時從表1中的白天和黑夜不同風頻數據可以看出,全年晝夜風向相反,白天和黑夜的主風向分別為東北風和西南風,可能受到山谷風的影響[52]。

圖2 研究區風向風速玫瑰圖Fig.2 The wind rose map in the study area

表1 不同時間不同風向的風頻

2.2 研究區總體通量源區特征

利用Kljun模型計算出30 min時間間隔的90%的通量貢獻率的通量足跡,并繪制出其所對應的上風距離分布圖和累計頻率概率圖。從圖中可以明顯的看出,通量貢獻區最遠點的分布比較集中,距離通量塔觀測點最遠的距離為7000m,主要的通量測量值來自于東北和西南方向,結合全年的風向風速圖2可以發現,風向和風速是影響通量貢獻區的主要因素。

對全年的貢獻率在10—80%的源區分布進行分析發現,在不穩定條件下的源區面積明顯小于穩定條件。當大氣處于不穩定條件時,其通量貢獻率在80%的源區面積為0.89 km2,源區長度在323.62—839.62 m間,0—90°方向上源區長度達到最大,此時大氣做垂直劇烈運動,物質不穩定輸送速度很快,渦流更易形成和傳播,EC系統所測得的通量值來自上風向更短距離,導致其源區面積較小;而當大氣處于穩定條件時,其通量貢獻率在80%的源區面積為4.88 km2,源區長度在672.69—2390.49 m之間,在180—270°方向上源區長度達到最大,空氣的湍流運動較弱,物質上下垂直擴散速度比較緩慢,EC系統可以測得來自較遠距離的渦流,因此所代表的源區面積較大[16,24,53]。

2.3 不同季節的通量源區足跡氣候態特征

圖6表示在40m觀測高度Kljun模型評價浙江鳳陽山針闊混交林分別在春、夏、秋、冬季節在大氣條件穩定和不穩定條件通量貢獻率為80%時的源區分布圖。可以看出,通量貢獻區存在季節差異,不同季節通量貢獻區的面積和主要分布方向都有差異,穩定條件下的源區面積要明顯大于不穩定大氣條件下,與全年的源區分布規律相似,且源區的方向和形狀均與風向風速圖相一致。春季在不穩定和穩定條件的源區面積分別為0.64 km2和4.89 km2,源區長度分別在277.08—707.41 m和622.73—1887.58 m之間,分別在0—90°和180—270°方向源區長度達到最大。夏季在不穩定和穩定條件下的源區面積分別為1.07 km2和2.58 km2,源區長度分別在442.04—1822.53 m和329.40—1880.01 m之間,源區長度最大值均出現在180—270°方向。秋季在不穩定和穩定條件下的源區面積分別為1.42 km2和5.37 km2,源區長度分別在359.69—1212.98 m和441.73—2717.48 m之間,源區長度最大值出現在0—90°和180—270°風向。冬季不穩定和穩定條件下的源區面積分別為0.70 km2和5.54 km2,源區長度分別在122.43—756.86 m和650.90—2468.66 m之間,源區長度最大值風向為0—90°和180—270°風向。明顯可以看出,源區長度的最大值出現方向主要在東北(0—90°)和西南(180—270°)方向,這完全符合全年源區的分布規律。在大氣條件不穩定的條件下,源區長度分布不超過2000 m,源區范圍從大到小排列為:秋季、夏季、冬季和春季,可以看出生長季的貢獻區在大氣不穩定條件下較大;當大氣條件穩定時,源區長度分布不超過3000 m,源區范圍從大到小排列為:冬季、秋季、春季和夏季,此時生長季的源區分布最小。

圖3 通量足跡最遠點分布圖Fig.3 The distribution of the flux footprint furthest point

圖4 通量足跡最遠點累計頻率圖Fig.4 The probability chart of the flux footprint furthest point

圖5 全年源區分布Fig.5 Flux footprint through 2017(0,0)點為觀測點,東南西北方向分別為x正軸、y負軸、x負軸和y正軸

圖6 不同季節的源區分布Fig.6 Flux footprint during different seasons

圖7 側風積分函數圖Fig.7 Crosswind integral function diagram

從圖7中可以看出通量貢獻率峰值均位于傳感器附近,不超過150m,且隨著大氣條件從不穩定到穩定,其在上風向的傾斜更少。在不穩定條件下,其湍流強度較高導致化合物向上傳輸和較短的傳輸距離與傳輸時間。在穩定條件下側風向距離是明顯高于不穩定條件下,該結論與Kim[42]等研究結果一致,與前文通量貢獻區在不穩定條件下的源區分布范圍小于大氣穩定條件的結果一致。Zm/L<0,表示大氣不穩定狀態,當Zm/L>0,表示大氣穩定狀態[22]。從表2可以看出在春、夏、秋、冬季,當大氣條件不穩定時,Z0、h、sigmav和u*值均大于穩定條件下的,在穩定條件下,夏季的z0、h、sigma和u*最高,冬季的L值最高;在不穩定條件下,夏季的z0、sigmav、u*和L的絕對值最高,冬季的h最高。

2.4 不同土地利用類型的通量貢獻率

下墊面的土地利用類型對通量值有重要影響,所以能將兩者進行結合,可以對源區不同植被對通量觀測值的貢獻有更加直觀的認識。本研究選取1120×1180 m的研究區范圍利用ART Footprint Tool[6,54- 56]進行定量化分析,將該范圍分割成3304個20×20m的方格,確定其格內的土地利用類型并求出其對通量值的貢獻率,最后疊加得到不同植被類型在全年對通量貢獻的百分占比(圖8)所選范圍全年總貢獻率達90%。從圖中可以看出源區貢獻從大到小依次為針闊混交林、闊葉林、建筑用地及道路、杉木林、毛竹林、柳杉林、黃山松林。結合圖1可以看出通常土地類型的面積越大,其占比就較大。毛竹林的面積是柳杉林的兩倍,但是其貢獻百分比卻是柳杉林的3.7倍,可以看出毛竹林在對該研究區的碳預算功能有重要作用。該區域內的建筑用地主要為生態定位站的實驗樓及標準氣象場,其本身的CO2通量理論上可以忽略不計。

表2 側風積分函數輸入參數

圖8 不同土地利用類型的源區貢獻率占比Fig.8 Proportion of source area contribution by land use type

3 討論

通量貢獻區主要受到風向和風速的影響,但是自然條件中的風向風速幾乎時刻在發生變化,為了更好的對貢獻區進行分析,本文分析了不同大氣層結條件和不同季節下的貢獻區分布。如結果(圖6)中所示,通量源區的位置在不同季節有相應的變化,主要受到風向風速的影響,因為夏季的風向風速變化幅度較大,所以其對應的源區分布范圍無論是在大氣穩定還是不穩定的條件下,都分布較廣,冬季時的大氣較為穩定從圖6中的冬b可以看出,該季節在西南方向的源區分布范圍明顯較其他季節大,而在大氣穩定的條件下夏季的源區范圍較小,主要是夏季植物茂盛生長,導致湍流形成比較迅速,從而導致源區的范圍較小,該結論與朱明佳[53]等的研究結果一致,在大氣不穩定的條件下源區范圍從大到小排列為:夏季、秋季、春季和冬季,可以看出生長季的貢獻區在大氣不穩定條件下較大,該結論與龔笑飛[16]等研究結論一致;無論從全年還是季節的尺度上來看,在大氣不穩定的情況下,其源區分布范圍均小于大氣穩定條件,主要是由于當大氣條件不穩定時,空氣上下交換較快,垂直運動強烈,擴散也較快,最終使通量值所代表的源區范圍較小,然而在大氣穩定條件下,空氣垂直產生湍流的速度減緩,擴散減弱,源區范圍可以追溯到上風方向更遠處,該結論與前人[5,16,24- 25,31- 32,57- 58]研究結論均一致。本研究區的源區距離最遠可達到7000m(圖3),這可能也與下墊面、冠層高度和儀器的觀測高度有關,本研究中觀測高度為40m,冠層高度大約15m,超過了冠層高度的2.5倍。與魏遠[5]等人所得出源區最遠點距離為3500m結果不同,主要其儀器高度為25m,冠層高度為16m,由于儀器測量高度不同從而導致源區最遠點的距離有差異;而顧永劍[17]等對濕地生態系統進行源區分析得出其源區長度不超過300m,主要由于濕地的植被均為較低矮的草類植被,儀器觀測高度僅為4.8m,因此造成其源區長度小于本研究。而下墊面會對摩擦速率產生影響,從而對通量貢獻區產生影響,在未來的研究中,特別是在地形及其復雜的情況下,探索研究出能夠考慮到下墊面,地形變化等方面的源區模型并加以運用。本研究發現建筑用地及道路的源區分布貢獻率也會達到5.48%,僅次于常綠闊葉林。雖然這部分對CO2通量理論上可以忽略不計,但人為活動產生的CO2排放會低估對該區域的森林碳匯功能。

4 結論

本文利用2017年全年浙江鳳陽山針闊混交林生態系統通量塔觀測資料,在保證數據質量的前提下,對其所測通量值所代表的源區分布進行分析,從而對研究區通量觀測值的空間代表性作出解釋,為其他森林生態系統CO2通量研究和類似下墊面的通量源區評估提供服務和參考,研究發現:

研究區內全年盛行東北風和西南風,風向在東北方向的占總風頻的34.16%,風向在西南方向的占總風頻的46.96%;而風向是影響源區分布的主要因素,研究區的源區分布主要也在東北和西南方向,從源區的季節分布也可以看出,各季節的源區分布方向和形狀基本與該季節的風向一致;當通量貢獻率達到90%和80%時,通量觀測值分別來自距通量塔7000m和3000m范圍,且各季節的側風向積分函數峰值均位于傳感器附近,不超過150m。在大氣層結不穩定條件下,湍流強度較高;穩定時,源區范圍大于大氣不穩定條件,主要由于大氣不穩定條件直接影響了大氣運動的強度,湍流形成速度與傳播距離;季節尺度上,當大氣條件不穩定時,源區范圍從大到小排列為:秋季、夏季、冬季和春季;當大氣條件穩定時源區范圍從大到小排列為:冬季、秋季、春季和夏季源區貢獻從大到小依次為針闊混交林、闊葉林、建筑用地及道路、杉木林、毛竹林、柳杉林、黃山松林。

猜你喜歡
大氣模型研究
一半模型
大氣的呵護
軍事文摘(2023年10期)2023-06-09 09:15:06
FMS與YBT相關性的實證研究
遼代千人邑研究述論
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
視錯覺在平面設計中的應用與研究
科技傳播(2019年22期)2020-01-14 03:06:54
EMA伺服控制系統研究
3D打印中的模型分割與打包
大氣古樸揮灑自如
主站蜘蛛池模板: 亚洲AV电影不卡在线观看| 久久久精品无码一二三区| 日本a∨在线观看| 亚洲国产中文综合专区在| 在线播放国产99re| 国产成人资源| 亚洲人成网站日本片| 欧美午夜在线观看| 亚洲首页在线观看| 国产乱子伦视频在线播放| 9啪在线视频| 国产91精品久久| 国产精品免费电影| 国产欧美日韩专区发布| 国产靠逼视频| 97色婷婷成人综合在线观看| 999国内精品久久免费视频| 无码'专区第一页| 亚洲欧美h| 在线观看国产小视频| 久青草免费视频| 澳门av无码| 思思99热精品在线| 全部无卡免费的毛片在线看| 狼友av永久网站免费观看| 日韩区欧美区| 日韩无码视频播放| 日本a∨在线观看| 亚洲不卡av中文在线| a级毛片免费播放| 亚洲男女天堂| 国产福利大秀91| 久久这里只有精品2| 国产福利在线观看精品| 日韩欧美91| 国模视频一区二区| 91福利一区二区三区| 91青青草视频| 亚洲精品国产综合99久久夜夜嗨| 欧美激情视频二区三区| 日韩在线2020专区| 久久频这里精品99香蕉久网址| 国产靠逼视频| 国产精品欧美在线观看| 日韩高清无码免费| 欧美国产三级| 国产超薄肉色丝袜网站| 伊人蕉久影院| 多人乱p欧美在线观看| 久久综合九九亚洲一区| 狠狠做深爱婷婷久久一区| 99热这里都是国产精品| 日韩无码视频网站| 永久免费精品视频| 亚洲侵犯无码网址在线观看| 狠狠色噜噜狠狠狠狠奇米777| 久久精品这里只有精99品| 国产尤物在线播放| 老汉色老汉首页a亚洲| 色网站在线视频| jizz在线观看| 国产毛片久久国产| AV无码无在线观看免费| 69精品在线观看| 伊人久热这里只有精品视频99| 中文字幕日韩视频欧美一区| 亚洲国产综合精品一区| 国产精品99在线观看| 久久精品最新免费国产成人| 亚洲欧美另类久久久精品播放的| 欧美日韩高清在线| 免费无码AV片在线观看国产 | 欧洲av毛片| 老熟妇喷水一区二区三区| 国产日本欧美在线观看| 欧美精品影院| 精品亚洲欧美中文字幕在线看| 小13箩利洗澡无码视频免费网站| 国产jizzjizz视频| 成人精品视频一区二区在线| 波多野结衣中文字幕久久| 国产精品男人的天堂|