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

潮汐應變對長江口北槽枯季湍流混合與層化的影響

2014-10-13 08:14:58熊龍兵胡國棟施慧燕
海洋工程 2014年4期
關鍵詞:差異

熊龍兵,浦 祥,時 鐘,胡國棟,施慧燕

(1.上海交通大學船舶海洋與建筑工程學院海洋工程國家重點實驗室,上海 200030;2.長江水利委員會水文局長江口水文水資源勘測局,上海 200136)

混合(mixing)與層化(stratification)是受淡水影響的河口海岸水域的基本物理特征,控制著垂向上的能量交換(如勢能、湍動能、熱能等)與物質輸移(如泥沙、營養鹽等),反之亦然,因此受到了海洋學家的廣泛關注。英國Fleming[1]在Tay estuary的觀測分析發現了河口鹽度層化現象。美國Pritchard[2]根據Chesapeake Bay的實測資料顯示,河口鹽淡水的混合在重力作用下往往會形成表層凈向海、底層凈向陸的縱向環流形式,即河口環流。美國Hansen and Rattray[3]研究河口鹽度平衡關系時采用了一個層化參數,即水體表、底層鹽度差與斷面平均鹽度的比值。英國Simpson and Hunter[4]、Simpson等[5]在研究陸架海鋒面時引入了層化指數(h/u3)并采用了一個定量反映層化強度的物理量,即勢能差異(potential energy anomaly)。在此基礎之上,Simpson等[5]導出了勢能差異方程(potential energy anomaly equation),并按照影響水體混合與層化的物理機制對其作了簡單分解。鑒于最初的勢能差異方程并未考慮河口地區徑流淡水注入的影響,Simpson等[6]對其作了進一步完善并考慮密度的水平對流作用,從而導出了河口環流及純潮流的勢能差異變化率公式,并提出潮汐應變(tidal straining)的概念。

所謂潮汐應變指的是潮流垂向流速剪切與水平密度梯度相互作用(如圖1),它傾向于使河口水體在一個潮周期內呈現混合與層化的交替變化:落潮時,水流向海,底摩擦效應導致垂向上形成明顯的流速剪切并使得密度等值線發生傾斜而產生層化(圖1(a));漲潮的過程與之相反,向陸水流在垂向上的反向流速剪切使得水體又趨于混合(圖1(b))。落潮期間,若潮汐應變造成的勢能差異的變化大于潮汐攪動,水體將呈現為應變致周期性層化(SIPS)狀態。Simpson等[6]利用這一概念解釋了Liverpool Bay混合與層化的潮相變化規律,此后,潮汐應變在許多河口的研究中得到應用,如:美國York River estuary(Sharples等[7]),德國Rhine ROFI(Simpson and Souza[8];Fisher 等[9]),美國 Hudson River estuary(Nepf and Geyer[10]),泉州灣(劉浩等[11]),黃河口(Wang 等[12]),長江口(李霞等[13])。

需要說明的是,英國Simpson等[6]所提出的潮汐應變概念是縱向一維的,適用于縱向密度梯度顯著、對流過程一維特征突出的河口,如Liverpool Bay。此后,一些學者基于觀測資料(Souza and Simpson[14];Lacy等[15])或數學模型(Scully等[16])研究發現,河口地區橫向上的流速剪切與密度梯度相互作用也會形成橫向潮汐應變,并對水體的混合與層化產生影響。在此基礎之上,德國Burchard and Hofmeister[17]和荷蘭de Boer等[18]采用海水密度的對流擴散方程分別導出了三維勢能差異方程,從而進一步完善了混合與層化的勢能差異理論。

圖1 湍流混合與潮汐應變概念圖 (基于Simpson等[6],Fig.2,p.128)Fig.1 Schematic of turbulent mixing and tidal straining(based on Simpson et al[6],Fig.2,p.128)

1 數學模型的設置與驗證

1.1 TELEMAC-3D簡介

TELEMAC-3D是由法國國家水力學實驗室(Laboratoire National d'Hydraulique)聯合歐洲多家科研機構共同開發的一個有限元水動力數學模型,并已廣泛應用于河流、湖泊、海岸河口等自然水域(Kopmann and Markofsky[35];Marques 等[36];Bedri等[37])。模型垂向上采用 σ 坐標變換,保證了計算網格在淺水區域依然具有較高的垂向分辨率。該坐標變換的基本思想源于大氣數值模擬(Phillips[38],p.184,paragraph 2,line 3),原始公式中σ取值由地面的1向上減小至0,而TELEMAC-3D中的σ的取值則由底床的0向上增至自由表面的 1,其坐標變換式為(Hervouet[39]):

式中:h為靜水深(m),η為自由表面水位(m)。

根據流體不可壓縮假定、Boussinesq假定和近似、靜壓假定即得到TELEMAC-3D的三維淺水方程,此外還包括鹽度對流擴散方程及海水狀態方程,其中,海水狀態方程為 (Hervouet[39],eq.(2.63),p.18):

式中:參考密度ρref=999.972 kg/m3,T為溫度,相應的參考溫度為Tref=4℃,S為鹽度。此方程的適用范圍:0℃ <T<40℃,0 psu<S<32 psu。這些基本控制方程需引入適當的湍流模型方可構成封閉系統,這里采用 k-ε 湍流模型,其 Cartesian 坐標系下的形式如下(Launder and Spalding[40];Rodi[41],eqs.(2.48-49),p.28):

式中:υ0,A0分別為分子渦動粘性系數和分子渦動擴散系數;且υ0=A0=1.0×10-6m2/s(Hearn[44],p.251)。

1.2 模型配置

文中模型范圍包括整個長江口,而主要研究區域則位于實測資料較全的北槽水域。就數學模擬的精度而言,為了盡量減小開邊界條件的誤差對北槽水域的干擾,模型開邊界應遠離該區域并綜合考慮潮流界與外海潮汐的影響。這里河流邊界設在江陰上游120°E,類似的文獻中也有報道(Xue等[45];Ge等[32]),外海的東邊界取在124°E,而南北邊界分別取在30.5°N和32.5°N,由此確定的模型計算區域如圖2所示。長江口北槽水域的地形資料來自長江水利委員會長江口水文水資源勘測局2010年第四季度的實測數據,而計算區域內的岸線及其余地形資料均來自“中國電子海圖2009(EMDAS Inc.版本:6920)”。北槽導堤的堤頂高程為吳淞基面以上2.0 m(陳志昌、樂嘉鉆[46]),經基準面換算后約為平均海面以上0.3 m(Ge等[32])。

水平方向上采用非結構三角網格離散,網格尺寸從外海開邊界上的約6 000 m平滑過渡至口內約400 m(如圖2)。此外,由于北槽航道兩側建有導堤和丁壩,為了更好地反映這些人工建筑物的影響,模型在北槽附近作了局部網格加密處理,其最小尺寸約為50 m。最終在水平方向上生成65 429個網格、33 762個節點,垂向上則以均勻分層的方式設置10個σ層。模型計算時間為2009年12月15日~2010年1月12日,共計28天。采用全隱式對流格式后,時間步長設為30 s。

圖2 長江口計算區域、水下地形及非結構三角網格Fig.2 Computational domain,bathymetry and unstructured triangular mesh within the Changjiang River estuary

1.3 邊界條件

河流邊界的徑流量采用大通水文站2009年12月~2010年1月實測徑流平均值12 000 m3/s,鹽度設為0 psu。由于東海前進潮波主要從東南方向傳入長江口(趙寶仁等[47]),故將東邊界和南邊界作為潮波入射邊界,設置 8 個主要分潮(M2、S2、N2、K2、K1、O1、P1、Q1),分潮調和常數由全球大洋模型 TPXO7.2 的中國海局部模型(http://volkov.oce.orst.edu/tides/YS.html)提供并按下式(Doodson[48],p.242)計算相應時刻的水位:

式中:N為分潮數目;fi為分潮的交點因子;Hi,ωi分別分潮的振幅和圓頻率;Vi,ui則分別為分潮的天文初相角和交點修正角;gi為分潮的Greenwich遲角。北部開邊界設為Thompson型輻射邊界(Thompson[49])。鹽度邊界條件直接由《渤海黃海東海海洋圖集(水文)》[50]多年平均的1月鹽度場插值得到。

自由表面考慮風應力作用,采用長江口南槽東站2010年1月的風速、風向觀測資料,流行風向為北偏西35°,平均風速約 5.5 m/s,自由表面風應力系數由下式計算(Hervouet[39],p.14):

式中:W即為相應風速。底部考慮摩擦作用,以對數流速分布的Nikuradse公式(Nikuradse[51];Hervouet[39],p.52)計算底部拖曳力系數:

式中:κ為Karman常數,一般取0.4;H為總水深;ks為底部粗糙度,它是定性反映底床粗糙程度的綜合參數,與床面泥沙粒徑、地形起伏變化、植被覆蓋情況等有關。CD的數值主要通過ks調節,該參數對水流的模擬結果有顯著影響。長江口地區底床表層沉積物類型主要以粘土質粉砂和砂為主,其中值粒徑在0.02~0.15 mm之間(劉紅[52],p.73),這里經調試率定后取ks=2.4×10-3m。

1.4 初始條件

水位、流速的響應過程迅速,因此均采用“冷啟動”方式。鹽度的響應過程較為緩慢,為了提高模型效率,文中的鹽度模擬采用“熱啟動”方式,即以《渤海 黃海 東海海洋圖集(水文)》實測的多年平均1月鹽度場插值得到平面上的初始鹽度場,并沿垂向各σ層復制。由于缺少溫度實測資料,模型將1月平均溫度統一設為6℃。

1.5 模型驗證

利用水利部長江水利委員會水文局長江口水文水資源勘測局2010年1月的水文實測資料對模型進行驗證,包括:橫沙、北槽中、牛皮礁三個潮位站的潮位過程;CSW、CS8兩個水文站大潮(2010年1月1~2日)和小潮(2010年1月9~10日)期間水深及表、中、底層流速、流向、鹽度過程。各驗證站在模型區域內的分布情況如圖3所示。

圖3 2010年枯季長江口北槽3潮位站(HS-橫沙,NPM-北槽中,NPJ-牛皮礁)及2水文觀測站(CSW,CS8)的位置Fig.3 Locations of three tidal stations(HS-Hengsha Station,NPM-North Passage(M)Station,NPJ-Niupijiao Station)and two hydrological gauging stations(CSW and CS8)within the North Passage of the Changjiang River estuary in the dry season of 2010

1.5.1 潮位驗證

模型對計算區域內的3個潮位站2009年12月30日~2010年1月11日的潮位過程進行了驗證(如圖4)。實測的潮位時間序列只有各潮周期的高潮位與低潮位數據。通過與模擬的連續潮位時間序列進行對比,3站潮位過程的模擬值與實測值總體吻合良好并且能夠較好地反映日潮不等現象,僅部分時段有所偏差,其中,大潮與中潮的計算偏差可能與底摩擦系數有關,而小潮的偏差除了與底摩擦系數有關外,還可能受到小潮期間風浪的影響。此外,模型對CSW及CS8站大潮和小潮水深時間序列的驗證(如圖5)也獲得了良好的精度。

圖4 長江口3潮位站2009年12月30日~2010年1月12日潮位的模擬值與實測值對比Fig.4 Comparisons between modeled and measured tidal elevations at three tidal stations within the Changjiang River estuary from 30 December 2009 to 12 January 2010

1.5.2 流速驗證

模型對北槽深水航道附近的CSW及CS8站大潮(2010年1月1日~2日)和小潮(2010年1月9日~10日)期間的分層流速、流向(表、中、底層)進行了驗證(如圖5)。CSW及CS8站分別位于北槽中段和北槽下段,實測的分層流速在漲潮期間隨水位的上升而逐漸增大,二者存在一定的相位差,落潮期間也有相似的規律,潮波形態并非嚴格的駐波。潮流流速基本表現為由表層向底層逐漸減小的規律。大潮至小潮期間,CSW站與CS8站的表層落潮流速明顯大于漲潮流速,而中層和底層的漲、落潮流強度相當,2站均呈現明顯的往復流特征。由圖5可知,CSW及CS8站模擬的各層流速、流向與實測值均吻合良好,模擬結果能夠較為合理地反映長江口北槽水流的基本特征。流速的偏差主要出現在落潮時刻,而流向在轉流時刻也存在一定偏差,這可能與北槽地形概化、口內底摩擦系數偏小等因素有關。

圖5 長江口北槽CSW站(左)及CS8站(右)大潮(2010年1月1~2日)和小潮(2010年1月9~10日)期間水深及表、中、底層流速大小與相應流向的模擬值與實測值對比Fig.5 Comparisons between modeled and measured water depth,surface,middle and bottom current speeds and their corresponding current directions at CSW(left panel)and CS8(right panel)within the North Passage of the Changjiang River estuary during spring tide(1st~ 2ndJanuary 2010)and neap tide(9th~ 10thJanuary 2010),respectively

1.5.3 鹽度驗證

模型對北槽深水航道附近的CSW及CS8站大潮(2010年1月1日~2日)和小潮(2010年1月9日~10日)期間的分層鹽度(表、中、底層)進行了驗證(如圖6)。CSW及CS8站的鹽度實測資料顯示,表、中、底層鹽度受潮汐影響呈明顯的周期性變化,具體表現為:漲潮期間,鹽度逐漸增大;落潮期間,鹽度逐漸減小。垂向上的鹽度分布基本呈現為由表層向底層逐漸增大的層化狀態或表層與底層基本相同的混合狀態。由圖6可知,大潮期間,CSW站與CS8站的表、中、底層鹽度的模擬值與實測值吻合相當好,偏差主要出現在CSW站的落潮時刻;小潮期間,CSW站表層和底層鹽度的模擬值與實測值有所偏差,且第二個潮周期的后半潮存在一定的相位差異,CS8站表層與中層鹽度的模擬值與實測值吻合良好,底層鹽度的模擬值在后半潮周期有所偏差但相位一致。鹽度模擬的偏差除了與水流模擬精度有關之外,也受初始鹽度場與鹽度邊界條件的影響。

1.5.4 北槽垂向鹽度場與河口環流

小潮期間,從北槽向陸端至北槽向海端的潮平均鹽度與水流沿航道的分布如圖7所示。此時由鹽度造成的層化較為明顯并貫穿于整個北槽航道,并且彎道附近的層化強度更高。北槽水體垂向上出現的明顯層化在重力作用下形成了表層凈向海、底層凈向陸的縱向環流形式,即河口環流(Pritchard[2],Fig.5,p.251),它對沿航道方向的輸移過程有著重要影響。表層受徑流下泄的影響,環流強度可達0.6 m/s,而下層平均約為0.2 m/s。大潮期間,由于水體垂向上的層化總體較弱,斜壓效應并不明顯,因此不會出現明顯的河口環流。

圖6 長江口北槽CSW站(左)、CS8站(右)大潮(2010年1月1~2日)和小潮(2010年1月9~10日)期間表、中、底層鹽度的模擬值與實測值對比Fig.6 Comparisons between modeled and measured surface,middle and bottom salinities at CSW(left panel)and CS8(right panel)within the North Passage of the Changjiang River estuary during spring tide(1st~2ndJanuary 2010)and neap tide(9th~10thJanuary 2010),respectively

圖7 模擬的長江口北槽2010年枯季小潮期間潮平均的流速矢量與鹽度縱向分布Fig.7 Longitudinal distributions of modeled tidally-averaged current vectors and salinity within the North Passage of the Changjiang River estuary during neap tide in the dry season of 2010

1.5.5 北槽平面流場與鹽度場

小潮期間漲急與落急時刻的表層流場與鹽度場如圖8所示。漲急時刻(圖8(a)),北槽導堤和丁壩均被淹沒,堤頂處存在非常明顯的越堤流,流速可達2.0 m/s以上且流向偏北。主航槽的流速顯著大于壩田區及航道外側,丁壩的束流作用明顯。受水流影響,主航槽鹽水入侵更為明顯而壩田區鹽度較低,高鹽水的上溯范圍較大;落急時刻(圖8(b)),北槽導堤和丁壩露出水面,水流歸槽,主航槽下泄流速較大而壩田區流速明顯減小,且存在尺度不一的平面環流,受其影響,主航槽與壩田區的鹽度分布趨于均勻,而高鹽水的上溯范圍也已明顯減小。此結果表明,模型已在一定程度上體現了導堤和丁壩對水流、鹽度的影響。

圖8 模擬的長江口北槽附近2010年枯季小潮期間漲急時刻和落急時刻的表層流速矢量與鹽度水平分布Fig.8 Horizontal distributions of modeled surface current vectors and salinities near the North Passage of the Changjiang River estuary(a)at maximum flood and(b)maximum ebb during neap tide in the dry season of 2010,respectively

2 混合與層化的機制分析

上一節對長江口三維水動力數學模型的驗證與分析顯示,此模型水流和鹽度的模擬均獲得了良好的精度,模擬結果合理可靠。在此基礎上,本節將利用模擬得到的潮位、流速和鹽度,根據勢能差異理論進行相關的定量計算,并結合湍動能耗散率的分布對北槽下段CS8站及毗鄰的北槽水域枯季混合與層化的特征及其物理機制展開分析和探討。

2.1 計算公式

2.1.1 勢能差異公式

河口水體垂向上混合與層化的變化過程,從能量角度而言,體現的是勢能與動能相互轉化的過程,而勢能差異理論則是研究這一重要物理過程的有效途徑。基于英國Simpson and Hunter[4]和Simpson等[5],Simpson[34](p.23)提出了水體“勢能差異”的計算公式:

2.1.2 勢能差異變化率公式

當不考慮表面熱交換、降雨等因素的影響時,由式(10)即可得到水體的勢能差異變化率 (Simpson等[6],p.127):

式中積分上限也已調整至自由表面η,此即為廣泛應用的勢能差異方程。由于水域特征與動力機制的差異,實際應用中,不同學者針對勢能差異方程作了不同程度的簡化。英國Simpson et al[6]根據該方程,假定密度的變化僅受x向對流過程的影響,且水平密度梯度沿垂向不變,即得到經典的潮汐應變公式(Simpson et al[6],eq.(1),p.127;de Boer et al[18],eq.(3),p.2):

在河口地區,影響水體垂向上混合與層化的另一重要物理機制為潮汐、風等外部動力的“攪動”所造成的湍流混合作用。與潮汐應變不同,攪動作用起著持續減弱水體層化的作用。德國Burchard and Hofmeister[17]通過推導三維勢能差異方程給出了潮汐與風共同攪動引起的勢能差異變化率的理論計算式(Burchard and Hofmeister[17],eq.(14),p.681):

式中:B即為浮力生成率,已根據k-ε湍流模型的形式調整了它的正負號。此公式中的垂向渦動擴散系數Az是潮流與風共同作用的結果,因此,式(13)表示的是潮汐與風共同攪動的作用。

由于北槽航道及導堤、丁壩的存在,模擬結果顯示,沿航道流速遠大于跨航道流速,且沿航道的鹽度梯度明顯占優,一維特征明顯。故僅考慮沿航道方向的潮汐應變和潮汐攪動作用,由式(12)~(13)得到本文的勢能差異方程(Burchard and Hofmeister[17],p.684):

式中的x方向即為沿航道方向,根據此公式可計算水體總的勢能差異變化率。

2.1.3 Simpson數

為了了解水體混合與層化的狀態,可以采用無量綱的Simpson數(Si),該參數由美國Monismith等[53](p.147)提出,也被稱為水平 Richardson數(Geyer and Ralston[54],p.39),其基本形式為:

2.1.4 梯度Richardson數

采用無量綱的梯度Richardson數對水體的穩定性作出判斷,該參數表征了水體勢能與動能之比,其計算公式為 (Holzman[55],eq.(1),p.13):

2.2 北槽水域混合與層化的大、小潮變化特征

本節根據式(10)計算了北槽及鄰近水域大潮平均與小潮平均的勢能差異空間分布,以初步認識該水域混合與層化的基本特征,如圖9所示。

圖9 模擬的長江口北槽附近2010年枯季大潮和小潮的潮平均勢能差異分布Fig.9 Distributions of modeled tidal mean potential energy anomalies near the North Passage of the Changjiang River estuary over(a)a spring tide and(b)a neap tide in the dry season of 2010,respectively

大潮期間(圖9(a)),長江口北槽水域總體層化較弱,北槽中段的勢能差異(約30 J/m3)略高于北槽上段和下段(約10 J/m3),且較強的層化主要出現在主航槽。壩田區水體的勢能差異約0 J/m3,幾乎為完全混合狀態;小潮期間(圖9(b)),長江口北槽水域的勢能差異不僅明顯高于大潮,高勢能差異的水域較大潮也有明顯擴大,其空間范圍從橫沙水道附一直延伸至北槽口門以外的海域。北槽中段的勢能差異(約90 J/m3)依然略高于兩端(約50 J/m3),且較強的層化仍然出現在主航槽附近。壩田區水體的勢能差異接近0 J/m3,混合依然較好。

由此可知,枯季長江口北槽水域的層化特征隨時間和空間有明顯變化:小潮期間的層化范圍和強度均明顯超過大潮;主航槽水域的層化強度始終高于壩田區,且北槽彎道附近的層化強度也高于北槽兩端。

2.3 北槽下段混合與層化的物理機制

為了初步了解混合與層化的物理機制的空間分布,分別選取大潮與小潮期間的某一典型落急時刻,根據式(12)及(13)計算北槽下段潮汐應變和潮汐與風共同攪動引起的勢能差異變化率(如圖10)。

大潮期間,潮汐應變在主航槽靠近丁壩的水域作用較強,其勢能差異變化率約為10×10-4W/m3,并且在彎道和口門附近的水域有明顯的局部極大值(約100×10-4W/m3),而壩田區的潮汐應變作用很弱(約-10×10-4W/m3)。潮汐與風共同攪動在靠近北槽彎道的主航槽水域作用較弱(約10×10-4W/m3),強攪動主要發生壩田區外側和丁壩附近的淺水區。向下靠近口門的主航槽水域,攪動作用明顯增強,其勢能差異變化率可達80×10-4W/m3,而相鄰壩田區的攪動作用則迅速減弱;小潮期間,潮汐應變在整個主航槽水域都很強,且堤頂附近的極大值特征十分明顯(約100×10-4W/m3),這可能與丁壩頂端的局部流場、鹽度場有關,壩田區潮汐應變引起的勢能差異變化率略小于0 W/m3。潮汐與風共同攪動在從彎道至口門的整個北槽下段主航槽水域已顯著減弱(約10×10-4W/m3),而壩田區外側和丁壩附近的淺水區仍然有較強的攪動作用。

由此可知,就北槽下段水域某一落急時刻而言,潮汐應變和潮汐與風共同攪動的強度在主航槽與壩田區有非常顯著的差異,受導堤和丁壩的影響明顯。不同的位置,導致混合與層化的主要物理機制往往不同;即使同一位置,大潮與小潮期間導致混合與層化的主要物理機制也有變化。

圖10 模擬的長江口北槽下段2010年枯季大潮(上)和小潮(下)落急時刻由潮汐應變(左)及潮汐與風共同攪動(右)引起的勢能差異變化率分布Fig.10 Distributions of modeled time derivative of the potential energy anomaly contributed by tidal straining(left panel)and combined tidal and wind stirring(right panel)within the lower reach of North Passage of the Changjiang River estuary at maximum ebb during spring tide(upper panel)and neap tide(lower panel)in the dry season of 2010,respectively

2.4 北槽下段CS8站混合與層化的特征和機制

北槽下段CS8站2010年1月1日~11日水深、水深平均的沿航道流速、勢能差異、潮汐應變及潮汐與風共同攪動引起的勢能差異變化率、總的勢能差異變化率以及潮平均Si數等參數的時間序列如圖11(a)~(f),包括了大潮(2010年1月1~2日)、中潮(2010年1月5~6日)和小潮(2010年1月9~10日)三個時段。

對比圖11(a)與11(b)可知,潮波形態在大潮至中潮介于行波與駐波之間,而小潮的駐波特征明顯。由圖11(c)可知,CS8站大潮至中潮的勢能差異在0 J/m3與30 J/m3之間波動,水體呈完全混合與弱層化的交替變化,平均勢能差異約10 J/m3。小潮的勢能差異始終大于30 J/m3,水體一直處于層化狀態,勢能差異平均值增至約60 J/m3。此外,勢能差異還存在明顯的潮周期波動。

圖11(d)顯示,潮汐應變引起的勢能差異變化率從大潮至小潮呈總體上升的趨勢,具有非常明顯的潮周期變化。其中,大潮至中潮,潮汐應變引起的勢能差異變化率在漲、落潮的平均值分別約為 -20×10-4、10×10-4W/m3。小潮期間由于斜壓效應引起的河口環流增強,潮汐應變已明顯增大,其勢能差異變化率漲、落潮的平均值分別約為 -10×10-4、30×10-4W/m3。潮汐與風共同攪動引起的湍流混合作用也具有明顯的潮周期變化,其勢能差異變化率從大潮至小潮是逐漸減小的,平均值分別約為30×10-4、10×10-4W/m3。圖11(e)中由潮汐應變和潮汐與風共同攪動引起的總的勢能差異變化率與勢能差異的時間變化過程(圖11(c))基本對應,其在大潮至中潮主要受潮汐與風共同攪動的控制,而小潮主要受潮汐應變的控制。

潮平均的Si數時間序列如圖11(f)。大潮至中潮,Si數略大于0.1且較為穩定,而小潮的Si數明顯增大并超過1.0。Si數逐漸增大反映了潮汐應變相對于潮汐與風共同攪動的湍流混合作用逐漸增強的規律。德國Becherer等[56]對Simpson等[6]的研究成果作了進一步分析,并總結得到Si數的上、下臨界值分別為0.84和0.088:當0.088<Si<0.84時,屬于SIPS;當Si<0.088時,屬于完全混合;當Si>0.84時,屬于持續性層化。Si數上、下臨界值的量級約為1.0和0.1,據此可以認為,枯季北槽下段CS8站的層化狀態在大潮至中潮屬于SIPS,而小潮則是明顯的持續性層化。

圖11 長江口北槽水文觀測站CS8模擬和計算的大潮至小潮(2010年1月1~11日)各種物理參數的時間序列Fig.11 Time series of various physical parameters from spring tide to neap tide(1st~11thJanuary 2010)modeled and calculated at the hydrological station CS8 within the North Passage of the Changjiang River estuary,respectively

需要說明的是,大潮至中潮期間,會出現漲潮層化增強而落潮層化減弱的情況,這與Simpson等[6]的經典潮汐應變理論有所差異,除了與潮波形態 (Fisher等[9])、轉流時刻表底層水流流向差異有關之外,另一個重要因素可能是橫向潮汐應變。美國Lacy等[15]認為,在混合作用較為明顯的大潮期間,橫向潮汐應變的影響比較顯著。美國Scully等[16]通過數學模擬研究Hudson River estuary發現,大潮期間,由于橫向潮汐應變的影響,漲潮時的層化反而強于落潮。由于北槽水域缺乏橫向驗證資料,因此這里并未予以考慮。

2.5 北槽下段CS8站潮汐應變對湍動能耗散率的影響

上一節的分析表明,潮汐應變是導致北槽下段CS8站水體層化強度出現漲、落潮不對稱分布的主要物理機制。在此基礎之上,本節將根據該站大潮(2010年1月1日~2日)和小潮(2010年1月9日~10日)期間的水流un、湍動能耗散率ε及梯度Ri數時間序列(圖12(a)~(c))繼續探討層化規律與湍動能耗散率之間可能存在的關系,并分析潮汐應變的影響。

沿航道方向的垂向流速分布如圖12(a)所示。圖12(b)顯示了CS8站大潮和小潮期間的湍動能耗散率分布特征,二者隨時間都有明顯的M4周期性特征。大潮期間,底部以上0~9 m范圍內,湍動能耗散率在漲急和落急時刻有極大值(平均約10-4W/kg),而在轉流時刻有極小值(約10-7W/kg),漲、落潮不對稱性明顯。空間上,湍動能耗散率一般底部最大(約10-3W/kg),這主要是由底摩擦作用所導致的強水流剪切引起的。沿底部向上,隨著水流剪切逐漸減弱,湍動能耗散率也趨于減小,而中間水層的湍動能耗散率在漲潮和落潮后期減小非常明顯(約10-9W/kg)。接近自由表面區域的湍動能耗散率較相鄰的下層水體反而有所增大,這是由于文中在自由表面考慮了定常風的影響,風應力增強了表層的湍流混合作用。小潮期間,CS8站湍動能耗散率的分布規律與大潮類似,但中間水層的低耗散率區較大潮有明顯擴大。

圖12 長江口北槽水文觀測站CS8模擬和計算的大潮(2010年1月1~2日)和小潮(2010年1月9~10日)物理參數的垂向時間序列,湍動能耗散率ε由k-ε湍流模型估算得到Fig.12 Time series of vertical profiles of physical parameters during spring tide(1st~2ndJanuary 2010)and neap tide(9th~10thJanuary 2010)modeled and calculated at the hydrological gauging station CS8 within the North Passage of the Changjiang River estuary(TKE dissipation rate ε is estimated from the k-ε turbulence model)

圖12(c)為大潮和小潮期間水體垂向上梯度 Ri數的分布,其臨界值取為0.25(Taylor[57];Miles[58])。由圖可知,CS8站梯度Ri數的量級大潮約為10-3~100,小潮約為10-2~101,而 Wang等[25]和倪智慧等[26]估算的長江口外羽狀流水域梯度Ri數的量級約10-3~101,個別時刻可達102,與這里的計算結果基本接近。通過對比CS8站梯度Ri數(圖12(c))與湍動能耗散率(圖12(b))的分布,大潮和小潮期間,該站中間水層出現的低耗散區與穩定層化區(Ri>0.25)存在明顯的對應關系,即層化越強的區域,湍動能耗散率越小。這可用湍動能輸運方程來解釋,為簡便起見忽略湍動能對流擴散項,則式(3)變為如下形式(Stacey等[59],p.19):

各項定義見1.1節的k-ε湍流封閉模型。剪切生成率P起著持續輸入湍動能并增強湍流耗散的作用,浮力生成率B會隨水體層化狀態的變化而起到完全不同的作用:當出現穩定層化時(Ri>0.25),B消耗湍動能并將其轉化為勢能,從而減小ε;當出現不穩定層化或混合時(Ri<0.25),B提供湍動能并將勢能轉化為動能,從而增大ε,由此揭示了強層化區出現低耗散的原因。需要注意的是,河口水體懸沙濃度所導致的密度層化也會對垂向上的湍流混合過程產生影響(Winterwerp[43])。根據式(17),懸沙導致的層化使得浮力生成率B增大,從而使得底床附近的湍動能耗散率ε減小。CS8站2010年枯季大潮和小潮實測的懸沙濃度范圍分別為0.5~1.2、0.2~0.8 g/L,由懸沙導致的密度層化很弱且主要位于底層水體。定量分析結果顯示,CS8站近底層的剪切作用明顯占優(考慮懸沙后的底部浮力生成率B仍小于剪切生成率P約1個量級),懸沙所導致的層化對湍動能耗散率的分布(圖12(b))并無明顯影響,這里謹慎忽略之。

湍動能耗散率漲、落潮的不對稱分布和水體混合與層化的狀態有明顯聯系,并內在地受到潮汐應變的控制:根據圖12(b)并對比圖11(c)~(d)可知,大潮期間,由于潮汐應變使得漲潮期間層化增強而落潮期間層化減弱,從而導致強耗散在漲潮時受到抑制難以向上層擴散,而在落潮時,強耗散卻幾乎充斥整個垂向水體。小潮期間,潮汐應變使得漲潮時層化減弱而落潮時層化增強,從而導致強耗散在漲潮時向上層擴散范圍較大而落潮時向上層擴散受到抑制,且由于潮汐應變明顯增強,使得整個潮周期內都呈一定的層化狀態,強耗散一直被限制在中下層水體而難以擴散至表層。可見潮汐應變控制著水體層化強度在漲、落潮的變化,是導致CS8站湍動能耗散率漲、落潮不對稱分布的重要因素。

3 結語

1)北槽水域2010年枯季大潮平均的勢能差異范圍約0~30 J/m3,小潮平均的勢能差異范圍約0~90 J/m3,且較大的勢能差異始終位于主航槽,壩田區則一直接近0 J/m3。該水域混合與層化的特征隨時間、空間有明顯差異,具體表現為:時間上,小潮期間的層化明顯強于大潮;空間上,主航槽的層化始終強于壩田區,而北槽中段往往具有更強的層化。

2)落急時刻,就北槽下段而言,潮汐應變、潮汐與風共同攪動引起的勢能差異變化率范圍分別約為 -20×10-4~100×10-4W/m3、0~100×10-4W/m3。從大潮至小潮,潮汐應變總體增強而潮汐與風共同攪動總體減弱。空間上,潮汐應變、潮汐與風共同攪動的強度在壩田區內、外也存在較大差異,受導堤和丁壩的影響明顯。不同的位置,影響混合與層化的主要物理機制往往不同;即使同一位置,大潮與小潮期間影響混合與層化的主要物理機制也有變化。

3)北槽下段CS8站大潮至中潮的平均勢能差異約10 J/m3,Si數范圍為0.15~0.4,小潮的平均勢能差異增至約60 J/m3,相應的Si數范圍為0.9~1.5。該站混合與層化的物理機制受到的潮汐影響,在大潮和中潮期間,潮汐與風共同攪動造成的湍流混合作用占優,但落潮后期的潮汐應變依然較大,水體呈周期性的混合與弱層化交替變化,屬SIPS狀態;小潮期間,斜壓效應導致的河口環流增強使得潮汐應變占優,水體一直呈較為明顯的層化,屬持續性層化狀態。潮汐應變是控制該站水體層化強度漲、落潮不對稱分布的主要物理機制。

4)北槽下段CS8站梯度Ri數的量級范圍在混合較好的表、底層約10-3~10-2,而在層化較好的中間水層約100~101。湍動能耗散率的量級范圍大潮為10-3~10-9W/kg,小潮為10-5~10-10W/kg,其分布隨時間有明顯的M4周期性特征,表、底層分別由于風應力和底摩擦作用而具有較大的耗散率,中間水層出現的低耗散區則與該區域出現的穩定層化直接相關。湍動能耗散率的漲、落潮不對稱分布特征與水體層化狀態直接相關,并受潮汐應變的控制。大潮期間,潮汐應變使得漲潮層化增強而落潮層化減弱,從而導致強耗散在漲潮期間向上層擴散受到抑制而落潮期間向上層擴散較為充分;小潮期間漲、落潮的規律類似但相位相反。

[1] Fleming J.Observations on the junction of the fresh waters of rivers with the salt water of the sea[J].Transactions of the Royal Society of Edinburgh,1816,8:507-513.

[2] Pritchard D W.Estuarine Hydrography[M].In:Landsberg H E(Ed.),Advances in Geophysics,Academic Press,1952,1:243-280.

[3] Hansen D V,Rattray M J.New dimensions in estuary classification[J].Limnology and Oceanography,1966,11(3):319-326.

[4] Simpson J H,Hunter J R.Fronts in the Irish Sea[J].Nature,1974,250:404-406.

[5] Simpson J H,Allen C M,Morris N C G.Fronts on the continental shelf[J].Journal of Geophysical Research,1978,83(C9):4607-4614.

[6] Simpson J H,Brown J,Matthews J,et al.Tidal straining,density currents,and stirring in the control of estuarine stratification[J].Estuaries,1990,13(2):125-132.

[7] Sharp1es J,Simpson J H,Brubaker J M.Observations and modelling of periodic stratification in the Upper York River Estuary,Virginia[J].Estuarine,Coastal and Shelf Science,1994,38(3):301-312.

[8] Simpson J H,Souza A J.Semidiurnal switching of stratification in the region of freshwater influence of the Rhine[J].Journal of Geophysical Research,1995,100(C4):7037-7044.

[9] Fisher N R,Simpson J H,Howarth M J.Turbulent dissipation in the Rhine ROFI forced by tidal flow and wind stress[J].Journal of Sea Research,2002,48(4):249-258.

[10] Nepf H M,Geyer W R.Intratidal variations in stratification and mixing in the Hudson estuary[J].Journal of Geophysical Research,1996,101(C5):12079-12086.

[11]劉 浩,駱智斌,潘偉然.泉州灣水體結構的潮周期變化[J].臺灣海峽,2009,28(3):316-320.(LIU Hao,LUO Zhibin,PAN Wei-ran.Intratidal variation of the water column in Quanzhou Bay[J].Journal of Oceanography in Taiwan Strait,2009,28(3):316-320.(in Chinese))

[12] Wang X H,Wang H J.Tidal straining effect on the suspended sediment transport in the Huanghe(Yellow River)Estuary,China[J].Ocean Dynamics,2010,60(5):1273-1283.

[13]李 霞,胡國棟,時 鐘,等.長江口南支南港北槽枯季水體中混合、層化與潮汐應變[J].水運工程,2013(9):79-88.(LI Xia,HU Guo-dong,SHI Zhong,et al.Mixing,stratification and tidal straining in dry season within the North Passage of the Changjiang River estuary,South Branch,South Channel[J].Port& Waterway Engineering,2013(9):79-88.(in Chinese))

[14] Souza A J,Simpson J H.Controls on stratification in the Rhine ROFI system[J].Journal of Marine Systems,1997,12(1):311-323.

[15] Lacy J R,Stacey M T,Burau J R,et al.Interaction of lateral baroclinic forcing and turbulence in an estuary[J].Journal of Geophysical Research,2003,108(C3),3089,doi:10.1029/2002JC001392.

[16]Scully M E,Geyer W R,Lerczak J A.The influence of lateral advection on the residual estuarine circulation:a numerical modeling study of the Hudson River estuary[J].Journal of Physical Oceanography,2009,39:107-124.

[17]Burchard H,Hofmeister R.A dynamic equation for the potential energy anomaly for analysing mixing and stratification in estuaries and coastal seas[J].Estuarine,Coastal and Shelf Science,2008,77(4):679-687.

[18] de Boer G J,Pietrzak J D,Winterwerp J C.Using the potential energy anomaly equation to investigate tidal straining and advection of stratification in a region of freshwater influence[J].Ocean Modelling,2008,22(1):1-11.

[19] Rippeth T P,Fisher N R,Simpson J H.The cycle of turbulent dissipation in the presence of tidal straining[J].Journal of Physical Oceanography,2001,31(8):2458-2471.

[20]Simpson J H,Burchard H,Fisher N R,et al.The semi-diurnal cycle of dissipation in a ROFI:model-measurement comparisons[J].Continental Shelf Research,2002,22(11):1615-1628.

[21]Fischer E,Burchard H,Hetland R D.Numerical investigations of the turbulent kinetic energy dissipation rate in the Rhine region of freshwater influence[J].Ocean Dynamics,2009,59(5):629-641.

[22]毛漢禮,甘子鈞,藍淑芳.長江沖淡水及其混合問題的初步探討[J].海洋與湖沼,1963,5(3):183-206.(MAO Hanli,GAN Zi-jun,LAN Shu-fang.A preliminary study of the Yangtze diluted water and its mixing processes[J].Oceanologia et Limnologia Sinica,1963,5(3):183-206.(in Chinese))

[23]張重樂,沈煥庭.長江口咸淡水混合及其對懸沙的影響[J].華東師范大學學報:自然科學版,1988,4:83-88.(ZHANG Chong-le,SHEN Huan-ting.Mixing of salt-and fresh-water in the Changjiang estuary and its effects on suspended sediment[J].Journal of East China Normal University,Natural Science,1988,4:83-88.(in Chinese))

[24] Shen H T,Zhang C L.Mixing of salt water and fresh water in the Changjiang River estuary and its effects on suspended sediment[J].Chinese Geographical Science,1992,2(4):373-381.

[25] Wang K S,Cheng H,Dong L X.A hydrographic comparison of the two sides of the Changjiang plume fronts[C]//Proceedings of International Symposium on Biogeochemical Study of the Changjiang Estuary and its Adjacent Coastal Waters of the East China Sea.Beijing:China Ocean Press,1990:62-75.

[26]倪智慧,陳 輝,董禮先,等.長江口外羽狀流水體中的垂向混合與層化的觀測與分析[J].上海交通大學學報,2012,46(11):1862-1873.(NI Zhi-hui,CHEN Hui,DONG Li-xian,et al.Measurement and analysis of vertical mixing and stratification within the plume outside the Changjiang River estuary[J].Journal of Shanghai Jiao Tong University,2012,46(11):1862-1873.(in Chinese))

[27]匡翠萍.長江口鹽水入侵三維數值模擬[J].河海大學學報,1997,4(25):54-60.(KUANG Cui-ping.A 3-D numerical model for saltwater intrusion in the Changjiang Estuary[J].Journal of Hohai University,1997,4(25):54-60.(in Chinese))

[28]朱建榮,肖成猷,沈煥庭.夏季長江口沖淡水擴展的數值模擬[J].海洋學報,1998,20(5):13-22.(ZHU Jian-rong,XIAO Cheng-you,SHEN Huan-ting.Numerical model simulation of expansion of Changjiang diluted water in summer[J].Acta Oceanologica Sinica,1998,20(5):13-22.(in Chinese))

[29]周濟福,劉青泉,李家春.河口混合過程的研究[J].中國科學:A輯,1999,29(9):835-843.(ZHOU Ji-fu,LIU Qingquan,LI Jia-chun.Study on mixing processes in estuaries[J].Science in China,Series A,1999,29(9):835-843.(in Chinese))

[30]龔 政.長江口三維斜壓流場及鹽度場數值模擬[D].南京:河海大學,2002:171.(GONG Zheng.Three-dimensional baroclinic numerical model of current and salinity in Yangtze Estuary[D].Nanjing:Hohai University,2002:171.(in Chinese))

[31]羅小峰,陳志昌.長江口北槽近期鹽度變化分析[J].水運工程,2006(11):79-82.(LUO Xiao-feng,CHEN Zhi-chang.Analysis of salinity’s variation in the North Channel of the Yangtze Estuary[J].Port& Waterway Engineering,2006(11):79-82.(in Chinese))

[32] Ge J,Ding P,Chen C.Impacts of deep waterway project on local circulations and salinity in the Changjiang estuary,China[C]//Proceedings of the 32ndInternational Conference on Coastal Engineering.2011.

[33] Shi J Z,Lu L F.A short note on the dispersion,mixing,stratification and circulation within the plume of the partially-mixed Changjiang River estuary,China[J].Journal of Hydro-environment Research,2010,5(2):111-126.

[34] Simpson J H.The shelf-sea fronts:implications of their existence and behavior[J].Philosophical Transactions of the Royal Society of London,1981,A302:531-546.

[35] Kopmann R,Markofsky M.Three-dimensional water quality modeling with TELEMAC-3D[J].Hydrological Processes,2000,14:2279-2292.

[36] Marques W C,Fernandes E H L,Moller O O.Straining and advection contributions to the mixing process of the Patos Lagoon coastal plume,Brazil[J].Journal of Geophysical Research,2010,115(C6):1-23.

[37] Bedri Z,Bruen M,Dowley A,et al.A three-dimensional hydro-environmental model of Dublin Bay[J].Environmental Modeling& Assessment,2011,16(4):369-384.

[38] Phillips N A.A coordinate system having some special advantages for numerical forecasting[J].Journal of Meteorology,1957,14:184-185.

[39] Hervouet J M.Hydrodynamics of Free Surface Flows:Modeling with the Finite Element Method[M].John Wiley & Sons,Ltd.,2007:341.

[40] Launder B E,Spalding D B.The numerical computation of turbulent flows[J].Computer Methods in Applied Mechanics and Engineering,1974,3(2):269-289.

[41] Rodi W.Turbulence Models and Their Application in Hydraulics:A State-Of-The-Art Review[M].2ndedition,Taylor and Francis,1984:104.

[42] Uittenbogaard R E.The importance of internal waves for mixing in a stratified estuarine tidal flow[D].Delft:Delft University of Technology,1995:328.

[43] Winterwerp J C.Stratification effects by cohesive and noncohesive sediment[J].Journal of Geophysical Research,2001,106(C10):22559-22574.

[44] Hearn C.The Dynamics of Coastal Models[M].New York:Cambridge University Press,2008:488.

[45] Xue P,Chen C,Ding P,et al.Saltwater intrusion into the Changjiang River:a model-guided mechanism study[J].Journal of Geophysical Research,2009,114(C02006),doi:10.1029/2008JC004831.

[46]陳志昌,樂嘉鉆.長江口深水航道整治原理[J].水利水運工程學報,2005(1):1-7.(CHEN Zhi-chang,LE Jia-zuan.Regulation principle of the Yangtze River Estuary deep channel[J].Hydro-Science and Engineering,2005(1):1-7.(in Chinese))

[47]趙保仁,方國洪,曹德明.渤、黃、東海潮汐潮流的數值模擬[J].海洋學報,1994,16(5):1-10.(ZHAO Bao-ren,FANG Guo-hong,CAO De-ming.Numerical modeling of tides and currents in Bohai Sea,Yellow Sea,and East China Sea[J].Acta Oceanologica Sinica,1994,16(5):1-10.(in Chinese))

[48] Doodson A T.The analysis of tidal observations[J].Philosophical Transactions of the Royal Society of London,1928,A227:223-279.

[49] Thompson K W.Time-dependent boundary conditions for hyperbolic systems,II[J].Journal of Computational Physics,1990,89(2):439-461.

[50]海洋圖集委員會.渤海 黃海 東海海洋圖集(水文)[M].北京:海洋出版社,1992:524.(Committee of Marine Atlas.Marine Atlas of Bohai Sea,Yellow Sea,East China Sea(Hydrology)[M].Beijing:Ocean Press,1992:524.(in Chinese))

[51] Nikuradse J.Stromungsgesetze in Rauhen Rohren[R].Forschung auf dem Gebiete des Ingenieurwesens,1933:No.361(English translation,1950,Laws of flow in rough pipes.NACA TM 1292,Washington,D.C.,62 pp.).

[52]劉 紅.長江河口泥沙混合和交換過程研究[D].上海:華東師范大學,2009:119.(LIU Hong.Sediment mixing and exchange processes in the Yangtze Estuary[D].Shanghai:East China Normal University,2009:119.(in Chinese))

[53]Monismith S G,Burau J R,Stacey M T.Stratification dynamics and gravitational circulation in northern San Francisco Bay[M].San Francisco Bay:The Ecosystem,AAAS,1996:123-153.

[54] Geyer W R,Ralston D K.The dynamics of strongly stratified estuaries[M].In:Wolanski E and McLusky D(Eds.),Treatise on Estuarine and Coastal Science,2011:37-51.

[55] Holzman B.The influence of stability on evaporation[J].Annals of the New York Academy of Sciences,1943,44(1):13-18.

[56] Becherer J,Burchard H,Fl?ser G,et al.Evidence of tidal straining in well-mixed channel flow from micro-structure observations[J].Geophysical Research Letters,2011,38(17):L 17611.

[57] Taylor G I.Effect of variation in density on the stability of superposed streams of fluid[J].Proceedings of the Royal Society of London,1931,132(820):499-523.

[58] Miles J W.On the stability of heterogeneous shear flows[J].Journal of Fluid Mechanics,1961,10:496-508.

[59]Stacey M T,Rippeth T P,Nash J D.Turbulence and stratification in estuaries and coastal seas[M].In:Wolanski E and McLusky D(Eds.),Treatise on Estuarine and Coastal Science,2011:9-35.

猜你喜歡
差異
“再見”和bye-bye等表達的意義差異
英語世界(2023年10期)2023-11-17 09:19:16
JT/T 782的2020版與2010版的差異分析
相似與差異
音樂探索(2022年2期)2022-05-30 21:01:37
關于中西方繪畫差異及對未來發展的思考
收藏界(2019年3期)2019-10-10 03:16:40
找句子差異
DL/T 868—2014與NB/T 47014—2011主要差異比較與分析
生物為什么會有差異?
法觀念差異下的境外NGO立法效應
構式“A+NP1+NP2”與“A+NP1+(都)是+NP2”的關聯和差異
論言語行為的得體性與禮貌的差異
現代語文(2016年21期)2016-05-25 13:13:50
主站蜘蛛池模板: 国产91全国探花系列在线播放 | 亚洲日韩国产精品无码专区| 91精品人妻互换| 国产乱子伦无码精品小说 | 9966国产精品视频| 婷五月综合| 亚洲欧美在线综合一区二区三区| 国产大片黄在线观看| 国产成人精品午夜视频'| 一级香蕉人体视频| 亚洲无码高清视频在线观看| 看av免费毛片手机播放| 欧洲极品无码一区二区三区| 日本国产在线| 精品国产网站| 人人澡人人爽欧美一区| 国产内射一区亚洲| 国产一级在线观看www色 | 视频二区中文无码| 国产精品无码AⅤ在线观看播放| 免费播放毛片| 久久精品视频亚洲| 亚洲区欧美区| 精品午夜国产福利观看| 国产AV毛片| 成年人福利视频| 国产免费高清无需播放器| 国产真实乱人视频| 欧美视频在线播放观看免费福利资源| 国产香蕉在线视频| 野花国产精品入口| 国产精品播放| 99精品在线视频观看| 天天色综网| 99re在线视频观看| 日韩a在线观看免费观看| 美女视频黄频a免费高清不卡| 免费人成在线观看成人片| 国产成人精品亚洲日本对白优播| 欧洲高清无码在线| 国产成人无码播放| 亚洲视频二| 国产欧美网站| 国产99在线| 狠狠色丁香婷婷综合| 欧洲av毛片| 五月激情综合网| 国产精品19p| 国产精品部在线观看| 日韩色图在线观看| 98超碰在线观看| 国产一区二区影院| 国产91av在线| 99视频在线观看免费| av免费在线观看美女叉开腿| 国产一级视频久久| 一个色综合久久| 国产成人成人一区二区| 成人欧美日韩| 99久久精品免费看国产免费软件| 91伊人国产| 亚洲高清中文字幕| 亚洲人成网站在线播放2019| 波多野结衣一区二区三区四区视频 | 永久成人无码激情视频免费| 精品国产一区二区三区在线观看| 亚洲成人黄色网址| 日日噜噜夜夜狠狠视频| 国产亚洲高清在线精品99| 欧洲亚洲欧美国产日本高清| 欧美日韩午夜| 亚洲第一区在线| 国产精品人莉莉成在线播放| 日本三区视频| 欧美日韩在线亚洲国产人| 四虎永久在线视频| 亚洲动漫h| 国产精品第三页在线看| 在线观看国产黄色| 国产成人无码久久久久毛片| 任我操在线视频| 精品三级在线|