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

    Lagrangian methods for water transport processes in a long-narrow bay-Xiangshan Bay, China*

    2014-04-05 21:44:04LIANGShuxiu梁書秀HANSonglin韓松林SUNZhaochen孫昭晨
    關(guān)鍵詞:梁書松林

    LIANG Shu-xiu (梁書秀), HAN Song-lin (韓松林), SUN Zhao-chen (孫昭晨)

    State Key Laboratory of Coastal and Offshore Engineering, Dalian University of Technology, Dalian 116024, China, E-mail: sxliang@dlut.edu.cn

    HU Zhan-ming (胡展銘)

    National Marine Environmental Monitoring Center, Dalian 116023, China

    Lagrangian methods for water transport processes in a long-narrow bay-Xiangshan Bay, China*

    LIANG Shu-xiu (梁書秀), HAN Song-lin (韓松林), SUN Zhao-chen (孫昭晨)

    State Key Laboratory of Coastal and Offshore Engineering, Dalian University of Technology, Dalian 116024, China, E-mail: sxliang@dlut.edu.cn

    HU Zhan-ming (胡展銘)

    National Marine Environmental Monitoring Center, Dalian 116023, China

    (Received July 10, 2013, Revised December 25, 2013)

    A better understanding of water transport processes is highly desirable for the exploitation of the ocean resources and the protection of the ocean ecological system. In this paper, the Lagrangian methods are used to study the water transport processes in Xiangshan Bay in China, a typical semi-closed and narrow-shaped bay with complex coastline and topography. A high-resolution 3-D hydrodynamic model is developed and verified, and the results from the model agree well with the field data. Based on the hydrodynamic model, the Lagrangian residual current is computed by using the particle tracking method. A concept based on the dynamical systems theory, the Lagrangian coherent structures (LCSs), is introduced to uncover the underlying structures which act as the transport barriers in the flow. The finite-time Lyapunov exponent (FTLE) fields are computed from the hydrodynamic model results to extract the LCSs. The results indicate that the LCSs act as the internal structures of the Lagrangian residual current and the Lagrangian residual current displays the residual current speed and direction of different water regimes separated by the LCSs. The water masses with different transport characteristics can be identified and their exchange ability with other water masses can be estimated by combining the Lagrangian particle tracking with the LCSs methods. The comprehensive applications of these Lagrangian methods reveal the underlying structures and the inhomogeneous characteristics of the water transport in Xiangshan Bay.

    water transport processes, circulation, Lagrangian coherent structures (LCSs), Xiangshan Bay

    Introduction

    In recent years, with the rapid development of the harbor industry, the fish culture and the social economy in coastal areas, the coastal water pollution becomes a very serious issue due to the excess and unreasonable exploitations. The overabundance of pollution matters in the water, such as the organic matter, the DIP, the DIN and the COD, etc., may cause eutrophication and even red tide. The pollution in the water could be diluted to a low concentration or transported to open seas by the seawater movement. However, the pollution may accumulate or fluctuate in the bay because of complex topography, weak current conditions or other environment factors. Therefore, it is necessary to study the water transport characteristics in order to better control the water quality.

    The Lagrangian and Eulerian methods are widely applied to analyze the water transport problems, e.g., the Eulerian and Lagrangian residual current, the Lagrangian particle tracking, the dye tracer study, the drift experiment, etc.. The Eulerian residual current could be obtained by averaging the velocity at a fixed position for one or several tidal cycles or by using other filtering methods. It was widely used to analyze the observed tidal currents and numerical simulation results[1,2]. However, Longuet-Higgins[3]pointed out that the mass transport velocity is not solely controlledby the mean velocity, but also determined by the Stokes drift velocity. The sum of the Eulerian mean velocity and the Stokes drift velocity is called the first-order approximation of the Lagrangian residual velocity[4,5], which was successfully applied in weakly nonlinear systems[6]. In a strongly nonlinear tidal system, the Lagrangian residual current can be calculated as the net displacement over one or more tidal cycles by the Lagrangian particle tracking method. The trajectory analysis was also proved to be a valuable tool for determining the origin as well as the fate of specific water masses[7]. As an incorporation of the Lagrangian method, the Lagrangian coherent structures (LCSs) can be used to study some underlying flow structures that are not evident from the Eulerian field because of the chaotic horizontal transport present in many oceanic flows[8,9]. The LCSs of the flow field are defined as the attracting or repelling material curves in the finite-time Lyapunov exponent (FTLE) fields. This method was applied in several transport problems. For example, Olascoaga et al.[10]used the LCSs to study the spreading of plankton blooms on the West Florida Shelf. Lekien et al.[11]applied this method to predict the transport of contaminants and pollutions. It was also been used to study the causes of the low quality coastal water near the shoreline[12].

    The purpose of the present study is to study the water transport processes in a complex topography bay, Xiangshan Bay in China by using the Lagrangian methods, and to provide some food for thought for the choice of adequate tools in estimating the water transport. In this study, a 3-D high resolution hydrodynamic model is implemented to simulate the water movement in Xiangshan Bay. The Lagrangian residual current and the LCSs are computed based on the Lagrangian particle tracking method. Furthermore, particle tracking experiments were conducted to reveal the effect of the residual current and the LCSs on the water mass transport and the water exchange.

    1. Study area

    Xiangshan Bay is a semi-closed narrow-shaped bay with complex topography, located in the middle coast of Zhejiang Province between 121o25′E-122o00′E and 29o23′N-29o49′N, as shown in Fig.1(a). The bay is 62.8 km long from the bay mouth to the bay head, and its width varies from 20 km near the mouth to a minimum of 3 km in the middle of the bay. Liuheng Island locates at the mouth of the embayment which divides the entrance into Niubishan channel in the southeast and Fodu channel in the2northeast. The water area of Xiangshan Bay is 563 km, and the averaged depth is about 10 m, and the maximum depth is 70 m. Figure 1(b) is the contour map of the water depth D in Xiangshan Bay.

    Xiangshan Bay, the largest aquaculture base of Zhejiang Province, enjoys an excellent aquaculture condition. Since 1990, the marine cage culture saw a great development and became the main aquaculture mode in the bay. In 2005, the fish cages in the bay were increased to 65 000. The main farmed species are the yellow croaker, the weever, and the sciaenops ocellatus. Recent years witnessed a rapid development of the harbor construction and the harbor industry, such as the port logistics, the processing centers, the Wushashan and Guohua power station. Along with the rapid development of the cage culture and the harbor industry in the bay, a large amount of industrial and agricultural wastewater has been discharged into the bay. The eutrophication has become a serious problem and the red tide happens occasionally, as a serious constraint to the development of aquaculture and tourism. From previous studies[13], it is known that the water exchange capacity gradually drops from the mouth to the head of the bay. The period of the 90% water exchange is about 80 d near the bay head areas. A better understanding of the water transport is therefore of great importance in Xiangshan Bay. For the convenience of research, Xiangshan Bay is subdivided into four regions according to the hydrological and water quality investigations, as shown in Fig.1(a). These regions are the intersection area of two channels (I), the mouth of Xiangshan Bay (II), the middle and head of the bay (III), (IV).

    2. Methods

    2.1 3-D hydrodynamic model

    A 3-D finite volume coastal ocean model (FVCOM, 2.7.1) is selected for the Xiangshan Bay application, which was developed originally by Chen et al.[14], based on the unstructured-grid, the finite-volume, and the 3-D primitive equations. The model governing equations consist of the momentum, the continuity, the temperature, the salinity, and the density equations. In the horizontal direction, a non-overlapping unstructured triangular grid is used to resolve the dynamics in the regions with complex shorelines and the sigma-coordinate transformation is used in the vertical direction to accommodate the irregular variable bottom topography. The equations are solved numerically by the flux calculation over an arbitrarilysized triangular mesh using the finite-volume approach. The finite-volume approach is better suited to guarantee the mass conservation in both the individual control element and the entire computation domain. The shape of Xiangshan Bay is complex, with numerous islands in it. Therefore, the FVCOM model is one of the best choices.

    The Smagorinsky eddy parameterization method and the MY-2.5 model are used to determine the horizontal and vertical diffusions for the momentum, the temperature and the salinity. The horizontal eddy viscosity coefficient is calculated by

    whereu,vandωare the components of the velocity in the horizontal(x,y)and vertical(σ)directions, respectively,C is the constant parameter,?ζis the area of the individual tracer control element. Then the thermal vertical eddy diffusion coefficient, AH=PrAM,Pris the turbulence Prandtl number.

    The equations of the MY-2.5 turbulence closure model are given by

    where q2is the turbulent kinetic energy,lis the turbulent macroscale.gis the gravitational acceleration,andis a wall proximity function. The von Karman constant k =0.4.E1,E2,E3are the close constants,,represent the horizontal diffusion of the turbulent kinetic energy and the macroscale,,are the vertical eddy viscosity coefficient and the thermal vertical eddy diffusion coefficient,is the vertical eddy diffusion coefficient of the turbulent kinetic energy. Turbulent close parameters are given by

    2.2 Lagrangian residual current

    The Lagrangian residual current ULis defined as the net displacement vector of the water column over a periodTand is calculated by the Lagrangian particle tracking method. It is expressed as

    where xTand xoare the end and start positions of the particles over the periodTandT is equal to 25 h in this study. With the integrated periodT equal to 25 h, it includes approximately two tidal circles and the averaged short-term characteristics of the residual current can be depicted. The particles in the water are tracked by solving the equation,

    It can be solved by using the fourth order Runge-Kutta scheme, as derived from solving the discrete integral

    2.3 Lagrangian coherent structures

    where the second equality comes from taking the Taylor series expansion of the flow about point x. Since δx(t)is infinitesimal, the term ο[δx(t)2] can be considered negligible. Therefore, the magnitude of the distance is given by (using the standard vector L2-norm)

    When the convention thatM?denotes the adjoint (transpose) of M are used, the symmetric matrix Δ is written as

    3. Model configuration and verification

    The model is implemented by two nested grids with different resolutions. The larger grid covers almost the entire Zhejiang coastal areas and there are 165 450 triangle meshes in the computation area (Fig.2). The largest element scales 3 500 m. The model boundaries are determined by the observed water elevation data in 2009 from Ocean Environment Monitoring Centre of Ningbo (OEMCN). The detailed information could be found in Xiong et al.[16].

    The study area covers Xiangshan Bay as its center, where the high resolution grid is used. It includes Xiangshan Bay, Fodu channel and Niubishan Channel surrounded by white lines and coastlines as shown in Fig.2. There are 18 645 elements and 9 958 nodes in the study area. The element scales vary from 80 m for complex coastlines and islands to 400 m for offshore areas. Eleven sigma layers are adopted vertically,which results in a vertical resolution of between 0.3 m to 5 m in most areas. The boundary conditions are extracted from a large model. The meteorological parameters (the wind components, the air temperature at 2 m, the shortwave radiation, the cloudiness etc.) are obtained from the analysis of the National Center for Environmental Prediction (NCEP), with a bilinear interpolation in space and a linear interpolation in time. The net heat flux at the air-sea interface can be calculated with these parameters.

    4. Results

    4.1 The Lagrangian residual current in Xiangshan Bay

    The hydrodynamic model is baroclinic, with con-siderations of the tide, the wind, the river runoff, and the surface heat exchange. Hence the Lagrangian residual currents computed in this study consist of the tide-induced, river-, wind-, and density-driven currents. Figs.6(a) and 6(b) show the distribution of the Lagrangian residual current field at the surface and bottom layers. In the Lagrangian residual currents there are few eddies or circulations because of the narrow-long shape of Xiangshan Bay. The direction of the Lagrangian residual currents is mainly parallel with the shoreline in the fiord. The surface Lagrangian residual current from Fodu channel flows to the Xiangshan Bay fiord and sees an intense outgoing current in the Xiangshan Bay fiord. The outgoing current meets the current from Fodu channel at the mouth of the bay and flows to Niubishan channel from the south area of the Xiangshan Bay mouth. The Lagrangian residual speed is 0.03 to 0.18 m/s in the fiord and 0.12 m/s to 0.3 m/s in the region (I). This means that the water column net displacement at the mouth of the bay is larger than that in the fiord. In the vertical direction, the Lagrangian residual current at the bottom layer is distinctly different, especially in the fiord (Fig.6(b)). In most areas of region (I), the Lagrangian residual current speed decreases from surface to bottom. The bottom Lagrangian residual current from Fodu channel is divided into two branches. One flows to the Xiangshan Bay fiord and the other flows to Niubishan channel. In the fiord, it shows an incoming current and the current speed is about 0.03 m/s-0.12 m/s. At both surface and bottom layers, the Lagrangian residual current speed decreases gradually from the mouth to the head of the bay. The two-layered circulation in the vertical direction is mainly due to the density stratification in the fiord[1].

    4.2 The Lagrangian coherent structures in Xiangshan Bay

    In this section, the FTLE fields in Xiangshan Bay are calculated by using the surface velocity simulated by the 3-D hydrodynamic model. The FTLE fields reflect the Lagrangian characteristics of the flow fields in an interval time t-. The meaningful LCS can be obtained with t-t0close to a typical time scale of the flow[9]. The choice of t-t0is based on the flow characteristics in the study area. The tidal current in our study area is rectilinear and rotary current with a period of about 12 h. Sot-t0could be chosen as one or several tidal current periods. A series of experiments wheret-t0=6 h, 12 h, 24 h, 48 h, respectively, are performed to estimate the optimal t-t0. The results show that the internal structures are already clear and typical fort-t0=12 h. For a shorter interval time (t-t0=6h), the FTLE value is relatively small and the LCSs have not developed fully, because the particles are not separated from each other. For a longer interval time, the structures are similar but some small fragile structures grow. To display the typical structures clearly,t-t0=12 h is taken for further analyses.

    5. Discussions

    In the definitions of the Lagrangian residual current and the LCSs, the former represents the mean velocity of the water particles and the latter represent the separation rates of the adjacent water particles during a given interval time. The both methods are based on the Lagrangian particle tracking method and are dependent on the initial tidal current phase. The relationships between them are explored in this paper by comparing the results of the Lagrangian residual current and the LCSs, computed at the same time. The results are shown in Figs.6(a) and 7(c) and the initial tidal current phase is at a low slack water. The outgoing current in the fiord and the water from Foduchannel meet at the mouth area of the bay and the directions of them are opposite (see Fig.6(a)). This may cause the water particles flowing in different directions and increase the separation rate. As expected, there is a material line with a high separation rate appears at the dividing line of the two currents (see Fig.7(c)). The Lagrangian residual currents in different water regimes are separated by the LCSs with similar residual current speeds and directions. The LCS in the region (II) acts as the separatrix of the outgoing current in the fiord and the water from Fodu channel. The Lagrangian residual current in the north of the LCS flows to the fiord as an outgoing current for the south region. With the two Lagrangian methods, the water transport characteristics are revealed from different viewpoints. The LCSs act as the internal structure of the Lagrangian residual current.

    The pollution in Xiangshan Bay is mainly transported out of the bay in the ebb tidal current. The circulation at the mouth of the bay is crucial for the water exchange with open seas. From the Lagrangian residual current and the LCSs in the Section 5, it can be seen that the circulation is relatively complex in the regions (I) and (II). The existence of the LCSs governs the transport and it is verified numerically that the flux of the particles through these distinguished lines could be negligible[8]. Particle tracking experiments are conducted further to reveal the relationship among the Lagrangian residual current, the LCSs and the water exchange characteristics in different regimes.Two parts of artificial particles are assigned in the region (II), which are initially located in both sides of the LCS (see Fig.7(a)) as shown in the top left panel of Fig.8 and different colors are used to distinguish them.

    The particles are released at high tide and will be tracked for one tidal current period. The locations of the particles are given every 2 h in Fig.8. It is shown that almost no mixing occurs between the two groups of particles within 12 h. Analysis of Fig.8 shows that few particles initially in north of the LCS (gray color particles) are transported out of the region (II). In contrast, a large number of particles initially in south of the LCS (black color particles) are flushed out. It is of a particular concern that hardly any particles are transported out from Fodu channel, even for the particles in north of the LCS. After 12 h, the gray color particles return back to the region (II) and most of them are also in north of the LCS. In order to quantify the number of the particles that flush out of the region (II) and the time it takes, the particles are tracked until the percentage of particles retained in Xiangshan Bay (all the regions except (I)) takes a relatively stable value. Figure 9 shows the time series of the percentage P of the particles remained in Xiangshan Bay. It is shown that 87% of the particles originally located in north of the LCS still stay in Xiangshan Bay, but the percentage is only 16% for the particles located in south of the LCS. The result of another couple of particles which are continuously released for 12 h is shown in Fig.9 with the line C and D. It is also shown that the percentage of the particles retained in Xiangshan Bay for the particles originally located in the north of the LCS is larger than that in the south. It reveals that the water exchange rate is much higher in the south region of the LCS.

    6. Conclusions

    In the present paper, a high resolution, 3-D hydrodynamic model is developed to simulate the water movements. Based on the hydrodynamic model, the Lagrangian residual current, the Lagrangian particle tracking method and the LCSs are used to analyze the water transport processes in Xiangshan Bay. The main conclusions are as follows:

    (1) The distribution of the Lagrangian residual current shows that it is an intense outgoing current in the Xiangshan Bay fiord at the surface layer and an incoming current at the bottom layer. At both surface and bottom layers, the Lagrangian residual current speed decreases gradually from the mouth to the head of the bay. The LCSs are the moving separatrixes of different water masses and they evolve periodically with time as the ebb and flood tidal current oscillates. The period is the same as the period of the dominant tidal component in the study area.

    (2) The distributions of the Lagrangian residual current and the LCSs are both dependent on the initial tidal current phase. With the two Lagrangian methods, the water transport characteristics are obtained from different points of view. The LCSs act as the transport barriers in the bay and the internal structures of the Lagrangian residual current. The Lagrangian residual current displays the residual current speed and direction of different water regimes separated by the LCSs.

    (3) The water masses with different transport characteristics can be distinguished by the LCSs and their exchange ability with other water masses can be estimated by the Lagrangian particle tracking methods. It is shown that some obvious LCSs exist at the mouth area of the bay and the water exchange abilities are different among these water regimes separated by the LCSs. The water exchange rate is higher in the south of the LCSs area than in the north area.

    This research demonstrates that the LCSs from the dynamical system theory is a useful tool to understand the water transport and mixing processes. The comprehensive applications of the LCSs and other Lagrangian methods provide a rich body of knowledge of the water transport characteristics in the bay, which can serve as a scientific basis for the control and protection of water quality.

    [1] DONG Li-xian, SU Ji-lan. Salinity distribution and mixing in Xiangshangang Bay, I. Salinity distribution and circulation pattern[J]. Oceanologia et Limnologia Sinica, 2000, 31(2): 151-158(in Chinese).

    [2] MURPHY P. L., VALLE-LEVINSON A. Tidal and residual circulation in the St. Andrew Bay system, Florida[J]. Continental Shelf Research, 2008, 28(19): 2678-2688.

    [3] LONGUET-HIGGINS M. S. On the transport of mass by time-varying ocean currents[J]. Deep Sea Research and Oceanographic, 1969, 16(5): 431-447.

    [4] LIU G., LIU Z. and GAO H. et al. Simulation of the Lagrangian tide-induced residual velocity in a tide-dominated coastal system: A case study of Jiaozhou Bay, China[J]. Ocean Dynamics, 2012, 62(10-12): 1443-1456.

    [5] JIANG W., FENG S. Analytical solution for the tidally induced Lagrangian residual current in a narrow bay[J]. Ocean Dynamics, 2011, 61(4): 543-558.

    [6] MULLER H., BLANKE B. and DUMAS F. et al. Estimating the Lagrangian residual circulation in the Iroise Sea[J]. Journal of Marine Systems, 2009, 78(1): S17-S36.

    [7] JONSSON B., LUNDBERG P. and DOOS K. Baltic sub-basin turnover times examined using the Rossby Centre Ocean Model[J]. Ambio, 2004, 33(4-5): 257-260.

    [8] SHADDEN S. C., LEKIEN F. and MARSDEN J. E. Definition and properties of Lagrangian coherent structures from finite-time Lyapunov exponents in two-dimensional aperiodic flows[J]. Physica D: Nonlinear Phenomena, 2005, 212(3): 271-304.

    [9] HUHN F., Von KAMEKE A. and ALLEN-PERKINS S. et al. Horizontal Lagrangian transport in a tidal-driven estuary-Transport barriers attached to prominent coastal boundaries[J]. Continental Shelf Research, 2012, 39(1): 1-13.

    [10] OLASCOAGA M. J., RYPINA I. I. and BROWN M. G. et al. Persistent transport barrier on the West Florida Shelf[J]. Geophysical research letters, 2006, 33(22): L22603.

    [11] LEKIEN F., COULLIETTE C. and MARIANO A. J. et al. Pollution release tied to invariant manifolds: A case study for the coast of Florida[J]. Physica D: Nonlinear Phenomena, 2005, 210(1): 1-20.

    [12] FIORENTINO L. A., OLASCOAGA M. J. and RENIERS A. et al. Using Lagrangian Coherent Structures to understand coastal water quality[J]. Continental Shelf Research, 2012, 47(1): 145-149.

    [13] DONG Li-xian, SU Ji-lan. Numerical study of water exchange in Xiangshan Bay, II. Model application and water exchange study[J]. Oceanlogia et limnologia sinica, 1999, 30(5): 465-470(in Chinese).

    [14] CHEN C., LIU H. and BEARDSLEY R. C. An unstructured grid, finite-volume, three-dimensional, primitive equations ocean model: Application to coastal ocean and estuaries[J]. Journal of Atmospheric and Oceanic Technology, 2003, 20(1): 159-186.

    [15] BRANICKI M., WIGGINS S. Finite-time Lagrangian transport analysis: Stable and unstable manifolds of hyperbolic trajectories and finite-time Lyapunov exponents[J]. Nonlinear Processes in Geophysics, 2010, 17(1): 1-36.

    [16] XIONG Wei, LIU Bi-jin and SUN Zhao-chen et al. 3D numerical simulation of tide and tidal currents in sea adjacent to Ningbo and Zhoushan[J]. Journal of Waterway and Harbor, 2011, 32(6): 399-407(in Chinese).

    [17] HAN S., LIANG S. and SUN Z. Study of pollution transport based on a dynamical system theory in Xiangshan Bay, China[C]. ASME 2013 32nd International Conference on Ocean, Offshore and Arctic Engineering. Nantes, France, 2013, 1-8.

    10.1016/S1001-6058(14)60063-9

    * Project supported by the National Natural Science Foundation of China (Grant No. 51279028), the Foundation for Innovative Research Groups of the National Natural Science Foundation of China (Grant No. 51221961) and the Public Welfare Projects of China’s Oceanic Administration (Grant Nos. 200805086, 201105009).

    Biography: LIANG Shu-xiu (1972-), Female, Ph. D., Associate professor

    SUN Zhao-chen,

    E-mail: sunzc@dlut.edu.cn

    猜你喜歡
    梁書松林
    Implementation and application of PyNE sub-voxel R2S for shutdown dose rate analysis
    寫給舅舅的一封信
    回答
    沉痛悼念霍松林名譽(yù)會(huì)長
    中華詩詞(2017年3期)2017-11-27 03:44:22
    A GPU accelerated finite volume coastal ocean model*
    我的爺爺和長征
    霍松林詩詞選
    松林村
    點(diǎn)絳唇·詠風(fēng)蘭
    Experimental study on the time-dependent dynamic mechanical behaviour of C60 concrete under high-temperatures
    中文字幕人妻熟人妻熟丝袜美| 国产69精品久久久久777片| 麻豆av噜噜一区二区三区| 成人美女网站在线观看视频| 精品福利观看| 日本与韩国留学比较| 日韩欧美一区二区三区在线观看| 国产一区二区三区在线臀色熟女| 久久精品国产亚洲网站| 国产精品不卡视频一区二区| 日本五十路高清| 黑人高潮一二区| 精品久久久久久久末码| 午夜精品国产一区二区电影 | 在线播放国产精品三级| 国产精品乱码一区二三区的特点| 国产免费一级a男人的天堂| 少妇丰满av| 在线观看美女被高潮喷水网站| 成人性生交大片免费视频hd| 亚洲av成人精品一区久久| 欧美激情久久久久久爽电影| 久久精品夜夜夜夜夜久久蜜豆| 综合色丁香网| 99热这里只有是精品在线观看| 日日摸夜夜添夜夜添av毛片| 国产三级在线视频| 久久久久久大精品| 亚洲成a人片在线一区二区| 午夜视频国产福利| 老熟妇仑乱视频hdxx| 国产一区亚洲一区在线观看| 插阴视频在线观看视频| 丝袜喷水一区| 我的老师免费观看完整版| 高清日韩中文字幕在线| 亚洲四区av| 精品久久久久久成人av| 成人高潮视频无遮挡免费网站| 欧美+亚洲+日韩+国产| 国产av不卡久久| 国产亚洲精品久久久com| 人妻久久中文字幕网| 亚洲性久久影院| 亚洲精品一区av在线观看| 国产高清视频在线观看网站| 国产精品一区www在线观看| 亚洲人成网站高清观看| 成人毛片a级毛片在线播放| 深夜精品福利| 免费观看的影片在线观看| 精品久久久久久久久久久久久| 老熟妇乱子伦视频在线观看| 久久久久久久久大av| 国产精品久久视频播放| 老师上课跳d突然被开到最大视频| 久久久成人免费电影| 淫秽高清视频在线观看| 国产v大片淫在线免费观看| 免费av观看视频| 亚洲激情五月婷婷啪啪| 国国产精品蜜臀av免费| 非洲黑人性xxxx精品又粗又长| 国产精品野战在线观看| 国产男靠女视频免费网站| 男女边吃奶边做爰视频| 日韩av在线大香蕉| 亚洲av一区综合| 国产精品,欧美在线| 国产欧美日韩精品亚洲av| 搡老妇女老女人老熟妇| 色播亚洲综合网| 麻豆国产av国片精品| 女人被狂操c到高潮| av黄色大香蕉| 欧美绝顶高潮抽搐喷水| 菩萨蛮人人尽说江南好唐韦庄 | 少妇被粗大猛烈的视频| 男插女下体视频免费在线播放| 久久久久久久久久久丰满| 亚洲专区国产一区二区| 日本精品一区二区三区蜜桃| 国产人妻一区二区三区在| 亚洲丝袜综合中文字幕| 国语自产精品视频在线第100页| 全区人妻精品视频| 男女下面进入的视频免费午夜| 91在线精品国自产拍蜜月| 精品午夜福利视频在线观看一区| 国产老妇女一区| 国产精品一及| 少妇人妻精品综合一区二区 | 久久久久免费精品人妻一区二区| 最新中文字幕久久久久| 国产精品久久视频播放| 真人做人爱边吃奶动态| 国产三级在线视频| 国产成人a∨麻豆精品| 别揉我奶头 嗯啊视频| 在线观看66精品国产| 久久精品国产清高在天天线| 亚洲成人久久爱视频| 日日摸夜夜添夜夜添小说| 日韩一本色道免费dvd| 成人特级av手机在线观看| 久久人人爽人人片av| 色5月婷婷丁香| 国产午夜精品久久久久久一区二区三区 | 亚洲精品乱码久久久v下载方式| 蜜桃久久精品国产亚洲av| 成人永久免费在线观看视频| 蜜臀久久99精品久久宅男| 观看美女的网站| 国产91av在线免费观看| 国产极品精品免费视频能看的| 亚洲av成人av| 99热这里只有是精品在线观看| 亚洲最大成人中文| 免费不卡的大黄色大毛片视频在线观看 | 丝袜喷水一区| 国产三级在线视频| 亚洲国产精品成人综合色| 久久精品国产亚洲av天美| 亚洲专区国产一区二区| 女生性感内裤真人,穿戴方法视频| 99热6这里只有精品| 免费av观看视频| 五月伊人婷婷丁香| 一本一本综合久久| 亚洲最大成人中文| 亚洲av成人av| 欧美另类亚洲清纯唯美| 成年女人看的毛片在线观看| 国产一区二区激情短视频| 两个人视频免费观看高清| 欧美成人精品欧美一级黄| 国产伦精品一区二区三区四那| 亚洲精品一卡2卡三卡4卡5卡| 国产精品免费一区二区三区在线| 欧美色视频一区免费| 亚洲人成网站在线观看播放| 村上凉子中文字幕在线| 床上黄色一级片| 校园人妻丝袜中文字幕| 日韩欧美在线乱码| 国产高潮美女av| 免费电影在线观看免费观看| 久久久久国产网址| 成年女人看的毛片在线观看| 精品午夜福利在线看| 国产成人影院久久av| 亚洲精品乱码久久久v下载方式| 亚洲成av人片在线播放无| 小说图片视频综合网站| 亚洲欧美日韩无卡精品| 成人漫画全彩无遮挡| 午夜福利视频1000在线观看| 久久人人爽人人片av| 久久久国产成人精品二区| 天堂网av新在线| 精品一区二区三区视频在线| 亚洲国产色片| 国产精品久久电影中文字幕| 大又大粗又爽又黄少妇毛片口| 黄色日韩在线| 久久国产乱子免费精品| 草草在线视频免费看| 免费黄网站久久成人精品| 久久精品国产亚洲av天美| a级一级毛片免费在线观看| 大香蕉久久网| 亚洲18禁久久av| 桃色一区二区三区在线观看| 九九在线视频观看精品| 国产精品爽爽va在线观看网站| 欧美成人精品欧美一级黄| 亚洲成人av在线免费| 变态另类丝袜制服| 热99re8久久精品国产| 国产又黄又爽又无遮挡在线| 国产成人a区在线观看| av专区在线播放| 三级经典国产精品| 十八禁网站免费在线| 99热这里只有是精品在线观看| 免费av不卡在线播放| 久久精品影院6| 国产69精品久久久久777片| 在线观看66精品国产| 亚洲内射少妇av| 免费在线观看成人毛片| 午夜精品一区二区三区免费看| 日韩成人伦理影院| 午夜福利在线观看免费完整高清在 | 精品一区二区三区av网在线观看| 男人舔奶头视频| 99国产精品一区二区蜜桃av| 亚洲av成人精品一区久久| 久久精品人妻少妇| 天堂√8在线中文| 一本久久中文字幕| 久久久久性生活片| 中文在线观看免费www的网站| 亚洲色图av天堂| 九九久久精品国产亚洲av麻豆| 亚洲精品国产成人久久av| 亚洲精品影视一区二区三区av| 最近在线观看免费完整版| 男人舔女人下体高潮全视频| 一级毛片久久久久久久久女| 91在线精品国自产拍蜜月| 美女内射精品一级片tv| 免费黄网站久久成人精品| 网址你懂的国产日韩在线| 男女下面进入的视频免费午夜| 亚洲国产欧美人成| 偷拍熟女少妇极品色| 禁无遮挡网站| 午夜精品国产一区二区电影 | 国产高清激情床上av| 狂野欧美白嫩少妇大欣赏| 久久久久国产精品人妻aⅴ院| 变态另类丝袜制服| 国产综合懂色| eeuss影院久久| 久99久视频精品免费| 一个人看的www免费观看视频| 欧美性猛交黑人性爽| 一级av片app| 国产高清不卡午夜福利| 国产精品精品国产色婷婷| 丝袜喷水一区| 精品久久久久久久久亚洲| 国产成人影院久久av| 伦精品一区二区三区| 搡老岳熟女国产| 国产精品一区www在线观看| 91狼人影院| 丰满人妻一区二区三区视频av| 桃色一区二区三区在线观看| 日韩一本色道免费dvd| 性插视频无遮挡在线免费观看| 麻豆一二三区av精品| 成人三级黄色视频| 天天躁日日操中文字幕| 久久久国产成人精品二区| 看片在线看免费视频| 国产爱豆传媒在线观看| 国产伦精品一区二区三区视频9| 久久久精品欧美日韩精品| 51国产日韩欧美| 最近最新中文字幕大全电影3| 欧美最新免费一区二区三区| 午夜a级毛片| 久久久久国内视频| 午夜精品一区二区三区免费看| 亚洲精品亚洲一区二区| 国产精品嫩草影院av在线观看| 不卡视频在线观看欧美| or卡值多少钱| 国内精品久久久久精免费| 少妇人妻一区二区三区视频| 三级毛片av免费| 国产精品一区二区三区四区免费观看 | 国产精品日韩av在线免费观看| 内地一区二区视频在线| 美女免费视频网站| 色视频www国产| 国内精品美女久久久久久| 在线观看免费视频日本深夜| av在线播放精品| 久久久欧美国产精品| 亚洲精品色激情综合| 女的被弄到高潮叫床怎么办| 国产爱豆传媒在线观看| 麻豆国产av国片精品| 国产高清不卡午夜福利| 亚洲精品影视一区二区三区av| 91在线观看av| 久久久久久久午夜电影| 人妻丰满熟妇av一区二区三区| 在线观看一区二区三区| 久久精品国产99精品国产亚洲性色| 国产免费一级a男人的天堂| 中国美女看黄片| 亚洲精华国产精华液的使用体验 | 亚洲内射少妇av| 又爽又黄无遮挡网站| 变态另类成人亚洲欧美熟女| 亚洲av熟女| 美女被艹到高潮喷水动态| 国产精品久久久久久精品电影| 干丝袜人妻中文字幕| 国产精品综合久久久久久久免费| 99在线视频只有这里精品首页| 黄色日韩在线| 亚洲七黄色美女视频| 亚洲最大成人中文| av免费在线看不卡| 别揉我奶头 嗯啊视频| 精品人妻视频免费看| 亚洲成人久久爱视频| 亚洲欧美精品自产自拍| 最近手机中文字幕大全| 亚洲中文字幕日韩| 六月丁香七月| 国产白丝娇喘喷水9色精品| 中国国产av一级| 欧美最黄视频在线播放免费| 变态另类丝袜制服| 99精品在免费线老司机午夜| 天堂√8在线中文| 国产aⅴ精品一区二区三区波| 国内揄拍国产精品人妻在线| 精品日产1卡2卡| 色在线成人网| 国产精品伦人一区二区| 午夜影院日韩av| 可以在线观看的亚洲视频| 日本-黄色视频高清免费观看| 中国美女看黄片| 亚洲色图av天堂| 日本爱情动作片www.在线观看 | 一个人观看的视频www高清免费观看| 成人永久免费在线观看视频| 日韩成人av中文字幕在线观看 | 男女视频在线观看网站免费| 婷婷亚洲欧美| 亚洲精品影视一区二区三区av| 黄色一级大片看看| 欧美激情在线99| 日韩国内少妇激情av| 99九九线精品视频在线观看视频| 国产成人一区二区在线| 一个人免费在线观看电影| 国产精品人妻久久久久久| 精品久久国产蜜桃| 久久久久久久久中文| 成年女人永久免费观看视频| 99久久精品一区二区三区| 亚洲自偷自拍三级| 精品一区二区三区av网在线观看| 亚洲av成人av| 日韩一区二区视频免费看| 亚洲精品色激情综合| 丝袜喷水一区| 黑人高潮一二区| 日韩欧美国产在线观看| 亚洲av免费在线观看| 久久人妻av系列| 日韩欧美一区二区三区在线观看| 国产aⅴ精品一区二区三区波| 一进一出好大好爽视频| 日韩欧美国产在线观看| 亚洲精品色激情综合| 丝袜喷水一区| 国产av不卡久久| 亚洲人成网站高清观看| 国产成人福利小说| 看片在线看免费视频| 亚洲国产欧美人成| 91在线观看av| 嫩草影院新地址| 国产伦一二天堂av在线观看| av在线播放精品| 黄色配什么色好看| 丝袜美腿在线中文| 在线观看66精品国产| 97热精品久久久久久| 欧美日韩综合久久久久久| 日韩欧美精品免费久久| 亚洲欧美日韩东京热| 美女免费视频网站| 我的老师免费观看完整版| 露出奶头的视频| 亚洲欧美精品综合久久99| 国产精品久久久久久久电影| 久久人人爽人人爽人人片va| 亚洲精品456在线播放app| 免费无遮挡裸体视频| 国产一区亚洲一区在线观看| 老司机午夜福利在线观看视频| 亚洲人成网站在线播放欧美日韩| 成人漫画全彩无遮挡| 国内久久婷婷六月综合欲色啪| 又粗又爽又猛毛片免费看| 狠狠狠狠99中文字幕| 一级黄片播放器| 国产av麻豆久久久久久久| 99久久精品热视频| 男女那种视频在线观看| 国产精品一区二区三区四区久久| 国产精品电影一区二区三区| a级毛片免费高清观看在线播放| 免费一级毛片在线播放高清视频| 亚洲无线在线观看| 搡老熟女国产l中国老女人| 乱人视频在线观看| 国内精品一区二区在线观看| 99久久无色码亚洲精品果冻| 女同久久另类99精品国产91| 国产片特级美女逼逼视频| 久久热精品热| 女人被狂操c到高潮| 18禁黄网站禁片免费观看直播| 亚洲av不卡在线观看| 久久这里只有精品中国| 国产在线精品亚洲第一网站| 国产黄a三级三级三级人| 波多野结衣巨乳人妻| 99九九线精品视频在线观看视频| 两性午夜刺激爽爽歪歪视频在线观看| 六月丁香七月| av在线蜜桃| 最近的中文字幕免费完整| 亚洲精品一区av在线观看| 干丝袜人妻中文字幕| 少妇丰满av| 尤物成人国产欧美一区二区三区| 亚洲自偷自拍三级| 18+在线观看网站| 黄色视频,在线免费观看| 免费看光身美女| 天天躁日日操中文字幕| 成人特级黄色片久久久久久久| 99久久精品国产国产毛片| 老熟妇仑乱视频hdxx| aaaaa片日本免费| 无遮挡黄片免费观看| 欧美极品一区二区三区四区| 国产精品国产三级国产av玫瑰| 亚洲人成网站高清观看| 淫秽高清视频在线观看| 久久鲁丝午夜福利片| 一个人观看的视频www高清免费观看| 在线观看免费视频日本深夜| 久99久视频精品免费| 最好的美女福利视频网| 亚洲av不卡在线观看| 老熟妇仑乱视频hdxx| 色哟哟·www| 人妻久久中文字幕网| 中文字幕久久专区| 91午夜精品亚洲一区二区三区| 亚洲五月天丁香| 国产伦一二天堂av在线观看| 成年女人永久免费观看视频| 国产久久久一区二区三区| 欧美激情久久久久久爽电影| 91午夜精品亚洲一区二区三区| 色哟哟·www| 国产在线男女| 亚洲国产欧美人成| 精品久久久久久久人妻蜜臀av| a级毛片a级免费在线| 人人妻人人看人人澡| 欧美精品国产亚洲| 最近最新中文字幕大全电影3| 久久这里只有精品中国| 色尼玛亚洲综合影院| 一进一出好大好爽视频| 成人亚洲精品av一区二区| 久久精品国产自在天天线| 男人狂女人下面高潮的视频| 村上凉子中文字幕在线| 欧美xxxx性猛交bbbb| 蜜臀久久99精品久久宅男| 国产精品伦人一区二区| 日本色播在线视频| 啦啦啦啦在线视频资源| 最好的美女福利视频网| 国产人妻一区二区三区在| 亚洲激情五月婷婷啪啪| 久久久久国内视频| 亚洲经典国产精华液单| videossex国产| 真人做人爱边吃奶动态| 亚洲美女视频黄频| 在线免费十八禁| 日韩国内少妇激情av| 亚洲av成人av| 校园春色视频在线观看| 国产精品无大码| 成人二区视频| 亚洲成人中文字幕在线播放| 国产精品久久久久久久久免| 日韩精品青青久久久久久| 午夜免费激情av| 国产老妇女一区| 六月丁香七月| 成人亚洲欧美一区二区av| 国产精华一区二区三区| 插逼视频在线观看| 岛国在线免费视频观看| 中国美白少妇内射xxxbb| 露出奶头的视频| 真人做人爱边吃奶动态| 免费不卡的大黄色大毛片视频在线观看 | 亚洲五月天丁香| 婷婷亚洲欧美| 老司机影院成人| 国产精品无大码| 老师上课跳d突然被开到最大视频| a级毛片免费高清观看在线播放| 日韩精品中文字幕看吧| 别揉我奶头~嗯~啊~动态视频| 亚洲乱码一区二区免费版| 亚洲图色成人| 精品久久国产蜜桃| 国产毛片a区久久久久| 老司机影院成人| 深夜a级毛片| 国产高清三级在线| 黄色视频,在线免费观看| 三级男女做爰猛烈吃奶摸视频| 亚洲人与动物交配视频| 久久这里只有精品中国| 小蜜桃在线观看免费完整版高清| 六月丁香七月| 欧美又色又爽又黄视频| 日本三级黄在线观看| 婷婷色综合大香蕉| 级片在线观看| 哪里可以看免费的av片| 国产成人freesex在线 | 国产三级在线视频| 久久久a久久爽久久v久久| 欧美性感艳星| 九九在线视频观看精品| 麻豆乱淫一区二区| www.色视频.com| 麻豆一二三区av精品| 性插视频无遮挡在线免费观看| 国产乱人视频| 久久精品国产鲁丝片午夜精品| 亚洲国产欧洲综合997久久,| 成人高潮视频无遮挡免费网站| 精品一区二区三区视频在线| 久久精品久久久久久噜噜老黄 | 中文字幕人妻熟人妻熟丝袜美| 国产一区二区三区在线臀色熟女| 在线观看免费视频日本深夜| 麻豆av噜噜一区二区三区| 99热只有精品国产| 国产伦在线观看视频一区| 一区二区三区免费毛片| 午夜福利在线在线| 色尼玛亚洲综合影院| 日韩强制内射视频| 看十八女毛片水多多多| 能在线免费观看的黄片| 少妇人妻一区二区三区视频| 校园春色视频在线观看| 成人特级av手机在线观看| 伦理电影大哥的女人| 久久久久九九精品影院| 精品人妻视频免费看| 成人综合一区亚洲| 校园人妻丝袜中文字幕| 国产精品综合久久久久久久免费| 亚洲国产精品久久男人天堂| 中文字幕熟女人妻在线| 在线观看免费视频日本深夜| 干丝袜人妻中文字幕| 免费高清视频大片| 国产白丝娇喘喷水9色精品| 在线免费观看不下载黄p国产| 久久久久久久久大av| 亚洲精华国产精华液的使用体验 | 精品久久久噜噜| 一夜夜www| 国产爱豆传媒在线观看| 国内久久婷婷六月综合欲色啪| 欧美性猛交╳xxx乱大交人| 精品一区二区三区视频在线观看免费| 成人二区视频| 国产一区二区激情短视频| 欧美最新免费一区二区三区| 麻豆av噜噜一区二区三区| 欧美精品国产亚洲| 国产黄片美女视频| 秋霞在线观看毛片| 亚洲精品色激情综合| 日韩成人伦理影院| 丝袜美腿在线中文| 国产熟女欧美一区二区| 国产淫片久久久久久久久| 午夜福利18| 淫妇啪啪啪对白视频| 两个人视频免费观看高清| 麻豆久久精品国产亚洲av| 少妇熟女欧美另类| 日本精品一区二区三区蜜桃| 免费无遮挡裸体视频| a级毛色黄片| 国产一区亚洲一区在线观看| 人人妻人人看人人澡| av视频在线观看入口| 中文亚洲av片在线观看爽| 俺也久久电影网| 欧美三级亚洲精品| 久久久久九九精品影院| 日本免费一区二区三区高清不卡| 国产单亲对白刺激| 国产一区二区三区av在线 | 国产乱人偷精品视频| 中文亚洲av片在线观看爽| 国产美女午夜福利| 亚洲精品日韩av片在线观看| 成年av动漫网址| 国产私拍福利视频在线观看| 亚洲av成人精品一区久久| 哪里可以看免费的av片| 日韩成人伦理影院| 日韩亚洲欧美综合|