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

    Estimation of Turbulent Sensible Heat and Momentum Fluxes over a Heterogeneous Urban Area Using a Large Aperture Scintillometer

    2015-06-09 21:37:22SangHyunLEEJunHoLEEandBoYoungKIM
    Advances in Atmospheric Sciences 2015年8期

    Sang-Hyun LEE,Jun-Ho LEE,and Bo-Young KIM

    1Department of Atmospheric Science,Kongju National University,Gongju 314-701,South Korea

    2Department of Optical Engineering,Kongju National University,Gongju 314-701,South Korea

    Estimation of Turbulent Sensible Heat and Momentum Fluxes over a Heterogeneous Urban Area Using a Large Aperture Scintillometer

    Sang-Hyun LEE?1,Jun-Ho LEE2,and Bo-Young KIM2

    1Department of Atmospheric Science,Kongju National University,Gongju 314-701,South Korea

    2Department of Optical Engineering,Kongju National University,Gongju 314-701,South Korea

    The accurate determination of surface-layer turbulent fluxes over urban areas is critical to understanding urban boundary layer(UBL)evolution.In this study,a remote-sensing technique using a large aperture scintillometer(LAS)was investigated to estimate surface-layerturbulentfluxes overa highly heterogeneous urban area.The LAS system,with an opticalpath length of 2.1 km,was deployed in an urban area characterized by a complicated land-use mix(residential houses,water body,bare ground,etc.).The turbulent sensible heat(QH)and momentum fluxes(τ)were estimated from the scintillation measurements obtained from the LAS system during the cold season.Three-dimensional LAS footprint modeling was introduced to identify the source areas(“footprint”)of the estimated turbulent fluxes.The analysis results showed that the LAS-derived turbulent fluxes for the highly heterogeneous urban area revealed reasonable temporal variation during daytime on clear days,in comparison to the land-surface process-resolving numerical modeling.A series of sensitivity tests indicated that the overall uncertainty in the LAS-derived daytime QHwas within 20%–30%in terms of the influence of input parameters and the nondimensional similarity function for the temperature structure function parameter,while the estimation errors inτwere less sensitive to the factors of influence,except aerodynamic roughness length.The 3D LAS footprint modeling characterized the source areas of the LAS-derived turbulent fluxes in the heterogeneous urban area,revealing that the representative spatial scales of the LAS system deployed with the 2.1 km optical path distance ranged from 0.2 to 2 km2(a“micro-αscale”), depending on local meteorological conditions.

    heterogeneous surface,large aperture scintillometer,mesoscale modeling,surface energy balance,turbulent exchange

    1.Introduction

    Turbulent transfer of heat,moisture,and momentum continuously occurs at the interface between the atmosphere and underlying Earth surface.The atmospheric boundary layer (ABL)is formed as a consequence of the interactions on time scales of a few hours or less(Garratt,1992;Foken,2008). Therefore,accurate determination of the turbulent fluxes in the surface layer is of primary importance to understanding the spatial and temporal variability of atmospheric flow and thermal structure through the ABL,which sequentially influences the formation of the characteristic local climate and meteorology,as well as air quality.

    Generally,eddy-covariance(EC)systems and scintillometers have been widely used to determine surface-layer turbulent fluxes.The EC system measures surface-layer turbulentfluxeson the orderofhundredsofmetersorless,which has been successfully used to provide the local-scale turbulent fluxes for various natural and urban surfaces(e.g.,Offerle et al.,2006).The EC system directly measures atmospheric turbulent fluctuations,and thus the EC technique has great merit in producing accurate turbulent fluxes for relatively homogeneous flat surfaces.However,it is difficult to determine the sensor’s location formeasuring turbulentfluxes over highly heterogeneous surfaces,due to the uncertainty of surface source areas of turbulent fluxes observed at a particular location.Unlike the EC technique,the scintillometry technique uses Monin–Obukhov similarity theory to estimate area-averaged surface-layer turbulent fluxes.Thus,the turbulent fluxes can be determined at various spatial scales ranging from 100 m to a few kilometers,depending on the scintillometer type and optical path distance.For example,a small aperture scintillometer(SAS)uses a short optical path length (a few hundred meters),while a large aperture scintillometer (LAS)has a longer path length of up to 5 km.

    Many field experiments based on the scintillometry technique have been conducted to determine surface-layer turbulent fluxes over natural surfaces(e.g.,De Bruin et al.,1995, 2002;Chehbouni et al.,2000;Beyrich et al.,2002;Zeweldi et al.,2010;Evans et al.,2012)and urban areas(e.g.,Kanda et al.,2002;Lagouarde et al.,2006;Pauscher,2010;Ward et al.,2014),and showed the usefulness of LAS-derived turbulent fluxes through direct comparison with EC measurements(e.g.,McAneney et al.,1995;Hartogensis et al.,2003; Lagouarde et al.,2006;Zeweldi et al.,2010;Liu et al.,2011, 2013;Xu et al.,2013;Ward et al.,2014).For example, Lagouarde et al.(2006)showed the LAS-derived turbulent fluxes with a path length of 1.8 km compared well to the EC fluxes over the city of Marseille,France.Ward et al.(2014) compared two LAS systems with different path lengths of 2.8 and 5.5 km and an EC instrument for surface-layer turbulent flux estimation over a residential area in Swindon,UK.These studies demonstrated the applicability of the LAS system over homogeneous residential areas in determining surfacelayer turbulent fluxes in the urban environment.

    Most urban areas at spatial scales of a few km or less are,morphologically,highly heterogeneous due to the diversity of buildings,natural/artificial land-use mix,and topography.Even though previous studies have successfully determined the surface-layer turbulent fluxes over homogeneous urban areas using the EC system and/or LAS system(e.g., Lagouarde et al.,2006;Ward et al.,2014),there are few studies concerning turbulent flux estimation over heterogeneous urban surfaces.As discussed,the EC technique is difficult to apply for highly heterogeneous surfaces,whereas the scintillation method using the LAS system can provide surfacelayer turbulent fluxes integrated over the LAS-path distance, suffering less due to the surface heterogeneity than the EC system.In the present study,the aim was to investigate the ability of the LAS technique to estimate surface-layer turbulent fluxes over an urban area with a significant mix of different land-use types.First,surface-layer turbulent sensible heat and momentum fluxes in a highly heterogeneous urban area were estimated using the LAS system.Then,for verification of the LAS-derived turbulent fluxes,the results were compared to those of a land-surface process-resolving numerical model and radiative flux measurements.Three-dimensional LAS footprint modeling was also performed to determine the source areas of the estimated turbulent fluxes over the heterogeneous urban area.

    The remainder of the paper is structured as follows:Section 2 briefly describes the scintillation methodology for estimating surface-layer turbulent fluxes and the 3D LAS footprint modeling.The experimental site,LAS instrumentation, and land-surface process-resolving modeling are described in section 3.Section 4 presents the LAS-derived turbulent fluxes and footprints over the heterogeneous urban area,and discusses the uncertainty of the LAS-derived turbulent fluxes through a series of sensitivity tests.A summary and conclusions follow in section 5.

    2.Methodology

    2.1.Scintillation method

    When electromagnetic radiation propagates through the atmosphere,its intensity is affected by the atmospheric turbulence,leading to fluctuations in the refractive index of the air (e.g.,Frehlich and Ochs,1990).The path-integrated structure parameter of the refractive index〈C2n〉can be given(Wang et al.,1978)by

    where r(=x/LLAS)is the normalized distance along the LAS path(LLAS)between transmitter and receiver,and C2n(r)is the structure parameter at distance r.Here,x and LLASindicate distance from the transmitter and optical path length, respectively.W(r)is the LAS weighting function,which has a mathematical bell-shape curve with a maximum value in the middle of the optical path(e.g.,Chehbouni et al.,2000; Meijninger and De Bruin,2000;Timmermans et al.,2009; Geli et al.,2012):

    where K(=2π/λ)is the optical wavenumber,k is the turbulent spatial wavenumber,?n(k)(=0.033k-11/3)is the 3D Kolmogorov spectrum of the refractive index,J1is a Bessel function of the first kind for order one,with x1=K Dr/2 and x2=[K D(1-r)]/2,where D is the aperture diameter of the LAS transmitter/receiver.

    whereσ2lnIis the statistical variance of the logarithm of intensity fluctuation obtained from measured radiation intensity. Hereafter,〈C2n〉is denoted as C2nfor simplicity.

    Because the intensity fluctuation changes with air temperature and humidity,C2ncan be expressed as a combination of the structure parameters of temperature(C2T),humidity(C2Q), and their covariance fluctuations(C2TQ).The influence by atmospheric pressure fluctuation can be negligible(Hill et al., 1980).The scintillation is predominantly influenced by temperature fluctuation forvisible and infrared radiation(Wesely, 1976;Moene,2003),leading to a simplified relation between C2n(m-2/3)and C2T(K2m-2/3)as follows:

    where T is the air temperature(K),P is the atmosphericpressure(Pa),andβis the Bowen ratio,which represents a humidity correction in scintillation.

    where zeffis the LAS effective height(m),d is the zero-plane displacement height(m),and fTis a dimensionless universal function for the temperature structure function parameter. The L is the Obukhov length,defined by

    where k(=0.4)is the von Karman constant and g(=9.8 m s-2)is the gravitational acceleration.The non-dimensional universal function for the temperature structure function parameter can be formulated for unstable atmospheric conditions following Andreas(1988):

    The friction velocity(u?)is also estimated by the following similarity relation(e.g.,Panofsky and Dutton,1984;Foken, 2008):

    where U is the mean wind speed at zeff,andψmis the integrated universal function for momentum flux(e.g.,Paulson, 1970;Businger et al.,1971).The integrated universal function used in this study for the unstable atmospheric surface layer is(H¨ogstr¨om,1988)

    Given the similarity relations,the temperature scale(T?), friction velocity(u?),and the Obukhov length(L)can be solved iteratively.The turbulent sensible heat flux(QH; W m-2)and momentum flux(τ;kg m-1s-2)can be calculated as follows:

    whereρis the air density(kg m-3)and cpis the specific heat of the air at a constant pressure(J kg-1K-1).

    Additional parameters required for estimation of the turbulent fluxes are the LAS effective height(zeff),surface aerodynamic roughness length(z0),and zero-plane displacement height(d).First,the LAS effective height can be estimated iteratively from the following relation(Hartogensis et al., 2003):

    where z(r)is the optical beam height at distance r.This relation indicates that zeffis dependent on atmospheric stability, and thus it was estimated at every time increment(30 min)of the available data in this study.

    2.2.Three-dimensional LAS footprint modeling

    Many flux footprint models have been suggested based on analytical and numerical approaches and applied in the interpretation of EC flux measurements(e.g.,Schuepp et al., 1990;Horst and Weil,1992,1994;Schmid,1994;Horst, 1999;Hsieh et al.,2000;Kormann and Meixner,2001;Kljun et al.,2004).The 3D LAS footprint can be calculated as an extension of the flux footprint model;that is,a combination of the flux footprint calculated at each segment location along the optical path and the LAS spatial weighting function(e.g., Hoedjes et al.,2002;Meijninger et al.,2002;Hartogensis et al.,2003).The turbulent flux footprint relation used for the LAS footprint modeling in this study is briefly introduced below.

    The turbulent flux F(x,y,zm)measured at height zmcan be related to the footprint function f(x,y,zm)and the spatial distribution of the surface fluxes F0(x′,y′,z′=0)as follows (Horst and Weil,1992):

    Here,the flux footprint function(units:m-2)is determined by the product of the crosswind distribution function Dy(x,y) (units:m-1)and the crosswind-integrated footprint function fy(x,zm)(units:m-1)(van Ulden,1978;Horst and Weil, 1992)as

    where x is the downwind distance from the measurement location,and y is the crosswind distance from the centerline. The dispersion to the crosswind direction is assumed to be symmetric using a Gaussian distribution function(Pasquill, 1974):

    whereσyis the standard deviation of the plume in the y dimension,which relies on atmospheric turbulence intensity, upwind distance(x),and mean wind speed[u(x)]asσvx/u(x) for a short-range dispersion(Schmid,1994).Here,the standard deviation of lateral wind fluctuations,σv(=1.9u?),was calculated following Panofsky and Dutton(1984).On the other hand,the footprint model by Hsieh et al.(2000)was used to calculate the crosswind-integrated footprint:

    where D and P are specific constants that depend on atmospheric stability,and zuis a length scale defined as

    The constant parameters are D=0.28 and P=0.59 for unstable conditions(Hsieh et al.,2000).

    3.Measurements and land-surface processresolving numerical modeling

    3.1.Site description

    The field experiment was conducted in a small urban area in Gongju,South Korea,where there is a population of 117 000 people.The city is located in a basin surrounded by extensive agricultural and hilly regions.It is also located inland, about 50 km west of the Yellow Sea.Figure 1a shows the experimental site within the city.The urbanized area is largely composed of single(1–3 story)residential houses,and a river (named“Geum River”)crosses over the city.Relatively new buildings are located alongside the northern stretches of the river,showing the directions of the city’s sprawl.The vegetation fraction within the built-up region ranges from 5% to 10%,and hilly mountains surround the city in the south and east directions.The measurement site was highly heterogeneous in terms of its land cover and topography,which is a typical morphological configuration that can be found in many cities of South Korea.

    3.2.LAS instrumentation

    The LAS instruments were deployed at two buildings of the Kongju National University,located within the city,to investigate the ability in estimating turbulent sensible heat and momentum fluxes over the heterogeneous urban area(Fig. 1a).The receiver(marked A in Fig.1a)and transmitter (marked B in Fig.1a)were installed on the rooftops of buildings that are 49 m and 48 m above mean sea level,and the optical path extends 2.1 km from southwest to northeast,crossing the Geum River(A–B in Fig.1a).This configuration of the LAS instruments measured area-averaged turbulent sensible heat and momentum fluxes that were representative of the heterogeneous urban area(residential area,bare ground, vegetated area,water body,etc.).

    The LAS system used was the LAS MkII manufactured by Kipp&Zonen,with an aperture diameter of D=0.149 m.The light source of the LAS MkII transmitter operates at a near-infrared wavelength of 850 nm,at which the scintillations measured by C2nare primarily influenced by turbulent temperature fluctuations.Based on Kleissl et al.(2010), the threshold C2nvalue for signal saturation was estimated as 4.06×10-14m-2/3for our LAS experimental configuration. The LAS sampled C2nat a 1 Hz frequency,which were averaged to 15 s intervals for turbulent flux estimation.The air temperature,wind speed,and atmospheric pressure were measured at the same temporal frequency from separate me-teorological sensors,but which were connected to the LAS system at the same height of the receiver.In addition,a 7 m meteorological tower on the rooftop of an eight-story building was also operated at the location of 320 m north of the receiver along the optical path.The meteorological data measured from the meteorological tower were used for the calculation of the LAS-derived turbulent fluxes.

    3.3.Case selection and surface input parameters

    The field measurements were conducted for 29 days during 5–19 November and 5–18 December 2013,in which adverse meteorological conditions(e.g.,fog)affected the measurements for many days.In this study,two sequential clear days on 11 and 12 November 2013 were selected to investigate the capability of the LAS in estimating the turbulent fluxes over the heterogeneous surfaces.Furthermore,since the LAS cannot discriminate atmospheric stability from the measured C2nvalues,the daytime period(0800–1700 LST) data were analyzed for estimation of the LAS-derived turbulent fluxes in unstable atmospheric conditions.

    Figure 1b shows the topographical variation along the LAS beam path and the optical beam height relative to the ground.From Eq.(11),the LAS effective heights(zeff)were estimated as 25.1–27.2 m during the daytime periods of 11 and 12 November,which satisfied a“saturation condition”required for reasonable estimation of turbulent fluxes using the LAS system(Kipp&Zonen,2012).The Bowen ratio (β)was initially assigned as 1.The effective aerodynamic roughness length and displacement height were assigned as z0=0.3 m and d=2.1 m,respectively.The aerodynamic parameters were determined following the general rule of z0≈0.1h and d≈0.7h(Garratt,1992;Grimmond and Oke, 1999),in which h(~3 m)indicates the mean height of obstacles along the beam path.Generally,the determination of surface aerodynamic parameters can be obtained by applying a morphological method(e.g.,Raupach,1994;Macdonald et al.,1998)and/or surface-layer similarity relation using EC measurements taken in the atmospheric surface layer for homogeneous surfaces(e.g.,Grimmond and Oke,1999;Roth, 2000).However,accurate estimation of the surface aerodynamic parameters is a very difficult task for heterogeneous surfaces due to intrinsic limitations of the methodologies.No standard method that determines the surface aerodynamic parameters for heterogeneous surfaces exists,which inevitably leads to uncertainty in the LAS-derived surface-layer turbulent fluxes.In this study,a series of sensitivity tests were performed to quantify the influence of errors in the surface input parameters on the LAS-derived turbulent fluxes.

    3.4.Land-surface process-resolving numerical modeling

    The ability of the LAS equipment used in this study in estimating surface-layer turbulent fluxes has been validated against the EC system for homogeneous land surfaces (e.g.,Wilson et al.,2013).Due to the difference of the footprints of the EC and LAS-derived turbulent fluxes,direct comparison of the two turbulent fluxes is difficult for highly heterogeneous surfaces.In order to verify the LAS-derived turbulent fluxes over the heterogeneous field experimental area,land-surface process-resolving modeling was used in the present study.Specifically,we conducted very high resolution urban numerical modeling using the Weather Research and Forecasting(WRF)model(Skamarock et al., 2008).Given a proper experimental configuration,the WRF model has the capability to explicitly represent significant surface land-use heterogeneity over the field experiment region and to resolve associated land-surface physical(radiative/dynamic/hydrological)processes.

    The WRF model was configured with five nested domains,covering East Asia with a 16.2 km grid spacing in domain 1 and nesting down to the Gongju region with a 200 m grid spacing in domain 5(Fig.2).To realistically represent the significant heterogeneity of land-use and topography over the experimental region,very high resolution static datasets of the Moderate Resolution Imaging Spectroradiometer (MODIS)-based land-use(Kang et al.,2010;spatial resolution of 7.5′′)and Shuttle Radar Topography Mission(SRTM) topography(http://www2.jpl.nasa.gov/srtm/index.html;spatial resolution of 3′′)were implemented as input data of the WRF model.Physical processes applied in the simulation included:Goddard shortwave radiation(Chou and Suarez, 1994),Rapid Radiative Transfer Model(RRTM)longwave radiation(Mlawer et al.,1997),Lin cloud microphysics(Lin et al.,1983),and a turbulent kinetic energy(TKE)–based turbulence scheme(Janji′c,2002).Land-surface processes were parameterized following the Noah land surface model(Chen and Dudhia,2001)with a single-layer urban canopy model (Kusaka et al.,2001)for urban patches.

    The WRF simulation was conducted starting from 0000 UTC 9 November 2013 to 0000 UTC 13 November 2013. In order to realistically reproduce surface-layer meteorology and relevant turbulent fluxes,observed wind fields from the meteorological site were nudged during the simulation period.The WRF model calculates surface-layer turbulent sensible heat and momentum fluxes at every grid mesh.The WRF-simulated turbulent fluxes with a 30 min time interval were extracted along the LAS beam path,which were multiplied by the LAS weighting function for comparison with the LAS-derived fluxes.In this study,a WRF model sensitivity simulation was also conducted,by replacing the Geum River with a natural surface type(“shrubland”)to further interpret the LAS-derived turbulent fluxes(discussed in section 4.2).

    4.Results

    4.1.Meteorological conditions and LAS footprints

    The weather conditions on 11 and 12 November 2013 were fairly good,with little influence by fog and clouds.As shown in Fig.3,the prevailing winds were northeasterly or northwesterly during the daytime on both days,except for the morning hours on 12 November.The wind speed changed within the range of0.5–4 m s-1during the daytime,being relatively low in the morning hours(0.5–2 m s-1)and relatively high in the afternoon(2–4 m s-1),while the air temperatureshowed a typical temporal variation for both the days,ranging from 0?C to 9?C.

    Based on the meteorological conditions(Table 1),3D LAS footprint modeling was conducted to interpret the LAS-derived turbulentfluxes.In this study,the flux footprints were calculated every 10 m in discrete locations along the LAS optical path,and then the estimates were multiplied by the LAS weighting function to model the 3D LAS footprints.As shown in Fig.4,the calculated LAS footprints clearly revealed the source regions of the LAS-derived turbulent fluxes over the heterogeneous urban area.The source regions of the LAS-derived fluxes depended highly on the angle between the wind direction and the LAS beam path,as well as the meteorological conditions(e.g.,atmospheric stability,wind speed,friction velocity).When the wind direction was almost parallel to the LAS optical path(Figs.4a and c),the calculated footprint indicated that 90%of the daytime turbulent fluxes were originated from areas within a width of 100–300 m along the optical path(0.2–0.6 km2).It is also clear that the width of in fluence became narrower as the atmosphericstability decreased(more unstable).On the other hand,the daytime turbulent fluxes of 90%were attributed to larger areas,expanding up to about 1 km(~2 km2)upwind from the LAS optical path,when the wind direction was almost perpendicular to the optical path(Figs.4b and d).The results also show thatthe LAS-derived turbulentfluxeswere the least influenced by sources near the LAS transmitter/receiver,as discussed by Chehbouni et al.(2000).

    Table 1.Time-averaged meteorological variables used for the representative LAS(large aperture scintillometer)footprint modeling.The wind direction and speed were obtained from the meteorologicaltower,and the Obukhov length(L)and friction velocity(u?)were estimated using the scintillation method.

    By virtue of the 3D footprint analysis,it was found that the LAS system,operating with an optical path length of 2.1 km and effective height of about 26 m,can represent the areaaveraged turbulent fluxes over this heterogeneous urban area on a micro-αscale,i.e.,0.2–2 km.The results imply that the LAS’s spatial-integration capability will be very useful to evaluate modeled surface-layer turbulent fluxes over significantly heterogeneous surfaces on the spatial scale.

    4.2.LAS-derived sensible heat and momentum fluxes

    The measured Cn2ranged from 0.1 to 3.5×10-14m-2/3and 0.1 to 2.8×10-14m-2/3during the daytime periods (0800–1700 LST)on 11 and 12 November,respectively,thus satisfying the saturation-free condition described in section 3.2.The LAS-derived turbulent sensible heat and momentum fluxes were calculated every 15 s from the Cn2values following the procedure described in section 2,and then the fluxes were averaged over 30 min intervals.Figure 5 shows that the LAS-derived momentum fluxes(τ)on 11 and 12 November 2013,which were well correlated with the mean wind speed profiles shown in Fig.3,ranged from 0.05 to 0.35 N m-2during the daytime.The land-surface process-resolving WRF simulation also compared well to the estimation,indicating the LAS estimation is able to represent the spatial-integrated momentum fluxes over heterogeneous surfaces.The ratio of friction velocity to mean wind speed(u?/U)varied from 0.14 in the early morning to 0.11 at later times on 11 November, and from 0.22 during morning to 0.10 in the afternoon on 12 November.The u?/U difference between the morning and the afternoon is attributed to the temporal variation of atmospheric stability,following Eq.(8).

    Figure 6 shows that the LAS-derived daytime sensible heat flux(QH)varied between 25 and 160 W m-2on 11 November,and between 20 and 120 W m-2on 12 November.The LAS-derived QHquantities were different in magnitude and temporal variation between the two days,with a dramatic reduction at around 1400 LST 12 November.Temporal fluctuations in the LAS-derived turbulent fluxes re flected the temporal variability of the measured scintillation intensity for the daytime periods.The calculated temporal variability(1σ) in QHwas estimated to be as much as about 21%relative to the 30 min averaged values on 11 November,and 22%on 12 November.The temporal variability inτand its footprint dependency were not as large as in QH(Fig.5)because the momentum fluxes were not directly associated with the scintillation intensity.

    The LAS-derived heat fluxes were in good agreement with the land-surface process-resolving simulation results, especially during morning hours(0800–1200 LST).Further comparison to the net radiative fluxes measured at the meteorological site shows that the estimated daily mean QH(107 and 72 W m-2)corresponded to about 49%and 37%of the net radiation of 219 and 194 W m-2on 11 and 12 November 2013,respectively.Despite the intrinsic difficulty in the quantitative verification of the LAS-derived turbulent fluxes over heterogeneous surfaces,the comparison results of the LAS-derived QHto the WRF simulation and the flux ratio analysis show that the scintillation method can reasonably estimate area-averaged surface-layer turbulent fluxes.

    The distinctive features of the LAS-derived QHtemporal variation are that the peak QH(especially on 11 November)occurred earlier in the day and that morning QHwas enhanced compared to that in the afternoon(asymmetric distribution).Generally,a peak QHvalue in the surface-layer en-ergy balance is found at 1300–1500 LST on clear days(e.g., Lee and Park,2008;Lee,2011;Loridan et al.,2013).In addition,the maximum values tend to lag in time over urban areas compared to natural surfaces(Lee and Baik,2011).To further investigate the features in the LAS-derived QH,another WRF simulation was conducted by replacing the Geum River with a natural surface type(“shrubland”).The WRF simulation revealed that the peak time lag and the morning enhancement in the LAS-derived QHcan be partly attributed to the existence of the river,where signi ficantly larger QHvalues(~200 W m-2)were found in morning hours(0800–1000 LST)due to the large temperature difference between the water body(~14.0?C)and the overlying atmosphere in the cold season.

    4.3.Sensitivity analysis

    The LAS system produced surface-layer turbulent fluxes over this heterogeneous urban area within a reasonable uncertainty range(1σ).In this section,the results from several sensitivity tests are discussed to further quantify the influence of input parameters(surface aerodynamic parameters and Bowen ratio)and surface-layer similarity function for the temperature structure function on the LAS-derived turbulent fluxes.

    The surface aerodynamic parameters(z0and d)could be estimated by an average of the surface parameters of each homogeneous patch consisting of heterogeneous source areas, as described in section 2.1.However,in practice it is dif ficult to determine the surface parameters for highly heterogeneous surfaces because of the complexity of the morphology of buildings/obstacles and the topography.The uncertainty in the surface aerodynamic parameters also influences the determination of the LAS effective height(zeff).As a realistic ballpark guess of the surface parameters for the study area (e.g.,Garratt,1992;Foken,2008),an aerodynamic roughness length of 0.1 m and 1.0 m and an LAS effective height of about±20%(21 and 31 m)from the base estimation were applied in the sensitivity tests.Besides the surface aerodynamic parameters,the influence of the Bowen ratio was also considered in the sensitivity tests,with values ofβ=0.5 and β=2.0.

    Figure 7 shows that the estimatedτwas very sensitive to z0,leading to a reduction ofabout-34%for z0=0.3 m butan increase of about+80%for z0=1.0 m,compared to the base estimation on both days.The influence of the Bowen ratio on the LAS-derivedτwas negligible,while the change in zeffled to momentum-flux difference within±5%.On the other hand,as shown in Fig.8,the estimated QHwas more variable to the change in zeff,ranging between±17%on a daytime average.The change in z0led to the flux difference of between -6%and 14%on 11 November,which was lower than that of between-4%and 9%on 12 November.The difference inβled to the smallest(<4%)changes in QH.Overall,the sensitivity tests indicate that the uncertainty in the surface parameters may affect the LAS-derived QHwithin±20%.

    To quantify the uncertainty by the choice of the temperature structure function parameter,the LAS turbulent fluxes were estimated using four different non-dimensional similarity functions for the temperature structure function parameter.When calculating turbulent sensible heat flux using thescintillation method,previous studies have used various nondimensional similarity functions(e.g.,Meijninger and De Bruin,2000;Lagouarde et al.,2006;Evans et al.,2012;Geli et al.,2012).Along with the surface-layer structure function suggested by Andreas(1988),the function by De Bruin et al. (1993)has been popularly used,which is given for unstable conditions(z/L<0)as follows:

    Foken and Kretschmer(1990)and Thiermann and Grassl (1992)proposed different functional forms as follows:

    where k=0.4,and

    respectively.By comparing the non-dimensional similarity functions for unstable atmospheric conditions,it was found that the function by Andreas(1988)had higher fTthan that by De Bruin et al.(1993),except for near neutral conditions (-z/L<0.01),while it had relatively lower fTthan the functions proposed by Foken and Kretschmer(1990)and Thiermann and Grassl(1992)for a wide range of unstable atmospheric conditions(Fig.9).

    As shown in Fig.10,the QHfrom the Andreas function was lower by~14%(range:8%–17%)in daytime average than that by the De Bruin function,but larger by about 9%and 10%than those from the functions of Foken and Kretschmer(1990)and Thiermann and Grassl(1992),respectively.The estimated QHquantities were well correlated in time with each other,with a correlation coefficient of 0.99.

    The influence of the temperature similarity functions onτ were less than 2%in daytime mean flux(not shown).It is worth noting that Lagouarde et al.(2006)reported the EC-derived QHroughly ranged between the LAS-derived QHusing the De Bruin and Andreas functions for a reasonably homogeneous urban area.Overall,the results indicate that the uncertainty caused by the non-dimensional similarity function forthe temperature structure function parameteriswithin ±15%in the LAS-derived QH.

    5.Summary and conclusions

    This study tested the ability of the LAS technique in determining surface-layer turbulent sensible heat and momentum fluxes over a heterogeneous urban area,in comparison with land-surface process-resolving numerical modeling. The LAS system was deployed and operated over a highly heterogeneous urban area(Gongju,South Korea)during the cold season with an optical path length of 2.1 km.Threedimensional LAS footprint modeling was also performed to interpret the LAS-derived turbulent fluxes over the heterogeneous surfaces.The WRF model,which was configured to explicitly represent the land-use heterogeneity over the experimental site and associated physical processes,was introduced for verification of the LAS-derived turbulent fluxes.

    The LAS-derived turbulentfluxes(QHandτ)showed reasonable temporal variation during the daytime on the two clear days of 11 and 12 November 2013 when compared to the results of the land-surface process-resolving modeling and radiative flux measurements.A series of sensitivity tests showed that the overall uncertainty in the daytimeQHestimate had 20%–30%as an upper limit:within±20% in surface aerodynamic parameters and±15%in the nondimensional similarity function for the temperature structure function parameter.In addition,the estimated error inτwas very sensitive to the determination of the aerodynamic roughness length.It is noted that only two days of data were used in the analysis,which may also have contributed uncertainties to some extent.On the other hand,the footprint modeling analysis showed that the LAS-derived surface-layer turbulent fluxes estimated with the 2.1 km optical path distance were contributed by source areas within 0.2–2 km2(“microαscale”)depending on local meteorological conditions.Unlike the flux footprint modeling of the EC system,the LAS footprint(or source area)was significantly influenced by an incident angle of wind to the LAS beam path,as well as by atmospheric conditions such as atmospheric stability and wind speed.

    Recently,mesoscale meteorological and environmental models have undergone rapid development by the inclusion of urban canopy models(e.g.,Masson,2000;Kusaka et al., 2001;Martilli et al.,2002;Lee and Park,2008;Oleson et al., 2008;Porson et al.,2010;Lee,2011;Ryu et al.,2011)for better representation of urban land-surface processes(e.g., Lee et al.,2011).Along with numerical model development,quantitative validation of surface-layer turbulent fluxes modeled by 3D mesoscale meteorological models is very important to improving forecast skill,as well as understanding boundary layer structures.Even though previous studies have successfully determined the surface-layer turbulent fluxes over various homogeneous surfaces using the EC system and/orLAS system(e.g.,Lagouarde etal.,2006;Zeweldi et al.,2010;Evans et al.,2012;Ward et al.,2014),few studies have provided estimations over heterogeneous urban surfaces.The present study demonstrates that the LAS technique can provide realistic surface-layer turbulent fluxes over highly heterogeneous urban areas within a reasonable uncertainty range.In addition,it shows that LAS footprint modeling can provide valuable insight for model–measurement comparisons by the interpretation of sources’contributions to the LAS estimates.Furthermore,this study indicates that LAS-derived turbulent fluxes can be useful for the verification of mesoscale meteorological and environmental modeling conducted over heterogeneous surfaces.

    Acknowledgements.The authors gratefully acknowledge the two anonymous reviewers and the editor for their invaluable comments.We also thank Doo-Il LEE for carrying out the WRF simulation.This work was supported by the Korea Meteorological Administration Research and Development Program(Grant No.CATER 2012-3081).

    REFERENCES

    Andreas,E.L.,1988:Atmospheric stability from scintillation measurements.Appl.Opt.,27,2241–2246.

    Beyrich,F.,H.A.R.De Bruin,W.M.L.Meijninger,J.W.Schipper,and H.Lohse,2002:Results from one-year continuous operation of a large aperture scintillometer over a heterogeneous land surface.Bound.-Layer Meteor.,105,85–97.

    Businger,J.A.,J.C.Wyngaard,Y.Izumi,and E.F.Bradley,1971: Flux-profile relationships in the atmospheric surface layer.J. Atmos.Sci.,28,181–189.

    Chehbouni,A.,and Coauthors,2000:Estimation of heat and momentum fluxes over complex terrain using a large aperture scintillometer.Agricultural and Forest Meteorology,105, 215–226.

    Chen,F.,and J.Dudhia,2001:Coupling an advanced land surfacehydrology model with the Penn State-NCAR MM5 modeling system.Part I:Model implementation and sensitivity.Mon. Wea.Rev.,129,569–585.

    Chou,M.D.,and M.J.Suarez,1994:An Efficient Thermal Infrared Radiation Parameterization for Use in General Circulation Models.NASA Tech Memo,104606,Vol.3,85 pp.

    De Bruin,H.A.R.,W.Kohsiek,and B.J.J.M.van den Hurk, 1993:A verification of some methods to determine the fluxes of momentum,sensible heat,and water vapour using standard deviation and structure parameter of scalar meteorological quantities.Bound.-Layer Meteor.,63,231–257.

    De Bruin,H.A.R.,B.J.J.M.van den Hurk,and W.Kohsiek, 1995:The scintillation method tested over a dry vineyard area.Bound.-Layer Meteor.,76,25–40.

    De Bruin,H.A.R.,W.M.L.Meijninger,A.S.Smedman,and M.Magnusson,2002:Displaced-beam small aperture scintillometer test.Part I:The WINTEX data-set.Bound.-Layer Meteor.,105,129–148.

    Evans,J.G.,D.D.McNeil,J.W.Finch,T.Murray,R.J.Harding,H.C.Ward,and A.Verhoef,2012:Determination of turbulent heat fluxes using a large aperture scintillometer over undulating mixed agricultural terrain.Agricultural and Forest Meteorology,166–167,221–233.

    Foken,T.,2008:Micrometeorology.Springer,306 pp.

    Foken,T.,and D.Kretschmer,1990:Stability dependence of the temperature structure parameter.Bound.-Layer Meteor.,53, 185–189.

    Frehlich,R.G.,and G.R.Ochs,1990:Effects of saturation on the optical scintillometer.Appl.Opt.,29,548–553.

    Garratt,J.R.,1992:The Atmospheric Boundary Layer.Cambridge University Press,316 pp.

    Geli,H.M.E.,C.M.U.Neale,D.Watts,J.Osterberg,H.A. R.De Bruin,W.Kohsiek,R.T.Pack,and L.E.Hipps, 2012:Scintillometer-based estimates of sensible heat flux using lidar-derived surface roughness.Journal of Hydrometeorology,13,1317–1331.

    Grimmond,C.S.B.,and T.R.Oke,1999:Aerodynamic properties of urban areas derived from analysis of surface form.J.Appl. Meteor.,38,1262–1292.

    Hartogensis,O.K.,C.J.Watts,J.C.Rodriguez,and H.A.R. De Bruin,2003:Derivation of an effective height for scintillometers:La Poza experiment in northwest Mexico.Journal of Hydrometeorology,4,915–928.

    Hill,R.J.,S.F.Clifford,and R.S.Lawrence,1980:Refractiveindex and absorption fluctuations in the infrared caused by temperature,humidity,and pressure fluctuations.Journal of the Optical Society of America,70,1192–1205.

    Hoedjes,J.C.B.,R.M.Zuurbier,and C.J.Watts,2002:Large aperture scintillometer used over a homogeneous irrigated area,partly affected by regional advection.Bound.-Layer Meteor.,105,99–117.

    H¨ogstr¨om,U.,1988:Non-dimensional wind and temperatureprofiles in the atmospheric surface layer:A re-evaluation. Bound.-Layer Meteor.,42,55–78.

    Horst,T.W.,1999:The footprint for estimation of atmospheresurface exchange fluxes by profile techniques.Bound.-Layer Meteor.,90,171–188.

    Horst,T.W.,and J.C.Weil,1992:Footprint estimation for scalar flux measurements in the atmospheric surface layer.Bound.-Layer Meteor.,59,279–296.

    Horst,T.W.,and J.C.Weil,1994:How far is far enough?The fetch requirements for micrometeorological measurement of surface fluxes.J.Atmos.Oceanic Technol.,11,1018–1025.

    Hsieh,C.I.,G.Katul,and T.W.Chi,2000:An approximate analytical model for footprint estimation of scalar fluxes in thermally stratified atmospheric flows.Advances in Water Resources,23,765–772.

    Janji′c,Z.I.,2002:Nonsingular implementation of the Mellor-Yamada Level 2.5 Scheme in the NCEP Meso Model.NCEP Office Note,No.437,61 pp.

    Kanda,M.,R.Moriwaki,M.Roth,and T.Oke,2002:Areaaveraged sensible heat flux and a new method to determine zero-plane displacement length over an urban surface using scintillometry.Bound.-Layer Meteor.,105,177–193.

    Kang,J.-H.,M.-S.Suh,and J.-H.Kwak,2010:Land cover classification over East Asian region using recent MODIS NDVI data(2006–2008).Atmosphere,20,415–426.(Korean with an English abstract)

    Kipp&Zonen,2012:Instruction Manual.Delft,,Netherlands,86 pp.

    Kleissl,J.,O.K.Hartogensis,and J.D.Gomez,2010:Test of scintillometer saturation correction methods using field experimental data.Bound.-Layer Meteor.,137,493–507.

    Kljun,N.,P.Calanca,M.W.Rotach,and H.P.Schmid,2004: A simple parameterisation for flux footprint predictions. Bound.-Layer Meteor.,112,503–523.

    Kormann,R.,and F.X.Meixner,2001:An analytical footprint model for non-neutral stratification.Bound.-Layer Meteor., 99,207–224.

    Kusaka,H.,H.Kondo,Y.Kikegawa,and F.Kimura,2001:A simple single-layer urban canopy model for atmospheric models: Comparison with multi-layer and slab models.Bound.-Layer Meteor.,101,329–358.

    Lagouarde,J.P.,M.Irvine,J.M.Bonnefond,C.S.B.Grimmond,N.Long,T.R.Oke,J.A.Salmond,and B.Offerle, 2006:Monitoring the sensible heat flux over urban areas using large aperture scintillometry:Case study of Marseille city during the ESCOMPTE experiment.Bound.-Layer Meteor., 118,449–476.

    Lee,S.-H.,2011:Further development of the vegetated urban canopy model including a grass-covered surface parametrization and photosynthesis effects.Bound.-Layer Meteor.,140, 315–342.

    Lee,S.-H.,and S.-U.Park,2008:A vegetated urban canopy model for meteorological and environmental modelling.Bound.-Layer Meteor.,126,73–102.

    Lee,S.-H.,and J.-J.Baik,2011:Evaluation of the vegetated urban canopy model(VUCM)and its impacts on urban boundary layer simulation.Asia-Pacific Journal of Atmospheric Sciences,47(2),151–165.

    Lee,S.-H.,and Coauthors,2011:Evaluation of urban surface parameterizations in the WRF model using measurements during the Texas Air Quality Study 2006 field campaign.Atmos. Chem.Phys.,11,2127–2143.

    Lin,Y.L.,R.D.Farley,and H.D.Orville,1983:Bulk parameterization of the snow field in a cloud model.J.Climate Appl. Meteor.,22,1065–1092.

    Liu,S.M.,Z.W.Xu,W.Z.Wang,Z.Z.Jia,M.J.Zhu,and J. M.Wang,2011:A comparison of eddy-covariance and large aperture scintillometer measurements with respect to the energy balance closure problem.Hydrology and Earth System Sciences,15,1291–1306.

    Liu,S.M.,Z.W.Xu,Z.L.Zhu,Z.Z.Jia,and M.J.Zhu,2013: Measurements of evapotranspiration from eddy-covariance systems and large aperture scintillometers in the Hai River Basin,China.J.Hydrol.,487,24–38.

    Loridan,T.,F.Lindberg,O.Jorba,S.Kotthaus,S.Grossman-Clarke,and C.S.B.Grimmond,2013:High resolution simulation of the variability of surface energy balance fluxes across central London with urban zones for energy partitioning.Bound.-Layer Meteor.,147,493–523.

    Martilli,A.,A.Clappier,and M.W.Rotach,2002:An urban surface exchange parameterisation for mesoscale models. Bound.-Layer Meteor.,104,261–304.

    Macdonald,R.W.,R.F.Griffiths,and D.J.Hall,1998:An improved method for the estimation of surface roughness of obstacle arrays.Atmos.Environ.,32,1857–1864.

    Masson,V.,2000:A physically-based scheme for the urban energy budget in atmospheric models.Bound.-Layer Meteor., 94,357–397.

    McAneney,K.J.,A.E.Green,and M.S.Astill,1995:Largeaperture scintillometry:The homogeneous case.Agricultural and Forest Meteorology,76,149–162.

    Meijninger,W.M.L.,and H.A.R.De Bruin,2000:The sensible heat fluxes over irrigated areas in western Turkey determined with a large aperture scintillometer.J.Hydrol.,229,42–49.

    Meijninger,W.M.L.,A.E.Green,O.K.Hartogensis,W.Kohsiek, J.C.B.Hoedjes,R.M.Zuurbier,and H.A.R.De Bruin, 2002:Determination of area-averaged water vapour fluxes with large aperture and radio wave scintillometers over a heterogeneous surface-Flevoland field experiment.Bound.-Layer Meteor.,105,63–83.

    Mlawer,E.J.,S.J.Taubman,P.D.Brown,M.J.Iacono,and S.A. Clough,1997:Radiative transfer for inhomogeneous atmospheres:RRTM,a validated correlated-k model for the longwave.J.Geophys.Res.,102,16 663–16 682.

    Moene,A.F.,2003:Effects of water vapour on the structure parameter of the refractive index for near-infrared radiation. Bound.-Layer Meteor.,107,635–653.

    Offerle,B.,C.S.B.Grimmond,K.Fortuniak,and W.Pawlak, 2006:Intraurban differencesof surface energy fluxesin a central European city.J.Appl.Meteor.Climatol.,45,125–136.

    Oleson,K.W.,G.B.Bonan,J.Feddema,M.Vertenstein,and C.S.B.Grimmond,2008:An urban parameterization for a global climate model.Part I:Formulation and evaluation for two cities.J.Appl.Meteor.Climatol.,47,1038–1060.

    Panofsky,H.A.,and J.A.Dutton,1984:Atmospheric Turbulence: Models and Methods for Engineering Applications.John Wiley and Sons,New York,397 pp.

    Pasquill,F.,1974:Atmospheric Diffusion.2nd ed.,John Wiley& Sons,New York,425 pp.

    Paulson,C.A.,1970:The mathematical representation of wind speed and temperature profiles in the unstable atmospheric surface layer.J.Appl.Meteor.,9,857–861.

    Pauscher,L.,2010:Scintillometer measurements above the urban area of London.Diploma thesis,Dept.of Micrometeorology,University of Bayreuth,95 pp.

    Porson,A.,P.A.Clark,I.N.Harman,M.J.Best,and S.E.Belcher, 2010:Implementation of a new urban energy budget scheme in the MetUM.Part I:Description and idealized simulations. Quart.J.Roy.Meteor.Soc.,136,1514–1529.

    Raupach,M.R.,1994:Simplified expressions for vegetation roughness length and zero-plane displacement as functions of canopy height and area index.Bound.-Layer Meteor.,71, 211–216.

    Roth,M.W.,2000:Review of atmospheric turbulence over cities. Quart.J.Roy.Meteor.Soc.,126,941–990.

    Ryu,Y.-H.,J.-J.Baik,and S.-H.Lee,2011:A new single-layer urban canopy model for use in mesoscale atmospheric models. J.Appl.Meteor.Climatol.,50,1773–1794.

    Schmid,H.P.,1994:Source areas for scalars and scalar fluxes. Bound.-Layer Meteor.,67,293–318.

    Schuepp,P.H.,M.Y.Leclerc,J.I.Macpherson,and R.L.Desjardins,1990:Footprint prediction of scalar fluxes from analytical solutions of the diffusion equation.Bound.-Layer Meteor.,50,355–373.

    Skamarock,W.C.,and Coauthors,2008:A description of the advanced research WRF version 3.NCAR Technical Note, NCAR/TN-475+STR,113 pp.

    Tatarskii,V.I.,1961:Wave Propagation in a Turbulent Medium. McGraw-Hill,285 pp.

    Thiermann,V.,and H.Grassl,1992:The measurement of turbulent surface-layer fluxes by use of bichromatic scintillation. Bound.-Layer Meteor.,58,367–389.

    Timmermans,W.J.,Z.Su,and A.Olioso,2009:Footprintissues in scintillometry over heterogeneous landscapes.Hydrol.Earth Syst.Sci.,13(11),2179–2190.

    van Ulden,A.P.,1978:Simple estimates for vertical diffusion from sources near the ground.Atmos.Environ.,12,2125–2129.

    Wang,T.I.,G.R.Ochs,and S.F.Clifford,1978:A saturationresistant optical scintillometer to measure C2n.Journal of the Optical Society of America,69,334–338.

    Ward,H.C.,J.G.Evans,and C.S.B.Grimmond,2014:Multiscale sensible heat fluxes in the suburban environment from large-aperture scintillometry and eddy covariance.Bound.-Layer Meteor.,152,65–89.

    Wesely,M.L.,1976:The combined effect of temperature and humidity fluctuations on refractive index.J.Appl.Meteor.,15, 43–49.

    Wilson,K.M.,A.van Tol,and J.Mes,2013:The upgraded Kipp&Zonen LAS MkII large aperture scintillometer instrument specifications.Tubingen Atmospheric Physics Symposium“Scintillometers and Applications”,7–9 Oct 2013,Tubingen,Germany.

    Wyngaard,J.C.,Y.Izumi,and S.A.Collins Jr.,1971:Behavior of the refractive-index-structure parameter near the ground. Journal of the Optical Society of America,61,1646–1650.

    Xu,Z.W.,and Coauthors,2013:Intercomparison of surface energy flux measurement systems used during the HiWATERMUSOEXE.J.Geophys.Res.,118,13 140–13 157.

    Zeweldi,D.A.,M.Gebremichael,J.M.Wang,T.Sammis,J. Kleissl,and D.Miller,2010:Intercomparison of sensible heat flux from large aperture scintillometer and eddy covariance methods:Field experiment over a homogeneous semi-arid region.Bound.-Layer Meteor.,135,151–159.

    :Lee,S.-H.,J.-H.Lee,and B.-Y.Kim,2015:Estimation of turbulent sensible heat and momentum fluxes over a heterogeneous urban area using a large aperture scintillometer.Adv.Atmos.Sci.,32(8),1092–1105,

    10.1007/s00376-015-4236-2.

    25 October 2014;revised 16 December 2014;accepted 9 January 2015)

    ?Corresponding author:Sang-Hyun LEE

    Email:sanghyun@kongju.ac.kr

    又大又爽又粗| tube8黄色片| 国产精品 国内视频| 老司机影院毛片| 精品乱码久久久久久99久播| 丝袜喷水一区| 飞空精品影院首页| 菩萨蛮人人尽说江南好唐韦庄| 日韩大片免费观看网站| 日韩一卡2卡3卡4卡2021年| 亚洲人成电影免费在线| 亚洲av美国av| 每晚都被弄得嗷嗷叫到高潮| 老司机深夜福利视频在线观看 | 不卡av一区二区三区| 国产成人a∨麻豆精品| 亚洲精品国产av成人精品| 美女福利国产在线| 国产成人欧美在线观看 | 亚洲精品美女久久av网站| h视频一区二区三区| 久久久精品免费免费高清| 国产野战对白在线观看| tube8黄色片| 久久久国产欧美日韩av| 成人国产一区最新在线观看| 欧美激情久久久久久爽电影 | 天天影视国产精品| 男女免费视频国产| 18禁裸乳无遮挡动漫免费视频| 中文字幕高清在线视频| 丝袜喷水一区| 电影成人av| 国产精品免费大片| 国产精品久久久久成人av| 精品亚洲成a人片在线观看| 视频区欧美日本亚洲| 免费观看a级毛片全部| 久久久久久久久免费视频了| 久久久久久久国产电影| 国产精品av久久久久免费| 国产成人影院久久av| 最新的欧美精品一区二区| 亚洲美女黄色视频免费看| 国产熟女午夜一区二区三区| 国产97色在线日韩免费| 国产亚洲欧美精品永久| 真人做人爱边吃奶动态| 亚洲一区中文字幕在线| 欧美成人午夜精品| 999精品在线视频| 日韩 亚洲 欧美在线| 啦啦啦啦在线视频资源| 久久久精品区二区三区| 国产深夜福利视频在线观看| 12—13女人毛片做爰片一| 日本五十路高清| 天堂俺去俺来也www色官网| 亚洲熟女毛片儿| 国产精品久久久久久精品古装| 动漫黄色视频在线观看| 久久久久久亚洲精品国产蜜桃av| 国产精品久久久久成人av| 国产在线观看jvid| 永久免费av网站大全| 日本猛色少妇xxxxx猛交久久| 一边摸一边抽搐一进一出视频| 天天影视国产精品| 亚洲精品久久久久久婷婷小说| 亚洲精品国产色婷婷电影| 亚洲欧美精品综合一区二区三区| 中文字幕人妻丝袜一区二区| 亚洲欧美成人综合另类久久久| 热re99久久国产66热| 国产在线一区二区三区精| 久久精品人人爽人人爽视色| 久久久久国产一级毛片高清牌| 老司机亚洲免费影院| 啦啦啦免费观看视频1| 最近最新中文字幕大全免费视频| 桃花免费在线播放| 久久99热这里只频精品6学生| 女警被强在线播放| 一区二区三区精品91| 成年女人毛片免费观看观看9 | 亚洲国产成人一精品久久久| 国产伦理片在线播放av一区| 亚洲av日韩在线播放| 嫁个100分男人电影在线观看| 久久久久久亚洲精品国产蜜桃av| 啦啦啦视频在线资源免费观看| 欧美日韩av久久| h视频一区二区三区| 老司机影院毛片| 欧美亚洲 丝袜 人妻 在线| 欧美少妇被猛烈插入视频| 下体分泌物呈黄色| 亚洲欧美色中文字幕在线| 成年人免费黄色播放视频| 久久久久国内视频| 欧美国产精品一级二级三级| 老司机深夜福利视频在线观看 | 三级毛片av免费| 久久ye,这里只有精品| 99国产精品99久久久久| 午夜两性在线视频| 国产淫语在线视频| 免费观看人在逋| 两性夫妻黄色片| 亚洲精品国产一区二区精华液| 我的亚洲天堂| 亚洲欧美清纯卡通| av电影中文网址| 91精品伊人久久大香线蕉| 国产一区二区 视频在线| 天堂8中文在线网| 美国免费a级毛片| 亚洲精品第二区| 狠狠精品人妻久久久久久综合| 飞空精品影院首页| av天堂在线播放| 美女午夜性视频免费| 久久毛片免费看一区二区三区| 人人妻人人爽人人添夜夜欢视频| 日韩,欧美,国产一区二区三区| bbb黄色大片| 搡老熟女国产l中国老女人| 国产真人三级小视频在线观看| 伊人亚洲综合成人网| 久久国产精品人妻蜜桃| 日韩制服丝袜自拍偷拍| 国产精品久久久久久人妻精品电影 | 母亲3免费完整高清在线观看| av片东京热男人的天堂| 欧美人与性动交α欧美精品济南到| 两性午夜刺激爽爽歪歪视频在线观看 | 欧美黑人精品巨大| 亚洲中文av在线| 黄频高清免费视频| 1024视频免费在线观看| 精品久久久久久电影网| 高清在线国产一区| 日日爽夜夜爽网站| 精品少妇内射三级| 91九色精品人成在线观看| 每晚都被弄得嗷嗷叫到高潮| 免费在线观看影片大全网站| 午夜免费鲁丝| 亚洲欧美精品综合一区二区三区| 国产精品影院久久| 1024香蕉在线观看| 国产一区二区激情短视频 | av天堂久久9| 老司机深夜福利视频在线观看 | 在线观看免费午夜福利视频| 欧美另类亚洲清纯唯美| 亚洲av片天天在线观看| 99久久99久久久精品蜜桃| 十分钟在线观看高清视频www| 久久ye,这里只有精品| 久久人人爽人人片av| 国产不卡av网站在线观看| 国产一区二区三区av在线| 啦啦啦在线免费观看视频4| 黄色 视频免费看| 久久中文看片网| 亚洲精品中文字幕一二三四区 | 老司机在亚洲福利影院| 午夜免费鲁丝| 老司机靠b影院| 成人免费观看视频高清| 操美女的视频在线观看| 女人爽到高潮嗷嗷叫在线视频| 新久久久久国产一级毛片| 国产av精品麻豆| e午夜精品久久久久久久| 久久精品国产综合久久久| 天天躁狠狠躁夜夜躁狠狠躁| 中文字幕人妻丝袜制服| 熟女少妇亚洲综合色aaa.| 久久国产精品影院| 99精品欧美一区二区三区四区| 亚洲美女黄色视频免费看| 色视频在线一区二区三区| 日韩中文字幕欧美一区二区| 午夜福利,免费看| 青草久久国产| 91精品伊人久久大香线蕉| 纯流量卡能插随身wifi吗| 精品少妇黑人巨大在线播放| 母亲3免费完整高清在线观看| 亚洲人成电影观看| 久久久久视频综合| 制服人妻中文乱码| 老司机影院成人| 超碰97精品在线观看| 在线观看免费视频网站a站| 女性生殖器流出的白浆| 国产成人影院久久av| 日本wwww免费看| 久久国产精品人妻蜜桃| 久久热在线av| 成年人午夜在线观看视频| 12—13女人毛片做爰片一| 少妇被粗大的猛进出69影院| 久久久国产精品麻豆| 肉色欧美久久久久久久蜜桃| 人人妻,人人澡人人爽秒播| 亚洲中文av在线| 国产国语露脸激情在线看| 精品高清国产在线一区| 人人妻,人人澡人人爽秒播| 99精国产麻豆久久婷婷| 少妇人妻久久综合中文| 99九九在线精品视频| 丁香六月天网| 51午夜福利影视在线观看| av又黄又爽大尺度在线免费看| 女性被躁到高潮视频| 国产精品香港三级国产av潘金莲| av超薄肉色丝袜交足视频| 高清在线国产一区| 久久久精品区二区三区| 亚洲成av片中文字幕在线观看| 后天国语完整版免费观看| 日日夜夜操网爽| 女人被躁到高潮嗷嗷叫费观| 性高湖久久久久久久久免费观看| 视频区图区小说| 欧美乱码精品一区二区三区| 精品亚洲成国产av| 51午夜福利影视在线观看| 多毛熟女@视频| 99热全是精品| 亚洲男人天堂网一区| 亚洲一区二区三区欧美精品| 国产精品亚洲av一区麻豆| netflix在线观看网站| 老司机影院成人| 91成人精品电影| 十八禁网站网址无遮挡| 久久性视频一级片| 欧美日韩亚洲国产一区二区在线观看 | 淫妇啪啪啪对白视频 | 99国产精品免费福利视频| 一区二区三区精品91| 一级毛片电影观看| 午夜福利视频精品| 黄频高清免费视频| 在线观看免费午夜福利视频| 国产av国产精品国产| 久久久久国产精品人妻一区二区| 国产亚洲av片在线观看秒播厂| 久久久久精品人妻al黑| 性色av乱码一区二区三区2| 最新在线观看一区二区三区| 午夜福利视频精品| 午夜视频精品福利| kizo精华| 亚洲一码二码三码区别大吗| 少妇的丰满在线观看| 91九色精品人成在线观看| 国产亚洲精品第一综合不卡| 亚洲九九香蕉| 美女扒开内裤让男人捅视频| 青青草视频在线视频观看| 99热网站在线观看| 一区二区三区四区激情视频| 亚洲国产成人一精品久久久| 黄片小视频在线播放| 电影成人av| 精品国产一区二区久久| 午夜91福利影院| 男人添女人高潮全过程视频| 国产一区二区三区在线臀色熟女 | 久久午夜综合久久蜜桃| 国产精品欧美亚洲77777| 91精品国产国语对白视频| 少妇 在线观看| 亚洲三区欧美一区| 国产成人系列免费观看| 成年美女黄网站色视频大全免费| 热99re8久久精品国产| 欧美另类一区| 精品国产国语对白av| 一本色道久久久久久精品综合| 99九九在线精品视频| 国产日韩欧美在线精品| 国产一区二区在线观看av| 热re99久久精品国产66热6| 亚洲av欧美aⅴ国产| 美国免费a级毛片| av有码第一页| 久久青草综合色| 宅男免费午夜| 精品国产乱码久久久久久小说| 亚洲av成人一区二区三| 色综合欧美亚洲国产小说| 久久久久网色| 国产在线免费精品| 老司机午夜福利在线观看视频 | 多毛熟女@视频| 99久久人妻综合| 日韩大码丰满熟妇| 亚洲成人手机| 亚洲成人国产一区在线观看| 亚洲伊人久久精品综合| 亚洲专区国产一区二区| 久久精品人人爽人人爽视色| 日本猛色少妇xxxxx猛交久久| 亚洲中文字幕日韩| 亚洲精品中文字幕一二三四区 | 色播在线永久视频| 热99久久久久精品小说推荐| 操出白浆在线播放| 亚洲欧美日韩另类电影网站| 91成年电影在线观看| 男女高潮啪啪啪动态图| 妹子高潮喷水视频| 女人久久www免费人成看片| 欧美激情高清一区二区三区| 亚洲一卡2卡3卡4卡5卡精品中文| 制服人妻中文乱码| 亚洲精品第二区| 欧美激情极品国产一区二区三区| 日日夜夜操网爽| 国产精品二区激情视频| 色视频在线一区二区三区| 久久亚洲精品不卡| 夫妻午夜视频| 久久久欧美国产精品| 91麻豆av在线| 日韩电影二区| 国产亚洲欧美精品永久| 男女无遮挡免费网站观看| 老司机午夜福利在线观看视频 | 宅男免费午夜| 欧美久久黑人一区二区| 一二三四社区在线视频社区8| 欧美久久黑人一区二区| 宅男免费午夜| 人人妻人人添人人爽欧美一区卜| 高清视频免费观看一区二区| 久久久久久久精品精品| 99re6热这里在线精品视频| 两个人看的免费小视频| 中文字幕人妻丝袜一区二区| 一级毛片女人18水好多| 国产精品.久久久| 午夜免费鲁丝| 美女高潮到喷水免费观看| 精品国产乱子伦一区二区三区 | 国产免费一区二区三区四区乱码| 成在线人永久免费视频| 一本大道久久a久久精品| 久久人妻福利社区极品人妻图片| 久久精品亚洲熟妇少妇任你| 欧美精品一区二区大全| 精品少妇一区二区三区视频日本电影| 黄片大片在线免费观看| 黄色 视频免费看| 亚洲天堂av无毛| 国产精品久久久久久精品古装| 制服诱惑二区| 亚洲全国av大片| av福利片在线| 精品亚洲乱码少妇综合久久| 最新在线观看一区二区三区| 国产真人三级小视频在线观看| 色综合欧美亚洲国产小说| 美女午夜性视频免费| 亚洲综合色网址| 久久狼人影院| 国产高清国产精品国产三级| 日本精品一区二区三区蜜桃| 久久青草综合色| 欧美激情极品国产一区二区三区| 热re99久久国产66热| 成人亚洲精品一区在线观看| 日韩免费高清中文字幕av| 十八禁网站网址无遮挡| 大码成人一级视频| 欧美国产精品一级二级三级| 亚洲精品国产av成人精品| 亚洲男人天堂网一区| 一级毛片电影观看| 无遮挡黄片免费观看| 久久久精品区二区三区| 久久久国产欧美日韩av| 亚洲一码二码三码区别大吗| 99九九在线精品视频| 亚洲,欧美精品.| 成人国语在线视频| 999精品在线视频| 亚洲国产日韩一区二区| 亚洲国产av新网站| 亚洲精品国产av成人精品| 国产欧美日韩精品亚洲av| 亚洲情色 制服丝袜| 99国产精品免费福利视频| 高潮久久久久久久久久久不卡| 日韩,欧美,国产一区二区三区| 久久 成人 亚洲| 国产成人精品在线电影| 亚洲专区字幕在线| 久久久国产精品麻豆| 丰满人妻熟妇乱又伦精品不卡| 日韩 亚洲 欧美在线| 日韩 亚洲 欧美在线| 亚洲精品粉嫩美女一区| 法律面前人人平等表现在哪些方面 | 久久久久网色| 久久久精品免费免费高清| 18禁黄网站禁片午夜丰满| 黄色a级毛片大全视频| 男女无遮挡免费网站观看| 久久久精品国产亚洲av高清涩受| 精品人妻1区二区| 爱豆传媒免费全集在线观看| 久热这里只有精品99| 中文字幕最新亚洲高清| 日本黄色日本黄色录像| 一区二区三区乱码不卡18| 国产福利在线免费观看视频| 亚洲三区欧美一区| 人人妻人人澡人人看| 男女午夜视频在线观看| 欧美中文综合在线视频| 国产精品二区激情视频| 9191精品国产免费久久| 狠狠狠狠99中文字幕| 精品福利观看| 男人操女人黄网站| 免费在线观看影片大全网站| 嫁个100分男人电影在线观看| 窝窝影院91人妻| 热re99久久国产66热| 新久久久久国产一级毛片| 黄频高清免费视频| 日日摸夜夜添夜夜添小说| 国产主播在线观看一区二区| 国产精品秋霞免费鲁丝片| 婷婷成人精品国产| 少妇裸体淫交视频免费看高清 | 一级毛片女人18水好多| 色94色欧美一区二区| 欧美 亚洲 国产 日韩一| 电影成人av| 999久久久精品免费观看国产| 亚洲av成人不卡在线观看播放网 | 男人添女人高潮全过程视频| 两个人看的免费小视频| 人人妻人人澡人人爽人人夜夜| 老司机深夜福利视频在线观看 | 男人舔女人的私密视频| 首页视频小说图片口味搜索| 丝瓜视频免费看黄片| 精品欧美一区二区三区在线| 欧美激情 高清一区二区三区| 日日摸夜夜添夜夜添小说| 男女床上黄色一级片免费看| 日韩大片免费观看网站| 成年女人毛片免费观看观看9 | 午夜福利乱码中文字幕| 日本wwww免费看| 每晚都被弄得嗷嗷叫到高潮| 久久久久国内视频| 国产精品免费视频内射| 国产精品久久久久久人妻精品电影 | 久久综合国产亚洲精品| 午夜福利在线免费观看网站| 午夜激情av网站| 国产av又大| 午夜免费成人在线视频| 国产激情久久老熟女| 亚洲男人天堂网一区| 色播在线永久视频| a级毛片在线看网站| 免费在线观看影片大全网站| 爱豆传媒免费全集在线观看| 亚洲,欧美精品.| 大片免费播放器 马上看| 满18在线观看网站| 欧美亚洲日本最大视频资源| 18禁裸乳无遮挡动漫免费视频| 久久毛片免费看一区二区三区| 啦啦啦中文免费视频观看日本| 国产有黄有色有爽视频| 国产视频一区二区在线看| 免费在线观看黄色视频的| 免费少妇av软件| 国产亚洲午夜精品一区二区久久| 涩涩av久久男人的天堂| 国产极品粉嫩免费观看在线| 欧美亚洲日本最大视频资源| 欧美精品一区二区免费开放| 精品久久久精品久久久| 国产免费av片在线观看野外av| 精品一区二区三区四区五区乱码| 丝袜在线中文字幕| 久久久国产成人免费| 国产一级毛片在线| 免费观看a级毛片全部| 亚洲精品中文字幕在线视频| 亚洲天堂av无毛| 欧美 日韩 精品 国产| 各种免费的搞黄视频| 欧美激情 高清一区二区三区| 精品国产国语对白av| 男女边摸边吃奶| 久久人人97超碰香蕉20202| 国产亚洲欧美精品永久| 在线天堂中文资源库| 免费在线观看视频国产中文字幕亚洲 | 91精品伊人久久大香线蕉| 国产精品久久久久久人妻精品电影 | 日韩欧美国产一区二区入口| 亚洲欧洲日产国产| 国产在线免费精品| 欧美日韩精品网址| 操出白浆在线播放| 欧美黑人精品巨大| 精品少妇内射三级| 1024视频免费在线观看| 嫩草影视91久久| 一级毛片电影观看| 精品久久久久久久毛片微露脸 | 欧美xxⅹ黑人| 另类精品久久| 大码成人一级视频| 日韩三级视频一区二区三区| 曰老女人黄片| 免费在线观看黄色视频的| 日韩精品免费视频一区二区三区| 啦啦啦视频在线资源免费观看| 俄罗斯特黄特色一大片| 一本大道久久a久久精品| 亚洲国产精品999| 久久国产精品大桥未久av| 搡老乐熟女国产| 两个人免费观看高清视频| 久久ye,这里只有精品| 久久久欧美国产精品| 熟女少妇亚洲综合色aaa.| 青春草视频在线免费观看| av网站在线播放免费| 两个人免费观看高清视频| 青春草亚洲视频在线观看| 亚洲男人天堂网一区| 亚洲 欧美一区二区三区| 成人亚洲精品一区在线观看| av电影中文网址| 高清黄色对白视频在线免费看| 午夜激情av网站| av福利片在线| av视频免费观看在线观看| 丝袜美足系列| 免费人妻精品一区二区三区视频| 免费黄频网站在线观看国产| 亚洲精品久久午夜乱码| av网站免费在线观看视频| 日韩制服丝袜自拍偷拍| 伦理电影免费视频| 少妇裸体淫交视频免费看高清 | 亚洲性夜色夜夜综合| 一区二区三区激情视频| 欧美一级毛片孕妇| 女人久久www免费人成看片| 国产在视频线精品| 满18在线观看网站| 97精品久久久久久久久久精品| 老司机靠b影院| 午夜福利在线观看吧| 亚洲中文字幕日韩| av福利片在线| 97在线人人人人妻| 一本色道久久久久久精品综合| 美女午夜性视频免费| 久久午夜综合久久蜜桃| 大香蕉久久网| 欧美乱码精品一区二区三区| 最近最新中文字幕大全免费视频| 秋霞在线观看毛片| av福利片在线| 亚洲国产毛片av蜜桃av| 欧美中文综合在线视频| 91精品三级在线观看| 99久久综合免费| 亚洲人成电影观看| 免费少妇av软件| 精品高清国产在线一区| 在线观看免费日韩欧美大片| 日韩视频在线欧美| 丝袜美足系列| 欧美成狂野欧美在线观看| 男女无遮挡免费网站观看| 99久久99久久久精品蜜桃| 电影成人av| 在线观看人妻少妇| 一区二区三区四区激情视频| 国产黄色免费在线视频| 天天影视国产精品| 99久久综合免费| 国产一区二区三区av在线| 久久精品aⅴ一区二区三区四区| 精品国产国语对白av| 男人操女人黄网站| 另类亚洲欧美激情| 国产精品熟女久久久久浪| 久久久精品国产亚洲av高清涩受| 一本一本久久a久久精品综合妖精| 国产精品久久久久久精品电影小说| 国产精品秋霞免费鲁丝片| 日本精品一区二区三区蜜桃| 水蜜桃什么品种好|