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

    基于Delft 3D模型的感潮河口示蹤模擬

    2020-11-05 11:16:30鄭向陽姜德娟王玉琳蔡永兵
    海洋科學(xué) 2020年10期
    關(guān)鍵詞:擴(kuò)散系數(shù)模型

    李 嘉, 鄭向陽, 張 華, 姜德娟, 王玉琳, 蔡永兵

    基于Delft 3D模型的感潮河口示蹤模擬

    李 嘉1, 2, 鄭向陽2, 3, 張 華2, 3, 姜德娟2, 3, 王玉琳1, 蔡永兵4

    (1. 揚(yáng)州大學(xué) 環(huán)境科學(xué)與工程學(xué)院, 江蘇 揚(yáng)州 225100; 2. 山東省海岸帶環(huán)境過程重點(diǎn)實(shí)驗(yàn)室, 山東 煙臺(tái) 264003; 3. 中國(guó)科學(xué)院煙臺(tái)海岸帶研究所, 山東 煙臺(tái) 264003; 4. 安徽科技學(xué)院 資源與環(huán)境學(xué)院, 安徽 鳳陽 233100)

    為研究感潮河口的水力特征及溶質(zhì)擴(kuò)散規(guī)律, 通過布放儀器監(jiān)測(cè)水文參數(shù), 分析研究區(qū)的水力特征; 采用熒光染色劑羅丹明B作為示蹤劑開展現(xiàn)場(chǎng)示蹤實(shí)驗(yàn), 研究示蹤劑的擴(kuò)散規(guī)律; 基于Delft 3D模型的水動(dòng)力模塊與示蹤模塊相耦合, 對(duì)研究區(qū)的動(dòng)力場(chǎng)和示蹤結(jié)果進(jìn)行模擬, 獲取了小清河口的水平湍流擴(kuò)散系數(shù)。結(jié)果表明: 調(diào)查期間(非汛期)小清河口的水力條件主要受潮汐控制, 潮汐對(duì)河口下游的作用更明顯。由于河流徑流量較小, 在潮汐作用下, 示蹤劑從小清河口向海傳輸?shù)乃俣容^慢, 導(dǎo)致示蹤劑在河口長(zhǎng)時(shí)間滯留。相關(guān)性分析(2)和均方根誤差(RMSE)結(jié)果表明示蹤劑遷移模擬結(jié)果的可靠性高?;诒狙芯康玫叫∏搴涌诘乃酵牧鲾U(kuò)散系數(shù)(,η)為60.12m2/s。本研究可為小清河口以及同類河口中水平湍流擴(kuò)散系數(shù)的估值提供參考, 對(duì)于評(píng)價(jià)污染物在同類河口中的傳輸行為具有指導(dǎo)意義。

    示蹤實(shí)驗(yàn); Delft 3D; 水平擴(kuò)散系數(shù); 感潮河口; 小清河

    擴(kuò)散是控制水體中污染物行為和分布的關(guān)鍵過程, 研究污染物擴(kuò)散過程是評(píng)估其在水環(huán)境中遷移的重要前提。擴(kuò)散系數(shù)與研究區(qū)的水力條件、研究區(qū)范圍和地形等因素有關(guān)。不同水域的擴(kuò)散系數(shù)通常存在較大的差異, 即使同一水域不同區(qū)段的擴(kuò)散系數(shù)也可能具有較大的變化范圍。例如, 梅江的縱向擴(kuò)散系數(shù)為18.2 m2/s[1]; 黃河不同河段的縱向擴(kuò)散系數(shù)為31~89 m2/s[2-4]; 逛蕩河不同河段的縱向擴(kuò)散系數(shù)為0.027 4~0.127 0 m2/s[5]。目前, 擴(kuò)散系數(shù)的獲取方式包括理論計(jì)算、經(jīng)驗(yàn)公式和示蹤實(shí)驗(yàn)[6]。顧莉和華祖林[6]總結(jié)發(fā)現(xiàn), 理論計(jì)算和經(jīng)驗(yàn)公式具有較大的局限性和較差的適用性; 相比之下, 通過示蹤實(shí)驗(yàn)可以較為準(zhǔn)確的獲取研究區(qū)的擴(kuò)散系數(shù), 同時(shí), 示蹤結(jié)果也可以為模型驗(yàn)證提供數(shù)據(jù)支撐[1]。

    示蹤實(shí)驗(yàn)是研究水體動(dòng)力學(xué)特征, 定量獲取示蹤劑彌散/擴(kuò)散系數(shù)的實(shí)驗(yàn)。其實(shí)驗(yàn)方法通常是在研究區(qū)上游通過特定的投放方式釋放示蹤劑, 根據(jù)現(xiàn)場(chǎng)情況在下游設(shè)置相應(yīng)的觀測(cè)站點(diǎn), 定時(shí)取樣, 分析示蹤劑濃度; 根據(jù)示蹤劑的遷移情況估算研究區(qū)的擴(kuò)散系數(shù)[5]。示蹤劑的投放方式可以分為單點(diǎn)源、多點(diǎn)源、連續(xù)釋放、瞬時(shí)釋放和有限時(shí)段釋放等[7]。示蹤劑分為反應(yīng)型示蹤劑和保守型示蹤劑, 前者在研究體系中會(huì)發(fā)生衰減; 后者在研究體系中較為穩(wěn)定。示蹤劑的選擇通常需滿足以下幾個(gè)條件: (1) 示蹤劑在研究水域中的背景濃度極低; (2) 示蹤劑在水體中的半衰期足夠長(zhǎng), 能保證在示蹤實(shí)驗(yàn)期間保持穩(wěn)定; (3)生態(tài)毒性低; (4) 檢測(cè)限低, 易于分析, 分析成本低。目前常用的幾類示蹤劑包括: 放射性物質(zhì)(如鐳、氡等同位素)、熒光染色劑(如羅丹明WT、羅丹明B、熒光素鈉和熒光黃等)、無機(jī)離子(如氯離子、溴離子和鋰離子等)。各類示蹤劑的優(yōu)缺點(diǎn)明顯, 常根據(jù)具體的實(shí)驗(yàn)條件進(jìn)行選擇。例如, 放射性物質(zhì)的環(huán)境背景值低, 但對(duì)實(shí)驗(yàn)條件要求高(實(shí)驗(yàn)人員需經(jīng)過專門培訓(xùn)), 且成本較高, 導(dǎo)致其應(yīng)用受限; 無機(jī)離子成本低、易操作, 但只適用于內(nèi)陸淡水水域, 且投加量較大; 相對(duì)而言, 熒光染色劑具有檢測(cè)限低、毒性弱、操作簡(jiǎn)單、不受研究水體限制等優(yōu)點(diǎn), 是示蹤實(shí)驗(yàn)中常用的一類示蹤劑。其中, 羅丹明類物質(zhì)(羅丹明WT、羅丹明B、磺酰羅丹明B等)是最常用的一類熒光示蹤劑。目前, 國(guó)內(nèi)外不少學(xué)者利用羅丹明類熒光染色劑在不同研究水域開展了示蹤實(shí)驗(yàn)研究[8-13]。

    Delft 3D模型是一套具有強(qiáng)大功能的水動(dòng)力-水質(zhì)模擬系統(tǒng), 主要由網(wǎng)格、水動(dòng)力、水質(zhì)、波浪、粒子示蹤、泥沙輸運(yùn)和生態(tài)7個(gè)模塊構(gòu)成。其自帶的Quickplot工具可以處理三維模擬結(jié)果, 可以將一維時(shí)間序列數(shù)據(jù)和二維空間數(shù)據(jù)進(jìn)行動(dòng)態(tài)顯示; Delft 3D模型配備網(wǎng)格編輯和生成功能, 可快速生成正交網(wǎng)格且靈活度高; 水動(dòng)力計(jì)算采用基于有限差分法中的顯、隱式交替數(shù)值積分求解, 穩(wěn)定性高, 計(jì)算速度較快; 此外, 模型的后處理能力較強(qiáng), 可以與ArcGIS和Matlab等軟件配套使用。Delft 3D模型已成為目前最先進(jìn)的水動(dòng)力數(shù)值模型之一, 被廣泛應(yīng)用于我國(guó)長(zhǎng)江口、海南島、鄱陽湖、三峽庫(kù)區(qū)和天津近海等地區(qū)[14-18]。

    研究表明, 超過90%的天然和人類活動(dòng)產(chǎn)生的物質(zhì)將隨河流輸送進(jìn)入海洋[19], 處于陸海交匯位置的河口區(qū)域作為陸源物質(zhì)向海輸送的必經(jīng)通道[20], 在陸源污染物向海輸運(yùn)過程中扮演著重要角色。由于受徑流和潮汐的雙重影響, 河口的水動(dòng)力條件特殊且復(fù)雜, 直接影響污染物在河口水體中的遷移和分布。在河口區(qū)域開展示蹤實(shí)驗(yàn), 不僅能揭示河口水動(dòng)力變化特征及其稀釋擴(kuò)散能力, 對(duì)于揭示污染物的遷移行為及其向海傳輸也具有重要意義。小清河是一條典型的陸源入海河流, 是萊州灣主要的污染來源[21-23], 其水質(zhì)好壞直接影響萊州灣區(qū)域經(jīng)濟(jì)的可持續(xù)發(fā)展和海洋資源的合理利用。因此, 本文選取小清河口作為研究區(qū), 通過釋放示蹤劑羅丹明B開展現(xiàn)場(chǎng)示蹤實(shí)驗(yàn), 采用多站位連續(xù)觀測(cè)的方式獲取示蹤劑的時(shí)空分布特征; 通過Delft 3D水動(dòng)力學(xué)模塊和粒子示蹤模塊耦合模擬示蹤劑的擴(kuò)散規(guī)律, 獲取實(shí)驗(yàn)期間小清河口的水平湍流擴(kuò)散系數(shù)。本研究有助于揭示小清河口的水動(dòng)力特征, 同時(shí)為污染物在該河口的遷移模擬提供了數(shù)據(jù)支撐。

    1 實(shí)驗(yàn)與方法

    1.1 研究區(qū)概況

    小清河下游自石村至入海口為感潮河段, 全長(zhǎng)約70 km。石村下游20 km處建有王道閘, 具有調(diào)水擋潮的功能。本研究選取王道閘至入??谧鳛檠芯繀^(qū)(圖1), 研究區(qū)全長(zhǎng)約43.6 km。小清河口河道寬度從王道閘至入??谥饾u增加, 河口呈喇叭狀向海擴(kuò)展。王道閘下游河道寬約100 m, 河口中游河道寬約150 m, 河口下游河道寬約250 m; 河道平均水深約為3 m。小清河口潮汐類型為不規(guī)則半日潮, 潮差約2 m, 受潮汐和徑流的雙重影響, 示蹤河段縱向流速季節(jié)性變化明顯, 豐水期河流徑流流速大于潮汐流速, 枯水期潮汐流速大于河流徑流流速[24]。

    圖1 研究區(qū)示意圖及觀測(cè)點(diǎn)分布

    1.2 示蹤實(shí)驗(yàn)

    2017年5月3—18日在小清河口開展現(xiàn)場(chǎng)示蹤實(shí)驗(yàn)。本文將示蹤劑泵入點(diǎn)設(shè)置在王道閘下游約1 000 m處, 5月4日上午9: 10(滿潮), 在河道中央利用蠕動(dòng)泵持續(xù)、勻速地釋放羅丹明B (RhB)溶液(濃度為20 g/L), 泵入速度約為1.6 L/min, 連續(xù)泵入總時(shí)間為5 h 10 min, RhB總投加量為10 kg。為了監(jiān)測(cè)示蹤劑的擴(kuò)散情況, 本文分別在泵入點(diǎn)下游11.5 km(觀測(cè)點(diǎn)1)、30.1 km(觀測(cè)點(diǎn)2)和43.6 km(觀測(cè)點(diǎn)3)處設(shè)置連續(xù)觀測(cè)站點(diǎn), 監(jiān)測(cè)頻率為1 h。

    實(shí)驗(yàn)開始前, 分別在3個(gè)觀測(cè)站位取背景河水, 用于配制標(biāo)準(zhǔn)曲線, 以去除河水基質(zhì)干擾。實(shí)驗(yàn)期間, 利用不銹鋼采水器在3個(gè)觀測(cè)站的河道中央定時(shí)采集表層水, 具體采樣情況如表1所示。水樣經(jīng)0.45 μm針式過濾器過濾后裝入15 mL塑料離心管, 擰緊瓶蓋, 置于加冰的保溫箱內(nèi)避光保存。樣品運(yùn)回實(shí)驗(yàn)室后, 立即用熒光光譜儀(F-7000, Hitachi, Japan)進(jìn)行分析。分析條件為: 激發(fā)波長(zhǎng)550 nm, 發(fā)射波長(zhǎng)580 nm。采用外標(biāo)法定量, 河水中RhB的檢出限為0.01 μg/L。

    為了獲取研究區(qū)的水文參數(shù), 以便為后續(xù)的模型模擬和驗(yàn)證工作提供基礎(chǔ)數(shù)據(jù), 本文在泵入點(diǎn)及3個(gè)觀測(cè)站位分別布放了相應(yīng)的觀測(cè)儀器, 各站位儀器布放情況如表1所示。

    表1 各站位儀器布放情況及采樣計(jì)劃

    注: “—”表示無采樣計(jì)劃

    1.3 示蹤模型構(gòu)建

    1.3.1 模型簡(jiǎn)介

    本文利用Delft 3D模型的FLOW模塊和PART模塊對(duì)示蹤實(shí)驗(yàn)結(jié)果進(jìn)行模擬。Delft 3D-FLOW是一個(gè)三維河流、海洋水動(dòng)力模型。由FLOW模塊計(jì)算得到的水動(dòng)力條件(流速、水位、密度、鹽度、垂直渦流和粘度等)可作為Delft3D其他模塊的輸入文件。FLOW模塊基于Boussinesq假定和靜水壓假定, 在垂向動(dòng)量方程中不考慮垂直加速度, 求解不可壓縮流體的納維-斯托克斯方程, 控制方程建立在正交曲線坐標(biāo)系()下[25]。其中, 二維水動(dòng)力數(shù)值模擬的連續(xù)性方程(1)和動(dòng)量方程(2)—(4)如下。

    垂向平均的連續(xù)性方程:

    為總水深;in和out分別為源和匯;為降水量;為蒸發(fā)量。

    垂向平均的動(dòng)量方程:

    方向:

    方向:

    分別為方向上的速度。

    ,η=at, (5)

    代表粒子的遷移時(shí)間(單位: s), 粒子被釋放時(shí),=0; 系數(shù)和(0<<1)由模型校準(zhǔn)所得。

    1.3.2 模型構(gòu)建

    1.3.2.1 構(gòu)建網(wǎng)格和地形數(shù)據(jù)

    在RGFGRID模塊中選擇笛卡爾坐標(biāo)系, 導(dǎo)入邊界文件, 分別在與邊界平行和垂直的方向手動(dòng)添加樣條曲線(spline), 生成粗略的網(wǎng)格; 然后在每個(gè)網(wǎng)格內(nèi)生成等數(shù)量的加密網(wǎng)格(圖2), 檢驗(yàn)網(wǎng)格正交性后, 保存為網(wǎng)格文件(.grd)。研究區(qū)的總網(wǎng)格數(shù)為12×1 472。小清河口的水深散點(diǎn)數(shù)據(jù)來自兩部分: 王道閘至羊口港段的水深為實(shí)測(cè)值, 羊口港以下河段的水深值取自海圖。利用Arcgis對(duì)散點(diǎn)水深數(shù)據(jù)進(jìn)行數(shù)字化處理, 導(dǎo)出坐標(biāo)值與深度值, 保存成水深文件(.xyz文件); 使用三角插值選項(xiàng)進(jìn)行插值, 然后進(jìn)行內(nèi)部擴(kuò)散和平滑處理, 將插值處理后的文件保存成地形文件(.dep)。結(jié)果如圖2所示。

    圖2 小清河口水深地形圖

    1.3.2.2 構(gòu)建水動(dòng)力模型

    將生成的邊界文件、網(wǎng)格文件和地形文件導(dǎo)入FLOW模塊, 設(shè)置模擬時(shí)間段為2017年3月20日— 2017年5月19日, 參照時(shí)間為2017年3月1日, 時(shí)間步長(zhǎng)設(shè)置為30 min; 其他模型參數(shù)均采用默認(rèn)值。由于缺乏模擬區(qū)上邊界的實(shí)測(cè)流量數(shù)據(jù), 本文基于泵入點(diǎn)的實(shí)測(cè)流速(5月3日至5月17日)與河道橫截面積估算河流的徑流量; 其他時(shí)段的河流徑流量取平均值8.3 m3/s。模型的開邊界位于觀測(cè)站3, 使用該站實(shí)測(cè)水深數(shù)據(jù)作為開邊界條件。

    1.3.2.3 構(gòu)建粒子示蹤模型

    將FLOW模塊生成的水動(dòng)力文件導(dǎo)入PART模塊。將模擬目標(biāo)設(shè)置為示蹤劑(Tracer), 設(shè)置粒子數(shù)量為10 000。釋放方式為連續(xù)釋放, 同時(shí)設(shè)置釋放點(diǎn)位置、釋放時(shí)間(5月4日09: 10: 00—14: 10: 00)與示蹤劑濃度(20 g/L); 調(diào)整擴(kuò)散系數(shù)和示蹤劑衰減常數(shù); 設(shè)置觀測(cè)站位; 最后設(shè)置輸出文件的時(shí)間范圍和步長(zhǎng)。

    1.4 模型驗(yàn)證

    利用示蹤實(shí)驗(yàn)期間在各觀測(cè)站位獲取的水文數(shù)據(jù)進(jìn)行流速、流向和水深的驗(yàn)證, 利用示蹤劑實(shí)測(cè)濃度驗(yàn)證示蹤模擬結(jié)果。分別采用相關(guān)性系數(shù)(R)和均方根誤差(RMSE)評(píng)估模型模擬結(jié)果。

    其中,sim為示蹤劑的模擬濃度值;obs為示蹤劑的實(shí)測(cè)濃度值;為數(shù)據(jù)量。

    2 結(jié)果與討論

    2.1 基于實(shí)測(cè)水文參數(shù)的小清河口水動(dòng)力特征分析

    河口水體受河水徑流和潮汐的雙重影響, 其水深和電導(dǎo)率呈現(xiàn)周期性變化。因此, 電導(dǎo)率和水深是反映河口水體運(yùn)動(dòng)的重要水文參數(shù)。由于海水的電導(dǎo)率遠(yuǎn)大于河水, 海水上溯可引起河口水體電導(dǎo)率增大, 同時(shí)伴隨水深升高。如圖3所示, 泵入點(diǎn)處的電導(dǎo)率總體呈現(xiàn)上升的趨勢(shì), 最大值為16 mS/cm; 且此處的電導(dǎo)率與水深變化無關(guān), 水深降低時(shí)電導(dǎo)率緩慢增大, 水深增加時(shí)電導(dǎo)率急劇增加。這表明潮汐不是控制此處水體電導(dǎo)率的單一因素。RhB分子中含有氯離子, 當(dāng)其溶于水后, 會(huì)電離出氯離子, 引起周圍水體電導(dǎo)率增加。在RhB釋放期間(09: 10—14: 10), 盡管時(shí)處落潮, 河水電導(dǎo)率仍隨RhB的泵入緩慢增加; 釋放結(jié)束后, 漲潮時(shí)(16: 10—19: 10), 高電導(dǎo)率的RhB水團(tuán)被推回泵入點(diǎn), 加之海水上溯, 導(dǎo)致此時(shí)水體電導(dǎo)率急劇增加。此外, 同樣是漲潮, 當(dāng)潮能較大時(shí)(03: 10—07: 10), 更多的RhB水團(tuán)被推回泵入點(diǎn)以及更多的海水上溯, 導(dǎo)致此時(shí)電導(dǎo)率增加的更明顯。這表明非汛期時(shí)潮汐對(duì)小清河口上游的水動(dòng)力也有較強(qiáng)的影響。

    圖3 不同站位的水文參數(shù)

    由于觀測(cè)站1沒有布放水位計(jì), 無法獲取水深實(shí)測(cè)值, 因此本文使用流速方向指示潮汐變化規(guī)律。沿河方向的流速有正負(fù)之分, 其中正值代表流向向海, 即落潮; 負(fù)值表示流向向岸, 即漲潮。當(dāng)流速值由正值逐漸變?yōu)?時(shí), 水深值最小; 當(dāng)流速值由負(fù)值逐漸變?yōu)?時(shí), 水深最大。如圖3所示, 該處的電導(dǎo)率與流速方向之間的關(guān)系并不固定。5月6日至5月10日, 落潮時(shí)(流速為正), 電導(dǎo)率增加; 漲潮時(shí)(流速為負(fù)), 電導(dǎo)率降低。表明此時(shí)該處的電導(dǎo)率主要受上游高電導(dǎo)率的河水(RhB電離)控制。5月10日后, RhB水團(tuán)遷移至觀測(cè)站1的下游, 漲潮時(shí)(流速為負(fù))電導(dǎo)率增加。如圖3所示, 觀測(cè)站2和3的電導(dǎo)率與其水深變化規(guī)律一致, 受潮能變化的影響兩個(gè)站位電導(dǎo)率的峰值呈現(xiàn)高低起伏的變化趨勢(shì), 表明這兩個(gè)站位的電導(dǎo)率完全由潮汐所控制。因?yàn)檫@兩個(gè)站位距離泵入點(diǎn)足夠遠(yuǎn)(>30 km), RhB水團(tuán)在遷移過程中不斷被稀釋, 當(dāng)其到達(dá)觀測(cè)站2和3時(shí), 并沒有對(duì)此處水體的電導(dǎo)率構(gòu)成影響。此外, 觀測(cè)站3處水體的電導(dǎo)率明顯高于觀測(cè)站2, 表明越接近河口下游, 潮汐的作用越明顯。結(jié)合4個(gè)站位水文參數(shù)的變化規(guī)律可發(fā)現(xiàn), 非汛期時(shí)小清河口的水動(dòng)力主要受潮汐控制, 尤其是河口下游區(qū)域。

    2.2 示蹤劑在小清河口的擴(kuò)散規(guī)律分析

    觀測(cè)站1相距泵入點(diǎn)11.5 km, 示蹤劑的質(zhì)量濃度峰值出現(xiàn)在5月6日至5月7日, 隨后RhB的濃度維持在一定水平波動(dòng), 5月12日RhB的濃度明顯降低(圖4)。觀測(cè)站2距離泵入點(diǎn)約30 km, RhB到達(dá)此處的時(shí)間大約為5月10日; 5月15日觀測(cè)站2的RhB濃度達(dá)到最大值; 5月18日, 即實(shí)驗(yàn)結(jié)束時(shí)仍能在觀測(cè)站2檢測(cè)到RhB(圖4)。因此實(shí)驗(yàn)期間, 示蹤劑并沒有到達(dá)觀測(cè)站3??傮w上, 示蹤劑在小清河口向海傳輸?shù)乃俣容^慢, 釋放14 d后, RhB仍在距釋放點(diǎn)約30 km的區(qū)域回蕩。這主要是因?yàn)槭聚檶?shí)驗(yàn)期間小清河正值枯水期, 河流徑流量較小, 在潮汐作用下, 示蹤劑向海擴(kuò)散受阻。由此可推斷出, 非汛期時(shí), 來自河流輸入的污染物容易在小清河口滯留, 對(duì)河口生態(tài)系統(tǒng)構(gòu)成巨大威脅。

    圖4 示蹤劑的時(shí)空分布規(guī)律

    2.3 小清河口水動(dòng)力模擬結(jié)果分析

    小清河口水動(dòng)力模擬時(shí)段為2017年3月20日0點(diǎn)至2017年5月19日0點(diǎn)。為了驗(yàn)證小清河口二維水動(dòng)力模型的準(zhǔn)確性, 本文分別對(duì)各站位的水深和流速的實(shí)測(cè)值與模擬值進(jìn)行了相關(guān)性分析。

    泵入點(diǎn)的水深實(shí)測(cè)值與模擬值如圖5所示, 二者的變化趨勢(shì)一致(<0.01), 相關(guān)系數(shù)(2)為0.878, 且水深實(shí)測(cè)值與模擬值在數(shù)量級(jí)上也較為接近。5月11日, 大潮引起水深明顯增加, 二維動(dòng)力學(xué)模型很好地模擬出這一過程。由于觀測(cè)站1缺少水深實(shí)測(cè)值, 故無法進(jìn)行對(duì)比。如圖5所示, 觀測(cè)站2的水深實(shí)測(cè)值與模擬值也呈現(xiàn)高度一致性(2=0.955,<0.01); 對(duì)于5月4日至5月7日觀測(cè)站2水深的不規(guī)則波動(dòng), 該模型也能很好的模擬。D站位的水深觀測(cè)值作為模型的開邊界條件, 因此模擬值與實(shí)測(cè)值也基本重合。

    圖5 各站位水深與流速的實(shí)測(cè)值和模擬值

    如圖5所示, 泵入點(diǎn)處的流速在東-西方向上以正值為主, 即為向海方向(東方向), 表明此處流向主要受河流徑流控制。泵入點(diǎn)處東-西向的流速實(shí)測(cè)值與模擬值呈高度一致性(2=0.956,<0.01)。觀測(cè)站1處的流速實(shí)測(cè)值與模擬值如圖5所示, 二維水動(dòng)力模型可以較好地模擬該站東-西方向的流速(2=0.668,<0.01)。觀測(cè)站2處的海流計(jì)發(fā)生故障, 僅記錄了5月4日至10日之間的流速數(shù)據(jù)(圖5)。結(jié)果表明, 觀測(cè)站2處的東-西向流速模擬值與實(shí)測(cè)值擬合度較高(2=0.667,<0.01)。

    總體上, 水深和流速的模擬結(jié)果與實(shí)測(cè)值擬合度較好。表明本文基于Delft 3D-FLOW模塊構(gòu)建的二維動(dòng)力學(xué)模型能夠提供較準(zhǔn)確的動(dòng)力場(chǎng), 該動(dòng)力學(xué)模型可以與獨(dú)立的示蹤模塊或水質(zhì)模塊進(jìn)行耦合, 進(jìn)而實(shí)現(xiàn)小清河口目標(biāo)物質(zhì)的遷移和歸趨模擬。

    2.4 小清河口示蹤模擬分析

    由于采樣期間, 示蹤劑未遷移至觀測(cè)站3。所以本文僅選取觀測(cè)站1和觀測(cè)站2的示蹤劑實(shí)測(cè)值來驗(yàn)證示蹤模型的模擬結(jié)果。兩個(gè)觀測(cè)站位的示蹤劑實(shí)測(cè)值與模擬值如圖6所示。觀測(cè)站1的示蹤劑模擬值與實(shí)測(cè)值較為接近, 優(yōu)化后的模型可以很好地再現(xiàn)5月6日至5月8日的濃度峰; 5月8日至5月15日的低濃度峰, 該模型也能較好地?cái)M合。相關(guān)性分析結(jié)果表明, 觀測(cè)站1的示蹤劑模擬值與實(shí)測(cè)值之間呈顯著正相關(guān)(<0.01), 相關(guān)系數(shù)(2)為0.717(圖6), 均方根誤差(RMSE)為0.176。如圖6所示, 優(yōu)化后的模型也能很好地?cái)M合觀測(cè)站2的示蹤劑濃度變化規(guī)律。該示蹤模型可以很好地再現(xiàn)示蹤劑峰值出現(xiàn)的時(shí)間, 模擬值在量級(jí)上與實(shí)測(cè)值也較為接近。相關(guān)性分析結(jié)果表明, 觀測(cè)站2的示蹤劑模擬值與實(shí)測(cè)值之間呈顯著正相關(guān)(<0.01), 二者之間的相關(guān)系數(shù)2=0.707(圖6), 均方根誤差(RMSE)為0.036。

    圖6 各站位示蹤劑的實(shí)測(cè)值和模擬值及相關(guān)性分析

    經(jīng)過模型參數(shù)調(diào)試, 小清河口示蹤劑的水平擴(kuò)散系數(shù)被確定為,η=0.12m2/s。在Delft 3D-PART模塊中, 認(rèn)為水平擴(kuò)散是由湍流引起的, 根據(jù)湍流理論, 擴(kuò)散系數(shù)隨時(shí)間增大[26]。盡管如此, 由于本文得出的水平湍流擴(kuò)散系數(shù)的冪指數(shù)僅為0.12, 隨著時(shí)間延長(zhǎng), 水平擴(kuò)散系數(shù)的增速極慢。在本研究中, 從釋放示蹤劑開始(5月4日09: 10)至實(shí)驗(yàn)結(jié)束(5月18日10: 00), 實(shí)驗(yàn)時(shí)長(zhǎng)約1 212 600 s。由本文得出的水平湍流擴(kuò)散參數(shù)的計(jì)算公式(,η= 60.12)可求得此時(shí)的水平擴(kuò)散系數(shù)約為32.2 m2/s。李林娟和童朝鋒[16]在基于Delft 3D模型模擬長(zhǎng)江口鹽度擴(kuò)散規(guī)律的研究中(模擬周期31 d), 將水平湍流擴(kuò)散系數(shù)的值設(shè)置為50~150 m2/s。羅家海[27]在計(jì)算深圳灣水平擴(kuò)散系數(shù)的案例中指出, 當(dāng)把潮流看作湍流時(shí), 水平擴(kuò)散系數(shù)的量級(jí)為102~103。此外, 匡國(guó)瑞等[28]也指出, 渤海的水平湍流擴(kuò)散系數(shù)在102m2/s范圍內(nèi)。通過對(duì)比可以發(fā)現(xiàn), 本研究所獲取的水平湍流擴(kuò)散系數(shù)的值較為合理。本文的結(jié)果可以為小清河口以及同類河口中水平湍流擴(kuò)散系數(shù)的估值提供數(shù)據(jù)參考, 同時(shí)對(duì)于評(píng)價(jià)污染物在同類河口中的傳輸行為具有指導(dǎo)意義。

    3 結(jié)論

    RhB被釋放后導(dǎo)致周圍水體電導(dǎo)率增加, 因此會(huì)干擾感潮河口部分區(qū)域水體的電導(dǎo)率變化規(guī)律。非汛期時(shí), 潮汐是控制小清河口水力特征的關(guān)鍵因素; 且越接近河口下游, 潮汐的作用越明顯, 表現(xiàn)為河口下游東-西向流速的絕對(duì)值較為接近。河口上游水體的流向以向海方向?yàn)橹? 表明此處河流徑流的作用強(qiáng)于潮汐。由于非汛期河流徑流較弱, 示蹤劑在小清河口向海傳輸?shù)乃俣容^慢。在潮汐的作用下, 釋放14 d后, RhB仍在距離釋放點(diǎn)下游30 km的區(qū)域回蕩。模型驗(yàn)證結(jié)果表明, 本文基于Delft 3D-FLOW模塊構(gòu)建的二維動(dòng)力學(xué)模型能夠準(zhǔn)確地模擬小清河口的動(dòng)力場(chǎng), 為后續(xù)的示蹤模擬提供了可靠的動(dòng)力學(xué)文件。通過Delft 3D-FLOW模塊和Delft 3D-PART模塊耦合, 本文成功模擬了不同觀測(cè)站位示蹤劑濃度的變化規(guī)律。示蹤劑模擬值與實(shí)測(cè)值呈顯著正相關(guān)(<0.01), 相關(guān)系數(shù)大于0.7。經(jīng)過模型參數(shù)調(diào)整, 小清河口水平湍流擴(kuò)散系數(shù)確定為,η= 60.12m2/s。通過與其他文獻(xiàn)報(bào)道值對(duì)比可得出, 在本次示蹤實(shí)驗(yàn)期間, 基于Delft 3D模型所獲取的水平湍流擴(kuò)散系數(shù)是合理的。

    [1] 吳群河. 梅江染料示蹤實(shí)驗(yàn)[J]. 中山大學(xué)學(xué)報(bào), 1990, 9(1): 25-29. Wu Qunhe. A tracing experiment in Meijiang River[J]. Journal of Sun Yatsen University, 1990, 9(1): 25-29.

    [2] 呂巍, 張建軍, 閆莉, 等. 黃河潼關(guān)河段污染物示蹤擴(kuò)散研究[J]. 人民黃河, 2016, 38(9): 59-63. Lv Wei, Zhang Jianjun, Yan Li, et al. Research of the pollutants tracer experiment in the Tongguan reach of the Yellow River[J]. Yellow River, 2016, 38(9): 59-63.

    [3] 肖翔群, 連煜, 胡國(guó)華, 等. 黃河孟津段擴(kuò)散系數(shù)的確定[J]. 水資源保護(hù), 1996, 12(1): 38-40. Xiao Xiangqun, Lian Yu, Hu Guohua, et al. Determination of diffusion coefficient in Mengjin reach of the Yellow River[J]. Water Resources Protection, 1996, 12(1): 38-40.

    [4] 邵慶山, 張紅軍. 多沙大河的擴(kuò)散示蹤研究[J]. 環(huán)境科學(xué)進(jìn)展, 1994, 2(3): 65-70. Shao Qingshan, Zhang Hongjun. Study on diffusion tracer of sandy river[J]. Progress in Environmental Science, 1994, 2(3): 65-70.

    [5] 于靖, 張華. 城市小型河流水動(dòng)力彌散和潛流交換過程[J]. 水科學(xué)進(jìn)展, 2015, 26(5): 714-721. Yu Jing, Zhang Hua. Hydrodynamic dispersion and hyporheic exchange in a small urban stream[J]. Advances in Water Science, 2015, 26(5): 714-721.

    [6] 顧莉, 華組林. 天然河流縱向離散系數(shù)確定方法的研究進(jìn)展[J]. 水利水電科技進(jìn)展, 2007, 27(2): 85-89. Gu Li, Hua Zulin. Advances in determination of longitudinal dispersion coefficient of natural rivers[J]. Advances in Science and Technology of Water Resources, 2007, 27(2): 85-89.

    [7] 周世良. 水?dāng)U散示蹤實(shí)驗(yàn)在環(huán)評(píng)中的應(yīng)用[J]. 海峽科學(xué), 2009(6): 3-4, 8. Zhou Shiliang. Application of water diffusion tracer test in environmental impact assessment[J]. Straits Science, 2009(6): 3-4, 8.

    [8] Hiatt M, Passalacqua P. Hydrological connectivity in river deltas: The first-order importance of channel-is-land exchange[J]. Water Resources Research, 2015, 51(4): 2264-2282.

    [9] Williams C F, Nelson S D. Comparison of Rhodamine-WT and bromide as a tracer for elucidating internal wetland flow dynamics[J]. Ecological Engineering, 2011, 37(10): 1492-1498.

    [10] Wolkersdorfer C. Tracer test in a settling pond: The passive mine water treatment plant of the 1 B mine pool, Nova Scotia, Canada[J]. Mine Water and the Environment, 2011, 30(2): 105-112.

    [11] 李培泉. 應(yīng)用羅丹明-B作示蹤劑研究海水的稀釋擴(kuò)散規(guī)律[J]. 海洋科學(xué), 1981, 5(2): 38.Li Peiquan. Study on the dilution and diffusion of seawater using Rhodamine B[J]. Marine Sciences, 1981, 5(2): 38.

    [12] 梁健. 污染擴(kuò)散的示蹤劑羅丹明-B的試驗(yàn)[J]. 重慶環(huán)境科學(xué), 1990, 12(6): 41-44. Liang Jian. Test of contamination diffusion using tracer rhodamine B[J]. Chongqing Environmental Science, 1990, 12(6): 41-44.

    [13] 聶艷華, 段文剛, 樹錦. 示蹤法定量分析水流連通問題[J]. 長(zhǎng)江科學(xué)院院報(bào), 2013, 30(2): 16-19. Nie Yanhua, Duan Wengang, Shu Jin. Quantitative analysis of flow connectivity by tracing method[J]. Journal of Yangtze River Scientific Research Institute, 2013, 30(2): 16-19.

    [14] 范翻平. 基于Delft3D模型的鄱陽湖水動(dòng)力模擬研究[D].南昌: 江西師范大學(xué), 2010. Fan Fanping. The study of hydrodynamic simulation of Poyang Lake based on Delft 3D model[D]. Nanchang: Jiangxi Normal University, 2010.

    [15] 龔文平, 李昌宇, 林國(guó)堯, 等. Delft 3D在離岸人工島建設(shè)中的應(yīng)用: 以海南島萬寧日月灣人工島為例[J]. 海洋工程, 2012, 30(3): 35-44. Gong Wenping, Li Changyu, Lin Guoyao, et al. Application of Delft 3D model for plan design of an artificial island: A case study for the artificial island construction in the Riyue Bay, Wanning City, Hainan Island[J]. The Ocean Engineering, 2012, 30(3): 35-44.

    [16] 李林娟, 童朝鋒. 基于Delft3D-Flow模型的長(zhǎng)江口鹽度擴(kuò)散規(guī)律模擬[J]. 人民長(zhǎng)江, 2016, 47(23): 107-111. Li Linjuan, Tong Chaofeng. Simulation of salinity diffusion of Yangtze River estuary based on Delft3D-flow model[J]. Yangtze River, 2016, 47(23): 107-111.

    [17] 陸仁強(qiáng), 何璐珂. 基于 Delft3D模型的近海水環(huán)境質(zhì)量數(shù)值模擬研究[J]. 海洋環(huán)境科學(xué), 2012, 31(6): 877-880. Lu Renqiang, He Luke. Study on numerical simulation of environmental quality in coastal water based on Delft 3D model[J]. Marine Environmental Science, 2012, 31(6): 877-880.

    [18] 黃慶超, 石巍方, 劉廣龍, 等. 基于Delft3D的三峽水庫(kù)不同工況下香溪河水動(dòng)力水質(zhì)模擬[J]. 水資源與工程學(xué)報(bào), 2017, 28(2): 33-39. Huang Qingchao, Shi Weifang, Liu Guanglong, et al. Modeling the hydrodynamics and water quality of Xiangxi River under different working conditions of Three Gorges Reservoir based on Delft 3D[J]. Journal of Water Resources & Water Engineering, 2017, 28(2): 33-39.

    [19] Abdel-moati A R. Behaviour and fluxes of copper and lead in the Nile River estuary[J]. Estuarine Coastal and Shelf Science, 1990, 30(2): 153-165.

    [20] Simpson J H, Vennell R, Souza A J. The salt fluxes in a tidally-energetic estuary[J]. Estuarine Coastal and Shelf Science, 2001, 52(1): 131-142.

    [21] 李嘉, 張瑞杰, 王潤(rùn)梅, 等. 小清河流域抗生素污染分布特征與生態(tài)風(fēng)險(xiǎn)評(píng)估[J]. 農(nóng)業(yè)環(huán)境科學(xué)學(xué)報(bào), 2016, 35(7): 1384-1391. Li Jia, Zhang Ruijie, Wang Runmei, et al. Distribution characteristics and ecological risk assessment of antibiotic pollution in Xiaoqing River watershed[J]. Journal of Agro-Environment Science, 2016, 35(7): 1384-1391.

    [22] Shen J Y, Luo X X, Zheng H, et al. Pollution and ecological risk characteristics of heavy metals in surface sediments in Xiaoqing River Estuary and adjacent sea areas[J]. Environmental Chemistry, 2017, 36(7): 1516-1524.

    [23] Shi Y L, Vestergren R, Xu L, et al. Characterizing direct emissions of perfluoroalkyl substances from ongoing fluoropolymer production sources: A spatial trend study of Xiaoqing River, China[J]. Environmental Pollution, 2015, 206: 104-112.

    [24] Zou T, Zhang H, Meng Q J, et al. Seasonal hydrodynamics and salt exchange of a shallow estuary in northern China[J]. Journal of Coastal Research, 2016, 74: 95-103.

    [25] Deltares. Delft3D-FLOW user manual version 3.14[R]. Netherlands: Delft, 2009.

    [26] Deltares. Delft3D-PART user manual version 2.13[R]. Netherlands: Delft, 2009.

    [27] 羅家海. 河口及海灣地區(qū)水平擴(kuò)散系數(shù)確定研究[J]. 海洋環(huán)境科學(xué), 1989, 8(3): 36-40. Luo Jiahai. Study on determination of horizontal diffusion coefficient in estuary and bay[J]. Marine Environmental Science, 1989, 8(3): 36-40.

    [28] 匡國(guó)瑞, 俞光耀, 張淮, 等. 現(xiàn)場(chǎng)海域水平、鉛直擴(kuò)散系數(shù)的推算[J]. 海洋學(xué)報(bào), 1994, 16(4): 23-34.Kuang Guorui, Yu Guangyao, Zhang Huai, et al. Estimation of horizontal and vertical diffusion coefficient in Bohai sea[J]. Acta Oceanologica Sinica, 1994, 16(4): 23-34.

    Simulation of tracer experiment in the tidal estuary based on the Delft 3D model

    LI Jia1, 2, ZHENG Xiang-yang2, 3, ZHANG Hua2, 3, JIANG De-juan2, 3, WANG Yu-lin1, CAI Yong-bing4

    (1. School of Environmental Science and Engineering, Yangzhou University, Yangzhou 225100, China; 2. Shandong Key Laboratory of Coastal Environmental Processes, Yantai 264003, China; 3. Yantai Institute of Coastal Zone Research, Chinese Academy of Sciences, Yantai 264003, China; 4. College of Resource and Environment, Anhui Science and Technology University, Anhui 233100, China)

    To reveal the hydraulic characteristics and the diffusion process of solutes from the tidal estuary, an experiment was carried out at the Xiaoqing River estuary using a field tracker using Rhodamine B and hydrological survey. A coupled model based on Delft3D-FLOW and Delft3D-PART was used to simulate the flow field of the estuary and the transport process of the tracer experiment. The horizontal turbulent diffusion coefficient (,η) was estimated based on the Delft-3D simulation. According to the results, it can be concluded that the hydraulic conditions of the Xiaoqing River estuary, especially the lower ranges of the tidal estuary, during the experimental period were determined by the tides. The velocity of longitudinal transport of the tracer in the Xiaoqing River estuary was slow due to the strong flow of the tides and the weak current of the river. Likewise, the tracer’s residence time would be long in the research area. The Delft3D-FLOW module provided a reliable hydrodynamic file for the simulation area. Correlation analysis (2) and root mean square error (RMSE) showed an excellent fit of the model. The horizontal turbulent diffusion coefficient (,η) was 60.12m2/s. This research could provide data and guidelines for subsequent water quality simulation work in the Xiaoqing River estuary.

    racer experiment; Delft 3D; horizontal turbulent diffusion coefficient; tidal estuary; the Xiaoqing River

    Nov. 12, 2019

    X143

    A

    1000-3096(2020)10-0023-10

    10.11759/hykx20191112002

    2019-11-12;

    2019-12-26

    國(guó)家自然科學(xué)基金項(xiàng)目(41671473); 山東省海岸帶環(huán)境過程重點(diǎn)實(shí)驗(yàn)室(中國(guó)科學(xué)院煙臺(tái)海岸帶研究所)開放基金(2019SDHADKFJJ12)

    [National Nature Science Foundation of China, No. 41671473; Shandong Key Laboratory of Coastal Environmental Processes, YICCAS, No. 2019SDHADKFJJ12]

    李嘉(1990-), 男, 山東煙臺(tái)人, 講師, 博士, 主要從事有機(jī)污染物遷移轉(zhuǎn)化模擬方面研究, E-mail: lijia3611@yzu.edu.cn; 張華,通信作者, 研究員, 主要從事水環(huán)境模擬研究, E-mail: hzhang@yic.ac.cn

    (本文編輯: 劉珊珊)

    猜你喜歡
    擴(kuò)散系數(shù)模型
    一半模型
    重要模型『一線三等角』
    重尾非線性自回歸模型自加權(quán)M-估計(jì)的漸近分布
    一類具有變擴(kuò)散系數(shù)的非局部反應(yīng)-擴(kuò)散方程解的爆破分析
    3D打印中的模型分割與打包
    基于Sauer-Freise 方法的Co- Mn 體系fcc 相互擴(kuò)散系數(shù)的研究
    上海金屬(2015年5期)2015-11-29 01:13:59
    FCC Ni-Cu 及Ni-Mn 合金互擴(kuò)散系數(shù)測(cè)定
    上海金屬(2015年6期)2015-11-29 01:09:09
    FLUKA幾何模型到CAD幾何模型轉(zhuǎn)換方法初步研究
    非時(shí)齊擴(kuò)散模型中擴(kuò)散系數(shù)的局部估計(jì)
    Ni-Te 系統(tǒng)的擴(kuò)散激活能和擴(kuò)散系數(shù)研究
    上海金屬(2013年4期)2013-12-20 07:57:07
    精品电影一区二区在线| 黄色成人免费大全| 操出白浆在线播放| 国产一区二区三区在线臀色熟女| 国产精品久久久人人做人人爽| 午夜福利,免费看| 国产av一区二区精品久久| 免费在线观看视频国产中文字幕亚洲| www国产在线视频色| 高清毛片免费观看视频网站| 欧美在线一区亚洲| 欧美日韩乱码在线| 老熟妇仑乱视频hdxx| 真人做人爱边吃奶动态| 免费搜索国产男女视频| 岛国视频午夜一区免费看| 乱人伦中国视频| 欧美日韩瑟瑟在线播放| 中文字幕色久视频| 午夜免费激情av| 90打野战视频偷拍视频| 午夜福利18| 侵犯人妻中文字幕一二三四区| 精品第一国产精品| 国产区一区二久久| 99国产精品一区二区蜜桃av| 国产精品日韩av在线免费观看 | 每晚都被弄得嗷嗷叫到高潮| 九色国产91popny在线| 欧美成人一区二区免费高清观看 | av在线播放免费不卡| 国产成人系列免费观看| 波多野结衣av一区二区av| 18禁观看日本| 777久久人妻少妇嫩草av网站| 久久这里只有精品19| 国产精品av久久久久免费| 欧美日本亚洲视频在线播放| 91字幕亚洲| 又黄又粗又硬又大视频| 久久青草综合色| 国产精品电影一区二区三区| 国产成人免费无遮挡视频| 一级毛片女人18水好多| 日韩大尺度精品在线看网址 | 天天躁狠狠躁夜夜躁狠狠躁| 国产不卡一卡二| 亚洲国产看品久久| avwww免费| 国产主播在线观看一区二区| 精品一品国产午夜福利视频| 精品国内亚洲2022精品成人| 嫁个100分男人电影在线观看| 亚洲精品久久成人aⅴ小说| 精品不卡国产一区二区三区| 日本在线视频免费播放| 欧美黄色淫秽网站| 母亲3免费完整高清在线观看| 日日干狠狠操夜夜爽| 母亲3免费完整高清在线观看| 麻豆av在线久日| www.www免费av| 一本久久中文字幕| 桃红色精品国产亚洲av| 亚洲色图 男人天堂 中文字幕| 免费看十八禁软件| 男女午夜视频在线观看| 国产成人精品久久二区二区免费| 老司机午夜福利在线观看视频| 露出奶头的视频| 两个人免费观看高清视频| 亚洲 国产 在线| 99久久久亚洲精品蜜臀av| 日本欧美视频一区| 久久精品国产综合久久久| av天堂久久9| 淫妇啪啪啪对白视频| 欧美中文综合在线视频| 人人妻人人澡人人看| 久久九九热精品免费| 很黄的视频免费| 午夜久久久在线观看| 久久人妻av系列| 啦啦啦观看免费观看视频高清 | 脱女人内裤的视频| 首页视频小说图片口味搜索| 久9热在线精品视频| av福利片在线| 欧美乱妇无乱码| 满18在线观看网站| 色综合站精品国产| 亚洲熟妇熟女久久| 在线天堂中文资源库| 一二三四社区在线视频社区8| 免费少妇av软件| 成人特级黄色片久久久久久久| 99国产精品一区二区三区| 韩国av一区二区三区四区| 欧美国产精品va在线观看不卡| 真人做人爱边吃奶动态| svipshipincom国产片| 极品人妻少妇av视频| 男女下面插进去视频免费观看| 少妇裸体淫交视频免费看高清 | 国产高清有码在线观看视频 | 国产精品久久视频播放| 日韩欧美免费精品| 女生性感内裤真人,穿戴方法视频| 男人舔女人下体高潮全视频| 亚洲国产精品999在线| 午夜福利18| 身体一侧抽搐| 国产精品亚洲一级av第二区| 两性午夜刺激爽爽歪歪视频在线观看 | 久久精品影院6| 日韩欧美国产一区二区入口| 此物有八面人人有两片| 久久久久久免费高清国产稀缺| 女同久久另类99精品国产91| 久久久久久久久中文| 美女 人体艺术 gogo| 麻豆av在线久日| 国产亚洲av嫩草精品影院| 国产亚洲欧美精品永久| 久久国产精品人妻蜜桃| 成人欧美大片| 九色亚洲精品在线播放| 麻豆久久精品国产亚洲av| 最好的美女福利视频网| 午夜免费观看网址| 国产亚洲精品久久久久5区| 一个人观看的视频www高清免费观看 | 一级作爱视频免费观看| 国产成人精品久久二区二区91| 亚洲无线在线观看| 欧美激情 高清一区二区三区| 一区二区日韩欧美中文字幕| 丁香欧美五月| 亚洲自拍偷在线| 国产xxxxx性猛交| 男女午夜视频在线观看| 亚洲伊人色综图| 国产欧美日韩一区二区三区在线| 一区二区三区激情视频| 亚洲精品中文字幕在线视频| 亚洲中文字幕日韩| 少妇 在线观看| 亚洲欧美精品综合一区二区三区| 超碰成人久久| 亚洲一区二区三区不卡视频| 热re99久久国产66热| 中文字幕av电影在线播放| 日韩大尺度精品在线看网址 | 国产免费av片在线观看野外av| av在线播放免费不卡| 国产不卡一卡二| 亚洲国产看品久久| 国产又色又爽无遮挡免费看| 啦啦啦免费观看视频1| 欧美色视频一区免费| 国产成人精品无人区| 免费高清在线观看日韩| av有码第一页| 亚洲av日韩精品久久久久久密| av有码第一页| 极品教师在线免费播放| 欧美日韩乱码在线| 成人免费观看视频高清| 国内毛片毛片毛片毛片毛片| 国产成人精品在线电影| 一个人免费在线观看的高清视频| 久久香蕉国产精品| 亚洲av电影在线进入| 搡老妇女老女人老熟妇| 久久人妻熟女aⅴ| 国产av精品麻豆| 午夜亚洲福利在线播放| 精品国产乱码久久久久久男人| 国产精品爽爽va在线观看网站 | 国产精品亚洲一级av第二区| 亚洲欧美激情在线| www.自偷自拍.com| av天堂在线播放| 亚洲国产精品999在线| 啦啦啦免费观看视频1| 一区福利在线观看| 欧美黑人欧美精品刺激| 老汉色av国产亚洲站长工具| 19禁男女啪啪无遮挡网站| 一二三四社区在线视频社区8| 久久久久久久午夜电影| 亚洲无线在线观看| av在线天堂中文字幕| 精品免费久久久久久久清纯| 亚洲av熟女| 亚洲五月婷婷丁香| 一区二区三区精品91| 欧美乱妇无乱码| 非洲黑人性xxxx精品又粗又长| 精品无人区乱码1区二区| 色综合欧美亚洲国产小说| 国产成人啪精品午夜网站| 最新美女视频免费是黄的| 久久青草综合色| 9191精品国产免费久久| 国产免费男女视频| 亚洲午夜理论影院| 一本综合久久免费| 亚洲午夜精品一区,二区,三区| 亚洲欧美日韩另类电影网站| 亚洲三区欧美一区| 日韩免费av在线播放| 午夜精品国产一区二区电影| 嫩草影视91久久| 精品国产一区二区三区四区第35| 国产av一区二区精品久久| 亚洲久久久国产精品| 亚洲av日韩精品久久久久久密| 18禁美女被吸乳视频| 男人操女人黄网站| 免费在线观看亚洲国产| 欧美丝袜亚洲另类 | 亚洲自拍偷在线| av视频免费观看在线观看| 亚洲精品美女久久久久99蜜臀| 国产一区二区三区在线臀色熟女| 亚洲aⅴ乱码一区二区在线播放 | 中文字幕精品免费在线观看视频| 一本综合久久免费| 丝袜在线中文字幕| 亚洲av成人不卡在线观看播放网| 亚洲成人精品中文字幕电影| 黄色女人牲交| 亚洲精品美女久久av网站| 一区在线观看完整版| 国产成人啪精品午夜网站| 国产激情欧美一区二区| 久久精品国产综合久久久| 免费无遮挡裸体视频| 国产午夜福利久久久久久| 999久久久国产精品视频| 亚洲精品久久成人aⅴ小说| 免费高清视频大片| 国产成人精品在线电影| 久久亚洲真实| 麻豆国产av国片精品| 最新在线观看一区二区三区| 天堂影院成人在线观看| 中文字幕久久专区| 欧美绝顶高潮抽搐喷水| 侵犯人妻中文字幕一二三四区| 好男人电影高清在线观看| 久久婷婷人人爽人人干人人爱 | svipshipincom国产片| 人妻丰满熟妇av一区二区三区| 韩国av一区二区三区四区| 少妇 在线观看| 深夜精品福利| 国产不卡一卡二| 久久精品影院6| 日韩av在线大香蕉| 美女 人体艺术 gogo| 91在线观看av| 亚洲中文字幕日韩| 亚洲欧美日韩无卡精品| 人人妻人人爽人人添夜夜欢视频| 日韩视频一区二区在线观看| 国产精品久久久久久人妻精品电影| 午夜福利视频1000在线观看 | 欧美亚洲日本最大视频资源| 婷婷精品国产亚洲av在线| 1024香蕉在线观看| 99国产综合亚洲精品| 国产亚洲精品久久久久5区| 久久久久久久久久久久大奶| 给我免费播放毛片高清在线观看| 久久久国产成人精品二区| 国产激情久久老熟女| 脱女人内裤的视频| 9191精品国产免费久久| 亚洲欧美精品综合一区二区三区| 久久久久久免费高清国产稀缺| 午夜福利18| 99热只有精品国产| 日韩欧美免费精品| 18美女黄网站色大片免费观看| 国产精华一区二区三区| 国产国语露脸激情在线看| 国产97色在线日韩免费| 国产片内射在线| 看黄色毛片网站| 侵犯人妻中文字幕一二三四区| 韩国精品一区二区三区| 国产蜜桃级精品一区二区三区| 欧美激情极品国产一区二区三区| 性色av乱码一区二区三区2| 国产精品久久视频播放| 美女 人体艺术 gogo| 亚洲国产精品久久男人天堂| 亚洲av成人一区二区三| 黄色成人免费大全| 12—13女人毛片做爰片一| 波多野结衣av一区二区av| 国产成年人精品一区二区| 国产av精品麻豆| 成人手机av| 欧美一级毛片孕妇| 亚洲精品av麻豆狂野| 在线观看免费视频日本深夜| 国产在线观看jvid| 日韩欧美免费精品| 国产国语露脸激情在线看| 国产成人av激情在线播放| 国产精品一区二区精品视频观看| 亚洲成人精品中文字幕电影| 色综合站精品国产| 99香蕉大伊视频| 国产精品永久免费网站| 99精品在免费线老司机午夜| 97人妻天天添夜夜摸| 午夜影院日韩av| 欧美日韩中文字幕国产精品一区二区三区 | 国产欧美日韩精品亚洲av| 亚洲精品久久国产高清桃花| 亚洲av电影不卡..在线观看| 女同久久另类99精品国产91| 亚洲熟妇中文字幕五十中出| 十八禁人妻一区二区| 国产激情久久老熟女| 亚洲一区二区三区不卡视频| 亚洲全国av大片| 人人妻人人爽人人添夜夜欢视频| 色播亚洲综合网| 国产精品久久电影中文字幕| 国产熟女xx| 午夜福利高清视频| 国产精品,欧美在线| 夜夜爽天天搞| 精品日产1卡2卡| www日本在线高清视频| 人人妻人人澡欧美一区二区 | 中文字幕最新亚洲高清| 亚洲国产精品合色在线| 欧美日本亚洲视频在线播放| videosex国产| 国产精品秋霞免费鲁丝片| 久久 成人 亚洲| 老熟妇仑乱视频hdxx| 中文字幕人妻丝袜一区二区| 1024视频免费在线观看| 免费久久久久久久精品成人欧美视频| 美国免费a级毛片| 婷婷精品国产亚洲av在线| 国产成人精品在线电影| 国产极品粉嫩免费观看在线| 日日夜夜操网爽| 午夜精品国产一区二区电影| 亚洲熟妇熟女久久| 久9热在线精品视频| 欧美一级a爱片免费观看看 | 亚洲国产精品成人综合色| 欧美中文日本在线观看视频| 成人永久免费在线观看视频| 国产亚洲精品第一综合不卡| 黄色视频,在线免费观看| 欧美色欧美亚洲另类二区 | 国产一区二区激情短视频| 欧美人与性动交α欧美精品济南到| 欧美大码av| 99久久久亚洲精品蜜臀av| 亚洲精华国产精华精| 91九色精品人成在线观看| 国产精品,欧美在线| 禁无遮挡网站| 亚洲一区二区三区不卡视频| 深夜精品福利| 欧美日韩亚洲国产一区二区在线观看| 一二三四在线观看免费中文在| 亚洲国产精品sss在线观看| 精品免费久久久久久久清纯| 两人在一起打扑克的视频| 色播在线永久视频| 嫩草影视91久久| 亚洲狠狠婷婷综合久久图片| 性少妇av在线| 亚洲最大成人中文| 国产一区在线观看成人免费| 最近最新中文字幕大全免费视频| 亚洲avbb在线观看| 国产精品免费视频内射| 精品国产美女av久久久久小说| 1024视频免费在线观看| 欧美大码av| 丝袜在线中文字幕| 成人免费观看视频高清| 久99久视频精品免费| 国产亚洲精品综合一区在线观看 | 午夜福利视频1000在线观看 | 色综合亚洲欧美另类图片| 日韩大码丰满熟妇| 国产在线精品亚洲第一网站| 99久久久亚洲精品蜜臀av| 变态另类丝袜制服| 日韩一卡2卡3卡4卡2021年| 黑人操中国人逼视频| 国产在线观看jvid| a级毛片在线看网站| 亚洲欧美日韩高清在线视频| 欧美一区二区精品小视频在线| 亚洲中文av在线| 三级毛片av免费| 两性午夜刺激爽爽歪歪视频在线观看 | 日本精品一区二区三区蜜桃| av福利片在线| 91在线观看av| 国产一区在线观看成人免费| 三级毛片av免费| 香蕉久久夜色| 免费少妇av软件| 麻豆av在线久日| 中文字幕久久专区| 伊人久久大香线蕉亚洲五| 日韩视频一区二区在线观看| 亚洲欧美精品综合久久99| 亚洲欧美精品综合一区二区三区| 人人妻人人爽人人添夜夜欢视频| 免费在线观看黄色视频的| 美国免费a级毛片| 国产又色又爽无遮挡免费看| 欧美成人一区二区免费高清观看 | 免费在线观看日本一区| 日本免费一区二区三区高清不卡 | 可以在线观看毛片的网站| 国产野战对白在线观看| 色播亚洲综合网| 高潮久久久久久久久久久不卡| 人人妻人人澡欧美一区二区 | 亚洲一区二区三区不卡视频| 精品一品国产午夜福利视频| 男人的好看免费观看在线视频 | 成人免费观看视频高清| 久久久国产成人精品二区| 一级黄色大片毛片| 亚洲av成人不卡在线观看播放网| 黄片小视频在线播放| 免费在线观看黄色视频的| 香蕉久久夜色| 精品一品国产午夜福利视频| 高清毛片免费观看视频网站| 国产成人精品久久二区二区免费| 成人三级黄色视频| 一边摸一边抽搐一进一小说| 啦啦啦韩国在线观看视频| 亚洲自拍偷在线| 亚洲精品久久国产高清桃花| 真人做人爱边吃奶动态| 国产又色又爽无遮挡免费看| 日韩高清综合在线| 热re99久久国产66热| 国产精品香港三级国产av潘金莲| 久久午夜亚洲精品久久| 亚洲午夜理论影院| 在线观看免费视频网站a站| 后天国语完整版免费观看| 搞女人的毛片| 亚洲男人天堂网一区| 亚洲成a人片在线一区二区| 欧美成人午夜精品| 国产精品日韩av在线免费观看 | 亚洲国产欧美日韩在线播放| 91av网站免费观看| 高清在线国产一区| 国产精品影院久久| 又大又爽又粗| 可以在线观看毛片的网站| 久久国产精品男人的天堂亚洲| 别揉我奶头~嗯~啊~动态视频| 女人爽到高潮嗷嗷叫在线视频| 国产欧美日韩一区二区三| 老司机午夜十八禁免费视频| 侵犯人妻中文字幕一二三四区| 90打野战视频偷拍视频| 多毛熟女@视频| 国产视频一区二区在线看| 成人免费观看视频高清| 在线播放国产精品三级| 久久久久亚洲av毛片大全| 久久久久九九精品影院| 91在线观看av| 亚洲av美国av| 成人三级做爰电影| 丝袜在线中文字幕| 久久精品成人免费网站| 国产亚洲精品一区二区www| 免费在线观看完整版高清| 日本vs欧美在线观看视频| 日韩中文字幕欧美一区二区| 国产精品av久久久久免费| 欧美久久黑人一区二区| 一级a爱片免费观看的视频| 欧美日韩一级在线毛片| 亚洲专区国产一区二区| 99国产综合亚洲精品| 91老司机精品| 国产不卡一卡二| 久久这里只有精品19| 亚洲成人免费电影在线观看| 男女之事视频高清在线观看| 在线播放国产精品三级| 欧美日本亚洲视频在线播放| 午夜精品在线福利| 国产精品国产高清国产av| videosex国产| 九色国产91popny在线| 亚洲国产精品999在线| 国产精品秋霞免费鲁丝片| 亚洲三区欧美一区| 制服丝袜大香蕉在线| 夜夜躁狠狠躁天天躁| 国产av一区二区精品久久| 国产区一区二久久| 一进一出抽搐gif免费好疼| 在线观看午夜福利视频| 波多野结衣巨乳人妻| 一进一出好大好爽视频| 久久精品国产99精品国产亚洲性色 | 国产熟女xx| 欧美激情极品国产一区二区三区| 高清在线国产一区| 亚洲片人在线观看| 久久久久久亚洲精品国产蜜桃av| 日本免费一区二区三区高清不卡 | 免费观看精品视频网站| 久久精品aⅴ一区二区三区四区| 夜夜夜夜夜久久久久| 午夜视频精品福利| 好男人在线观看高清免费视频 | 国产精品亚洲美女久久久| 国产一区二区激情短视频| 国产日韩一区二区三区精品不卡| 久久狼人影院| 女人精品久久久久毛片| 女人被躁到高潮嗷嗷叫费观| 亚洲精品久久国产高清桃花| 麻豆一二三区av精品| 日本精品一区二区三区蜜桃| 高清在线国产一区| 国产不卡一卡二| 亚洲成国产人片在线观看| 亚洲精品在线观看二区| 三级毛片av免费| 免费在线观看影片大全网站| 这个男人来自地球电影免费观看| 欧美 亚洲 国产 日韩一| 国产精品九九99| 天堂影院成人在线观看| 亚洲片人在线观看| 国产亚洲欧美精品永久| 成人三级黄色视频| 狂野欧美激情性xxxx| av中文乱码字幕在线| 99久久国产精品久久久| 身体一侧抽搐| 两性夫妻黄色片| 一区二区三区激情视频| 在线观看午夜福利视频| 国产乱人伦免费视频| 亚洲aⅴ乱码一区二区在线播放 | 美女高潮到喷水免费观看| 两个人看的免费小视频| 亚洲,欧美精品.| 九色国产91popny在线| 亚洲精品一区av在线观看| 午夜精品久久久久久毛片777| 亚洲成人久久性| 99精品久久久久人妻精品| 一二三四在线观看免费中文在| 90打野战视频偷拍视频| 色综合婷婷激情| 国产一区二区三区综合在线观看| 国产成人av激情在线播放| 欧美黑人欧美精品刺激| 亚洲欧美日韩高清在线视频| 国产一区二区三区综合在线观看| 国产成人av激情在线播放| 一二三四社区在线视频社区8| 国产精品一区二区免费欧美| 88av欧美| 日韩大尺度精品在线看网址 | 国产亚洲欧美精品永久| 久久亚洲精品不卡| 亚洲国产中文字幕在线视频| 丰满的人妻完整版| 欧美久久黑人一区二区| 亚洲精品美女久久av网站| 99国产精品99久久久久| 露出奶头的视频| 欧美日本视频| 免费高清在线观看日韩| 成人永久免费在线观看视频| 男女之事视频高清在线观看| 在线观看舔阴道视频| 琪琪午夜伦伦电影理论片6080| 亚洲午夜理论影院| 午夜亚洲福利在线播放| 久久久久国产精品人妻aⅴ院| 电影成人av| 九色亚洲精品在线播放| 午夜福利在线观看吧| 妹子高潮喷水视频| 免费在线观看黄色视频的| 亚洲精品在线观看二区| 十分钟在线观看高清视频www| 中文字幕人成人乱码亚洲影| 两个人看的免费小视频|