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

Validation of WRF model on simulating forcing data for Heihe River Basin

2011-12-09 09:36:12XiaoDuoPanXinLi
Sciences in Cold and Arid Regions 2011年4期

XiaoDuo Pan , Xin Li

1. Laboratory of Remote Sensing and Geospatial Science, Cold and Arid Regions Environmental and Engineering Research Institute, Chinese Academy of Sciences, Lanzhou, Gansu 730000, China

2. Graduate University of Chinese Academy of Sciences, Beijing 100049, China

Validation of WRF model on simulating forcing data for Heihe River Basin

XiaoDuo Pan1,2*, Xin Li1

1. Laboratory of Remote Sensing and Geospatial Science, Cold and Arid Regions Environmental and Engineering Research Institute, Chinese Academy of Sciences, Lanzhou, Gansu 730000, China

2. Graduate University of Chinese Academy of Sciences, Beijing 100049, China

The research of coupling WRF (Weather Research and Forecasting Model) with a land surface model is enhanced to explore the interaction of the atmosphere and land surface; however, regional applicability of WRF model is questioned. In order to do the validation of WRF model on simulating forcing data for the Heihe River Basin, daily meteorological observation data from 15 stations of CMA (China Meteorological Administration) and hourly meteorological observation data from seven sites of WATER(Watershed Airborne Telemetry Experimental Research) are used to compare with WRF simulations, with a time range of a whole year for 2008. Results show that the average MBE (Mean Bias Error) of daily 2-m surface temperature, surface pressure, 2-m relative humidity and 10-m wind speed were -0.19 °C, -4.49 hPa, 4.08% and 0.92 m/s, the average RMSE (Root Mean Square Error)of them were 2.11 °C, 5.37 hPa, 9.55% and 1.73 m/s, and the averageR(correlation coefficient) of them were 0.99, 0.98, 0.80 and 0.55, respectively. The average MBE of hourly 2-m surface temperature, surface pressure, 2-m relative humidity, 10-m wind speed,downward shortwave radiation and downward longwave were -0.16 °C, -6.62 hPa, -5.14%, 0.26 m/s, 33.0 W/m2and -6.44 W/m2, the average RMSE of them were 2.62 °C, 17.10 hPa, 20.71%, 2.46 m/s, 152.9 W/m2and 53.5 W/m2, and the averageRof them were 0.96, 0.97, 0.70, 0.26, 0.91 and 0.60, respectively. Thus, the following conclusions were obtained: (1) regardless of daily or hourly validation, WRF model simulations of 2-m surface temperature, surface pressure and relative humidity are more reliable, especially for 2-m surface air temperature and surface pressure, the values of MBE were small andRwere more than 0.96;(2) the WRF simulating downward shortwave radiation was relatively good, the averageRbetween WRF simulation and hourly observation data was above 0.9, and the averageRof downward longwave radiation was 0.6; (3) both wind speed and rainfall simulated from WRF model did not agree well with observation data.

forcing data; weather research and forecasting model; watershed airborne telemetry experimental research; Heihe River Basin

1. Introduction

The near-surface atmospheric elements including air temperature, pressure, relative humility, wind, precipitation and radiation are called forcing data to drive hydrological,land surface, and ecological models (Liet al., 2007). Cosgroveet al. (2003) emphasized the importance of forcing data, and pointed out that the land surface model relies heavily on these data to predict the land surface state when he introduced the North American land surface data assimilation system. Coupling or use of the General Circulation Model (GCM) output to drive the land surface model is an important source of forcing data for large-scale land surface models (Milleret al., 1994; Amel, 1999; Nijssenet al.,2001), however the spatial and temporal resolution can not meet the requirements of a regional scale land surface model and its coarse resolution ignored regional environmental impact assessment (Wilbyet al., 1999). Thus, a regional scale land surface model or land surface data assimilation system needs forcing data with a higher temporal-spatial resolution. Therefore, downscaling GCM output is necessary in order to prepare for a long series of high spatial and temporal resolution forcing data for regional scale or basin-scale land surface, hydrological and ecological simulation and assimilation.

Downscaling is a method of obtaining high-resolution climate or climate change information from relatively coarse-resolution GCMs (Wilby and Wigley, 2000). There are two main approaches for deriving information on local or regional scales from the global climate scenarios generated by GCMs (Wilby and Wigley, 1997, 2000; Wilbyet al., 2002): (1) numerical downscaling (also known as "dynamical downscaling") involving a nested Regional Climate Model (RCM) or (2) empirical statistical downscaling employing statistical relationships between the large-scale climatic state and local variations derived from historical data records. The major advantage of statistical downscaling is saving computation resources (von Storchet al., 2000; Fanet al., 2005, 2007). However, statistical downscaling requires considerable amount of observation data to derive the transfer function between the observed small-scale variables and the larger scale variables, and it is difficult to apply statistical downscaling for mountain areas such as much of western China. Also, there is no theoretical basis for its application when there is no significant correlation between large scale climatic factors and regional scale factors (Mearnset al., 1999). Dynamic downscaling uses large scale climatic background information to establish higher resolution meteorological distribution fields by GCMs (Giorgi and Mearns, 1999). Regional climate model is considered to be an effective dynamic downscaling method to derive high-resolution local climate information (Liet al., 2009). The model’s physical meaning is clear, it is built on the nonlinear interaction between large scale and regional scale, and can distinguish local features characterized by local topography, vegetation cover and soil hydrological processes (Giorgi and Bates,1989; Honget al., 1999), and can be used anywhere independent of observation data with an improved description of local features. With the rapid development of computer hardware and software, regional climate models have been widely promoted.

Kunstmann and Stadler (2005) and Kunstmannet al.(2008) coupled Mesoscale Model 5 (MM5) into a distributed hydrological model called Water Balance Simulation Model (WaSIM) for high resolution runoff simulation in the catchment of the Mangfall River and developed a decision support system in the Volta Basin, respectively. Jasperet al.(2002) compared the performance of five different high-resolution numerical weather prediction models (with resolutions between 2×2 km2and 14×14 km2) for the prediction of peak flows in the alpine Ticino-Toce watershed.Yuet al. (2002) linked a Hydrologic Model System (HMS)to RCM, which was designed to provide fine spatiotemporal output for hydrologic and other applications, to simulate a series of storm events passing over the Susquehanna River Basin, and to simulate various hydrologic processes in soil,land surface, and ground-water hydrology using observed and modeled storm events.

It was convenient to couple RCM into a land surface model, and the results from the coupled model were validated by observed land surface states, however, the validation of forcing data was not sufficient enough, especially for high temporal resolution observed data (Baoet al.,2006). Thus, the validation accuracy of forcing data from RCM is important before input to drive the land surface model. Due to the hourly observed data from WATER Project (detail introduction about WATER see Liet al.,2009), the forcing data from Weather Research and Forecasting Model (WRF) is validated by daily and hourly observed data in our study.

2. Data and WRF model

2.1. Study area

This research was carried out in the Heihe River Basin,China’s second largest inland river basin, which is located between 97°24′E-102°10′E and 37°41′N-42°42′N, and covers an area of approximately 140,000 km2(Figure 1).Landscapes are diverse, distributed from upper stream to down stream: glacier, frozen soil, alpine meadow, forest,irrigated crops, riparian ecosystem, desert, and gobi. The Heihe River Basin is located in the central Eurasian continent, far from the sea, surrounded by high mountains; its climate is mainly affected by the middle-high latitude westerlies and polar cold air mass circulation. Main climate characters are dry, scarce and concentrated precipitation,strong winds, abundant sunshine, strong solar radiation, and a large temperature difference between day and night.

2.2. Data

The model was validated at 15 surface meteorological stations maintained by the China Meteorological Administration (CMA) and seven observation-strengthened stations setup by WATER project. Table 1 lists the geographical characters at these stations. Data from CMA stations are used for daily validation; the WATER’s are for hourly validation. All CMA data in 2008 are used for comparison, the time span of WATER’s are listed in Table 2.

2.3. WRF model and its configuration and initialization

The WRF (Michalakeset al., 1998, 2001) model is a next-generation mesoscale numerical weather prediction system that serves both operational and research communities. This system consists of multiple dynamical cores,preprocessors for producing initial and lateral boundary conditions for simulations. WRF is built using software tools to enable extensibility and efficient computational parallelism. The use of WRF system has been reported in a variety of areas including storm prediction and research,air-quality modeling, wildfire, hurricane, tropical storm prediction, and regional climate and weather prediction(Michalakeset al., 2005).

Figure 1 Heihe River Basin

Table 1 Geographic features at validation stations

In this study, WRF is used for downscaling of weather and climate ranging from one kilometer to thousands of kilometers and used for deriving meteorological parameters required for hydrological models. The model uses terrain-following hydrostatic pressure coordinate system with permitted vertical grid stretching (Laprise, 1992). Arakawa-C grid staggering is used for horizontal discretization.The model equations are conservative for scalar variables.The detailed description of WRF is presented in Skamarocket al. (2008). For this study, two way nested computational domains of 41×55×28 and 101×121×28 grid points and horizontal resolutions of 25 km and 5 km, respectively, have been set. The first domain covers most of the Gansu Province, ranging from 32.6°N to 47.4°N in latitude and 92.4°E-107.6°E in longitude (Figure 2). The second domain covers the Heihe River Basin ranging from 37°N to 43°N in latitude and 96.6°E-103.4°E in longitude. The model is initialized by real boundary conditions using NCAR-NCEP’s Final Analysis (FNL) data (NCEP-DSS083.2, 2009) having a resolution of 1°×1° (111 km×111 km). A ratio of 1:5 is maintained between resolutions of the outer domain and FNL data to ensure reliable boundary conditions for the model. The WRF simulations have been carried out on the Dell R900, Ubuntu 9.10, g95 compiler with gcc.

Table 2 Time span of WATER station (whole year means year 2008)

Figure 2 Nesting domain configuration for numerical experiment (blue cross means CMA stations, red dot means WATER stations)

3. Results

WRF model simulation was validated with daily CMA and hourly WATER observation data for 2-m surface temperature, surface pressure, 2-m relative humidity and 10-m wind speed, were validated with hourly WATER observation data only for downward shortwave radiation and downward longwave radiation.

Figures 3, 4, and 5 indicate the daily and hourly 2-m temperature validation for WRF model in the Heihe River Basin. Results show that WRF simulation agrees well with observation data.Rvalue of daily validation between WRF simulation and each CMA station is more than 0.97, the absolute of MBE for every station is less than 2 °C,and MBE of Ejin, Mazongshan and Alashan are near zero withRmore than 0.98.Rvalue of hourly validation between WRF simulation and each WATER site is more than 0.96.

Figure 3 Daily validation of WRF 2-m air temperature

Figure 4 Hourly validation of WRF 2-m air temperature

Figure 5 Comparison of daily surface air temperature between WRF simulating and observed data in 2008 in Ejin station

Figures 6 and 7, with the same 2-m temperature, show that surface pressure from WRF model agrees well with observation data not only for daily but also for hourly validation.Rvalue between WRF simulation and each CMA station is more than 0.95, the absolute of MBE for every station is less than 20 hPa, and there are 6 stations (total is 15) whose MBE are less than 2 hPa. Also,Rvalue between WRF simulation and WATER sites are more than 0.95.There is an interval in theX-axis of hourly validation for AR site, where the observation data was checked, and found that when the pressure value is more than 700 hPa, the value type is "integer", though others are "float". Thus, the data in AR is not wrong but recorded in a coarse way.

Figure 8 shows that 2-m relative humidity of WRF simulation is higher than most of CMA stations’ observation data, but Figure 9 does not show the same pattern at all, theRvalues for daily validation are around 0.80, and are around 0.68 for hourly validation.

WRF simulation of 10-m wind speed yields better results in daily timescale than hourly ones (Figures 10 and 11),the average RMSE of daily 10-m wind speed is around 1.73 m/s, which is better than hourly’s (around 2.46 m/s).

There is no daily observation data for downward shortwave radiation and downward longwave radiation. Scatterplots of downward shortwave radiation show a better convergence than downward longwave radiation, especially for GT, ML and YK stations (Figures 12 and 13).

4. Discussion

Daily and hourly validation is summarized in Tables 3 and 4. The average MBE of daily 2-m surface temperature,surface pressure, 2-m relative humidity and 10-m wind speed were -0.19 °C, -4.49 hPa, 4.08% and 0.92 m/s; the average RMSE were 0.11 °C, 5.37 hPa, 9.55 % and 1.73 m/s,and the averageRwere 0.99, 0.98, 0.80 and 0.55, respectively.The average MBE (not including HZZ station, it is discussed later) of hourly 2-m temperature, surface pressure, 2-m relative humidity, 10-m wind speed, downward shortwave radiation and downward longwave radiation were -0.16 °C,-6.62 hPa, -5.14%, 0.26 m/s, 32.98 W/m2and -6.44 W/m2;the average RMSE were 2.62 °C, 17.10 hPa, 20.71%, 2.46 m/s, 152.911 W/m2and 53.53 W/m2; and the averageRwere 0.96, 0.97, 0.70, 0.26, 0.91 and 0.60, respectively.

4.1. Error discussion

HZZ station is excluded in the hourly validation scatter plots because of its larger error than others (Table 4). A reasonable explanation is that the large error is mainly caused by imprecise terrain information description in WRF mode,so more precise terrain data was used instead of WRF’s.Results show that WRF simulation with better terrain description agrees well with observation in HZZ (Figure 14).Enhanced detail about the terrain replacing procession will be discussed in another paper.

4.2. WRF simulation via reanalysis

NCEP II and Princeton data were used to compare with WRF simulation. NCEP II and Princeton date were statistically downscaled by assuming a constant temperature lapse rate (-0.65 °C/100 m) and a constant pressure lapse rate(-10 Pa/100 m), and assuming relative humidity does not change with elevation. Results (Figures 15 and 16) show that WRF simulation of 2-m temperature and surface pressure agree better with observation than NCEP II and Princeton data, but it is not obvious for 2-m relative humidity.

Figure 6 Daily validation of WRF air Pressure

Figure 7 Hourly validation of WRF air pressure

Figure 8 Daily validation of WRF 2-m relative humility

Figure 9 Hourly validation of WRF 2-m relative humility

Figure 10 Daily validation of WRF 10-m wind speed

Figure 11 Hourly validation of WRF 10-m wind speed

Figure 12 Hourly validation of downward shortwave (SWD) flux at ground surface

Figure 13 Hourly validation of downward longwave flux (LWD) at ground surface

?

Figure 14 Comparison analysis of WRF simulation between pre- and post-terrain replacing procession in June in WATER HZZ station(WRF-1 means pre-procession, WRF-2 means post-procession, T2: 2-m temperature, PSFC: surface pressure, RH2: 2-m relative humidity,WS10: 10-m wind speed, SWD: downward shortwave radiation, LWD: downward longwave radiation)

Figure 15 Comparison curves of 2-m temperature, surface pressure and 2-m relative humidity

Figure 16 MBE, RMSE and R between observed data and NCEP reanalysis, Princeton reanalysis, WRF simulating in WATER stations(a) for 2-m temperature; (b) for surface pressure; (c) for 2-m relative humidity.

5. Conclusions

Based on the WRF simulation daily and hourly validation from 15 CMA stations and 7 WATER stations observation, the following conclusions were obtained: (1) regardless of daily or hourly validation, WRF model simulations of 2-m surface temperature, surface pressure and relative humidity are more reliable, especially for 2-m surface air temperature and surface pressure, the MBE were small and theRs were more than 0.96; (2) the WRF simulating downward shortwave radiation was relatively good, the averageRbetween WRF simulation and hourly observation data was above 0.9, the averageRof downward longwave radiation was 0.6; (3) both wind speed and rainfall simulated from WRF model did not agree well with observation data.

This work was supported by grant from the National High Technology Research and Development Program (863) of China (Grant No.2009AA122104) and grants from the National Natural Science Foundation of China (No.40901202,No.40925004). The data used in this paper were obtained from the WATER experiment, which is jointly supported by the CAS Action Plan for West Development Program (Grant No.KZCX2-XB2-09) and Chinese State Key Basic Research Project (Grant No.2007CB714400). The input data for WRF model are from the Research Data Archive (RDA)which are maintained by the Computational and Information Systems Laboratory (CISL) at the National Center for Atmospheric Research (NCAR). NCAR is sponsored by the National Science Foundation (NSF). The original data are available from the RDA (http://dss.ucar.edu) in dataset number ds083.2. The authors thank two anonymous reviewers and the editor for their very helpful comments.

Amel N, 1999. The effect of climate change on hydrological regimes in Europe:A continental perspective. Global Environmental Change, 9: 5-23.

Bao Y, Lü SH, Lu DR, Hou RQ, 2006. Application of regional climate model(RegCM3) in Northwest China, I: Simulation of an arid extreme event.Journal of Glaciology and geocryology, 28(2): 164-175.

Cosgrove BA, Lohmann D, Mitchell KE, Houser PR, Wood EF, Schaake JC,Robock A, Marshall C, Sheffield J, Duan QY, Luo LF, Higgins RW, Pinker RT, Tarpley JD, Meng J, 2003. Real-time and retrospective forcing in the North American Land Data Assimilation System (NLDAS) project. Journal of Geophysical Research, 108(D22): 8842.

Fan LJ, Fu ZB, Chen DL, 2007. Estimation of local temperature change scenarios in North China using statistical downscaling method. Chinese Journal of Atmospheric Sciences, 31(5): 887-897.

Fan LJ, Fu ZB, Chen DL, 2005. Review on creating future climatic change scenarios by statistical downscaling techniques. Advances in Earth Sciences,20(3): 320-329.

Giorgi F, Bates T, 1989. The climatological skill of a regional model over complex terrain. Monthly Weather Review, 117: 2325-2347.

Giorgi F, Mearns LO, 1999. Introduction to special section: Regional climate modeling revisited. Journal of Geophysical Research, 104: 6335-6352.

Hong SY, Juang HMH, Lee DK, 1999. Evaluation of a regional spectral model for the East Asian monsoon case studies for July 1987 and 1988. Journal of Meteorology. Society of Japan, 77: 553-572.

Jasper K, Gurtz J, Lang H, 2002. Advanced flood forecasting in Alpine watersheds by coupling meteorological observations and forecasts with a distributed hydrological model. Journal of Hydrology, 267: 40-52.

Kunstmann H, Jung G, Wagner S, Clottey H, 2008. Integration of atmospheric sciences and hydrology for the development of decision support systems in sustainable water management. Physics and Chemistry of the Earth, 33:165-174.

Kunstmann H, Stadler C, 2005. High resolution distributed atmospheric-hydrologic modeling for Alpine catchments. Journal of Hydrology,314: 105-124.

Laprise R, 1992. The Euler equation of motion with hydrostatic pressure as independent coordinate. Monthly Weather Review, 120(1): 197-207.

Li X, Huang CL, Che T, 2007. Progress and prospect of Chinese land data assimilation system. Progress in Natural Science, 17(2): 163-173.

Li X, Li XW, Li ZY, Ma MG, Wang J, Xiao Q, Liu Q, Che T, Chen EX, Yan GJ,Hu ZY, Zhang LX, Chu RZ, Su PX, Liu QH, Liu SM, Wang JD, Niu Z,Chen Y, Jin R, Wang WZ, Ran YH, Xin XZ, Ren HZ, 2009. Watershed allied telemetry experimental research. Journal of Geophysical Research,114(D22103).

Mearns LO, Bogardi I, Giorgi F, Matyasovszky I, Palecki M, 1999. Comparison of climate change scenarios generated from regional climate model experiments and statistical downscaling. Journal of Geophysical Research,104(D6): 6603-6621.

Michalakes J, Chen S, Dudhia J, Hart L, Klemp J, Middlecoff J, Skamarock W,2001. Development of a Next Generation Regional Weather Research and Forecast Model. In: Zwieflhofer W, Kreitz N (eds.). Developments in Teracomputing: Proceedings of the Ninth ECMWF Workshop on the Use of High Performance Computing in Meteorology. World Scientific, Singapore:269-276.

Michalakes J, Dudhia J, Gill D, Henderson T, Klemp J, Skamarock W, Wang W,2005. The Weather Research and Forecast Model: Software Architecture and Performance. In: Zwiefhofer W, Mozdzynski G (eds.), Proceedings of the 11th ECMWF Workshop on the Use of High Performance Computing in Meteorology. World Scientific, Reading U.K.: 156-168.

Michalakes J, Dudhia J, Gill D, Klemp J, Skamarock W, 1998. Design of a next-generation regional weather research and forecast model. In: Zwieflhofer W, Kreitz N (eds.), 1999. Towards Teracomputing: Proceedings of the Eighth ECMWF Workshop on the Use of Parallel Processors in Meteorology. World Scientific, River Edge, New Jersey: 117-124.

Miller J, Russell G, Caliri G, 1994. Continenta1 scale river flow in climate models. Journal of Climate, 7: 914-928.

Nijssen B, O’Donnell GM, Lettenmaier DP, 2001. Predicting the discharge of global rivers. Journal of Climate, 14: 3307-3323.

Skamarock WC, Klemp JB, Dudhia J, Gill DO, Barker DM, Duda MG, Huang XY, Wang W, Powers JG, 2008. A description of the advanced research WRF Version 3. www.mmm.ucar.edu/wrf/users/docs/user_guide/ARWUsersGuide.pdf.

University Corporation for Atmospheric Research, National Center for Atmospheric Research, 2009. http://dss.ucar.edu/datasets/ds083.2/.

von Storch H, Hewitson B, Mearns L, 2000. Review of empirical downscaling techniques. In: Iversen T, H?iskar BAK (eds.). Regional climate development under global warming, Conference Proceedings RegClim Spring Meeting. General Technical Reports, 4: 29-46.

Wilby RL, Dawson CW, Barrow EM, 2002. SDSM-a decision support tool for the assessment of regional climate change impact. Environmental Modeling and Software, 17: 145-157.

Wilby RL, Hay LE, Leavesley GH, 1999. A comparison of down—scaled and raw GCM output:Implications for climate change scenarios in the San Juan River Basin, Colorado. Journal of Hydro1ogy, 225: 67-91.

Wilby RL, Wigley TM, 1997. Downscaling general circulation model output: a review of methods and limitations. Progress in Physical Geography, 21:530-548.

Wilby RL, Wigley TM, 2000. Precipitation predictors for downscaling: Observed and general circulation model relationships. International Journal of Climatology, 20(5): 641-661.

Yu Z, Barron EJ, Yarnal B, Lakhtakia MN, White RA, Plollard D, Miller DA,2002. Evaluation of basin-scale hydrologic response to a multi-storm simulation. Journal of Hydrology, 257: 212-225.

10.3724/SP.J.1226.2011.00344

*Correspondence to: Dr. XiaoDuo Pan, Assistant Professor of Cold and Arid Regions Environmental and Engineering Research Institute, Chinese Academy of Sciences. No.320, West Donggang Road, Lanzhou, Gansu 730000, China. Tel: +86-931-4967961;Email: panxiaoduo@lzb.ac.cn

23 January 2011 Accepted: 12 April 2011

主站蜘蛛池模板: 国内自拍久第一页| 亚洲欧美一区二区三区图片| 欧美成在线视频| 尤物国产在线| 欧美在线一级片| 91色在线观看| 在线观看国产一区二区三区99| 国产乱人乱偷精品视频a人人澡| 亚洲欧美h| 国产哺乳奶水91在线播放| 高清免费毛片| 青草视频久久| 呦女精品网站| 色网站在线免费观看| 韩日免费小视频| 欧美成人国产| 亚洲AV无码乱码在线观看代蜜桃| 欧美三級片黃色三級片黃色1| 日韩精品一区二区三区中文无码| 一本大道东京热无码av| 91av成人日本不卡三区| 成年片色大黄全免费网站久久| a级毛片在线免费| 国产人成网线在线播放va| 亚洲大尺度在线| 97成人在线视频| 欧美va亚洲va香蕉在线| 91最新精品视频发布页| 欧美日本二区| 国产地址二永久伊甸园| 欧美综合区自拍亚洲综合天堂| 九色91在线视频| 91美女视频在线观看| 中字无码精油按摩中出视频| 国产精品自拍露脸视频| 亚洲妓女综合网995久久| 中文字幕有乳无码| 99国产精品一区二区| 欧美区国产区| 无码一区中文字幕| 99在线视频免费| 久青草免费在线视频| 免费在线播放毛片| 女同国产精品一区二区| 欧美日韩动态图| 精品人妻AV区| 国产精品无码在线看| 尤物亚洲最大AV无码网站| 一级片免费网站| 日韩成人午夜| 97在线免费| 国产凹凸一区在线观看视频| AV老司机AV天堂| 欧美伊人色综合久久天天| a级毛片在线免费| 日韩欧美国产中文| 国产高潮流白浆视频| 国产视频 第一页| 亚洲乱伦视频| 在线观看国产精美视频| 免费A级毛片无码免费视频| 真人高潮娇喘嗯啊在线观看 | 国产成人精品亚洲日本对白优播| 四虎免费视频网站| 国产精品福利在线观看无码卡| 婷婷色婷婷| 亚洲精品桃花岛av在线| 伊人AV天堂| 91丝袜美腿高跟国产极品老师| 91小视频在线| 色综合久久久久8天国| 国产成人综合亚洲欧美在| 国产激情在线视频| 精品日韩亚洲欧美高清a| 97在线公开视频| 亚洲av无码片一区二区三区| 亚洲精品国产日韩无码AV永久免费网 | 久久视精品| 精品91视频| 人妻少妇乱子伦精品无码专区毛片| 三上悠亚精品二区在线观看| 免费视频在线2021入口|