渤海灣盆地臨清坳陷的數值模擬與異常沉降

王子揚邱亮王昆王彥舉蘇子俊向禹艮

渤海灣盆地臨清坳陷的數值模擬與異常沉降

王子揚1,邱亮2*,王昆1,王彥舉1,蘇子俊1,向禹艮1

1. 山東科技大學地球科學與工程學院,山東 青島 266590 2. 中國地質大學（北京）地球科學與資源學院, 北京 100083

為探索渤海灣盆地南部的臨清坳陷機理，本文利用回剝法和改進的有限元應變速率反演技術，基于3個剖面14口實際井和合成井的數據，獲取了拉伸因子、應變速率、計算沉降量等數據，揭示了基于觀測沉降量和計算沉降量井差值的異常沉降現象：中新世時期負差值的減速沉降，和上新世至第四紀時期正差值的加速沉降。平衡剖面恢復和區域熱流值計算則指示了中新世減速沉降的機制主要為熱源的持續流入，抑制了巖石圈的熱收縮；中新世以來的加速沉降則與后裂谷期的正斷層活動密切相關。

臨清坳陷; 地面沉降; 數值模擬

裂后沉降是指發生在斷裂構造控制的同裂谷階段之后的熱收縮沉降。前人研究認為，巖石圈伸展主導的裂谷沉降后，熱收縮作用主導了裂后沉降階段[1,2]。然而，最近的研究表明，野外觀測獲取的實際裂后沉降與麥肯齊熱收縮沉降模型的預測結果存在差異，這種差異被稱為加速沉降[3]。加速沉降在被動大陸邊緣盆地和大型裂谷盆地均有發現，如歐洲的羅克爾盆地和法羅設得蘭盆地，南海的馬來盆地、珠江口盆地，及東亞的松遼盆地和渤海灣盆地等陸相裂谷盆地中也有發現[4-6]。渤海灣盆地的研究表明，加速沉降現象集中于后期，并沒有貫穿整個裂后沉降階段[7]。那么，裂后階段是否存在減速沉降的現象，加速沉降和減速沉降及其轉變的的地球動力學機制是什么。

臨清坳陷位于太平洋板塊在東亞下方地幔過渡區的滯留板片西緣與重力異常交界處，是研究深部地質過程與上層地殼盆地演化之間相互作用的理想區域。本文選取位于渤海灣盆地西南部的臨清坳陷為研究對象，重點研究了后裂谷階段加速沉降和減速沉降的產生及其機制，利用包含14口真實或合成井的3個地震截面進行建模，從中獲得拉伸因子、應變速率、預測沉降和古熱流等數據，揭示臨清坳陷的后裂谷期構造演化。

1 研究概況

臨清坳陷是一個中、新生代復合型斷陷沉降區，在構造區劃上屬于渤海灣盆地二級負向構造單元，位于渤海灣盆地的西南收斂端，呈北北東向展布，北接黃驊坳陷、滄縣隆起，南至內黃隆起，西鄰太行山隆起，東至魯西南隆起[7,8]。臨清坳陷構造演化經歷了古近紀始新世～漸新世快速沉降與堆積的斷陷階段以及新近紀中新世一第四紀連續沉積的區域性坳陷階段[9]。區內主要發育中-新生代蓋層，新生代地層是中國重要的油氣富集層位，主要由古近紀沙河街組、東營組，新近紀館陶組和明化鎮組及第四紀地層組成[10]。

2 材料與方法

2.1 數據來源及參數

本次研究中用于沉降分析的數據來自臨清坳陷的3個地震剖面（從南到北依次為1、2和3）中的3口實井（標記為s）和11口合成井（標記為h）。基于地震剖面和實際井的分析，對模擬井的深度和巖性進行了推算。通過對收集到的地層巖性、厚度和孔隙度等基礎數據進行分解和回剝，重建觀測到的構造沉降厚度。地殼和巖石圈厚度及其他參數均來自前人研究（表1）；為了簡化計算，我們將古水深定為零。

表 1 參數及取值[8-10]

2.2 實際構造沉降量

用于計算構造沉降的軟件是斯倫貝謝公司開發的基于回剝法的PetroMod 2012.1。回剝法是一種基于Airy均衡模型的反演方法，包括以下幾種步驟：分解、古水深校正、沉積界面溫度校正、剝蝕校正、海面起伏和構造沉降提取[11]。在回剝過程中還應考慮跨越780萬年，從距今23.8百萬年至距今16.0百萬年（23.8～16.0 Ma）的剝蝕作用（即東營運動）。為了便于比較不同剖面的應變速率和斷層活躍率，假設東營運動期間沉積和剝蝕速率為零。

表 2 回剝法中使用的巖石力學參數[11]

2.3 模擬構造沉降量

確定預測沉降需要兩個步驟。第一步，通過輸入背剝曲線進行改良的一維應變速率反演，獲得應變總量和應變速率。第二步，采用改良的有限可拓模型，以為輸入參數，計算預測沉降。

2.3.1 一維應變速率反演前人研究提出了一種基于構造沉降史計算應變總量和應變速率的新算法（式1）[12,13]。

式中，()是應變速率對時間的函數，由于方程1難以得到解析解，可采用鮑威爾算法進行數值迭代計算。然而，鮑威爾算法作為一種共軛收斂加速方法，需要大量數據，而沉降曲線無法提供這些數據。此外，沉降曲線的線性依賴性可能導致搜索退化，導致不收斂問題[14]。本文對傳統方法進行了改進，采用牛頓二階插值迭代法求得數值解；將積分轉化為和，就可以求得的值（式2）。

式中，是中點回歸參數，G是中點應變速率，,和是鮑威爾和諧矩陣的初始取值，T是伸展起始時間，其它字符見表1。

2.3.2 改進的有限擴展模型有限擴展模型是Jarvis等提出的經典均勻伸展模型的改進版本。有限擴展模型不僅考慮了均勻伸展模型中包含的熱平流項，還在熱流方程中加入了熱擴散函數，以反映巖石圈結構中溫度隨時間變化的情況[15]。但是，來自地殼的放射性生熱在盆地發育中具有重要作用，在模擬盆地演化過程時，必須考慮放射性生熱。通過確定大陸地溫并使其等同于基礎熱流，可以用以深度為變量的負指數模型來模擬地殼放射性生熱。

式中，等式右側第一項代表水平熱流值，第二項代表熱擴散函數，第三項代表了放射性熱流值。

由于應變總量是沉降量的函數，因此對上式進行變形，得到的表達式：

式中，y是利用回剝法計算得出的實際構造沉降量，其它參數含義見表1。

3 構造沉降

3.1 凹陷南部剖面A

剖面A橫切整個坳陷盆地，西部受分別西傾和東傾的斷層共同控制，形成不對稱地塹，東部受控盆斷層控制，形成地塹構造，中央隆起充當兩側斷層的共同下盤。剖面西緣斷層活動不占優勢，古生代地層與斷層帶重疊，形成西部斜坡帶（圖1a）。

圖 1 橫切坳陷南部的構造剖面A

（a）根據地震剖面解釋的構造剖面，井sa1～sa6的位置；（b～g）sa1～sa6井的伸展總量、應變速率和觀測沉降（藍線）與理論沉降（橙色線）曲線。

剖面A的裂谷期地層為古近系沙河街組和東營組組成，平均應變總量約為1.75。裂谷期的應變速率變化指示了兩幕次的裂谷活動時期，對應形成了沙河街組三段～二段（Es3-Es2）和東營組（Ed）地層。Es3-2活動時期以Sa6井為中心，最大應變速率達～2×10-15s-1，值最大為1.80。Ed時期以Ma3井為中心，最大應變率是～1.3×10-15s-1，值最大為1.70。根據不同的應變率可以觀察到不同井構造沉降的變化。Es3-2應變率的遞增順序是Ma6、Ma5、Ma3、Ma1、Ma2、Sa4。Ed應變率的遞增順序為Ma3、Ma2、Ma1、Ma6、Sa4、Ma5（圖1b-g），兩幕活動的應變速率表現此消彼長的數值關系。

東營組和館陶組間的平行不整合分隔了裂谷期與后裂谷期。中新世的后裂谷期應變率為零，值保持不變，與McKenzie的理論模型一致。然而，上新世和第四紀的應變率和值明顯增加，表明該時期出現快速沉降（圖1b-g）。根據沉降的觀測值（藍線）與預測值（粉線）的區別，可以確定中新世以來的理論預測值與實際觀測值之間存在差異（圖1b-g）。沉降差在中新世期為-96～328 m，在上新世期間為+237～407 m。中新世觀測沉降值的傾角平緩，而上新世-第四紀則較為陡峭。

3.2 凹陷中部剖面B

凹陷中部是同裂谷期斷層活動最強烈的區域。中部巖陷的結構是典型的半地塹，以西傾的鏟狀蘭聊斷層為主。同裂谷期的兩幕活動模式（Es3-2和Ed）由應變速率圖清晰可見（圖2b-f），平均應變率和值分別為2.2×10-15s-1和1.78，Es3-2時期斷層活動較強烈；Ed期以弱斷層為特征，應變率只有0.5×10-15s-1。Es3-2和Ed的應變速率呈現正相關，變化模式與剖面A相反，整個凹陷中部的沉降受蘭聊斷層的完全控制。

圖 2 橫切坳陷中部的構造剖面B

后裂谷期的應變速率在中新世均為零，而在上新世-第四紀期間明顯增加，平均值為～0.4×10-15s-1；值在中新世保持不變，在上新世-第四紀明顯增加，平均值為1.82。中新世觀測沉降直線的傾角平緩，而上新世-第四紀更為陡峭，藍、粉線之間的區域代表觀測沉降量和計算沉降量之間的差值，中新世兩者的沉降差為-266～343 m，上新世則為+123～247 m。

3.3 凹陷北部剖面C

北部凹陷是由蘭聊斷層控制的復合半地塹，結合上盤發育的對立的次級斷層。這些次級結構以旋轉的平面正斷層為標志，形成多米諾骨牌斷層系統。單峰應變率模式表明，該地區出現了獨特的Es3-2斷層期，與南部和中部的裂谷發育模式明顯不同。Ed時期多為弱斷層，特點是應變率低，約為0.13×10-15s-1。

圖 3 橫切坳陷北部的構造剖面C

后裂谷階段的應變速率模式和應變總量模式與南部和中部凹陷類似。中新世時期應變速率為零應變總量保持不變，在上新世-第四紀期間增至平均值約為0.4×10-15s-1，應變總量持續增長至第四紀。中新世兩線之差為-185～218 m，在上新世-第四紀為+266～291 m（圖3b-d）。通過北部凹陷地層的厚度和應變率可知，在新近紀期間，沉積中心遷移到了西部洼陷。

4 討論

4.1 加速沉降與減速沉降

根據均勻伸展模型，巖石圈在同裂谷末期斷層活動逐漸停止，后裂谷初期熱收縮沉降啟動。應變速率ε和應變總量是評價同裂谷期斷層伸展強度的主要定量指標。當應變速率為零，應變總量保持不變，即表示正斷層活動的終止。由于總熱流隨時間呈現負指數型變化，熱沉降曲線因此表現為負指數模型[16-18]。

盡管研究區的中新世后裂谷期應變速率為零、應變總量值保持穩定，符合均勻伸展模型的預測，但觀測沉降量與計算沉降量存在較大差異，觀測沉降量明顯小于計算沉降量，該負差值范圍為-343～96 m，這指示了一種非常規的減速沉降模式。由于在23.8～16 Ma階段整個渤海灣盆地進入了一個剝蝕夷平期，該減速沉降起始于16 Ma。新近紀-第四紀階段，應變速率>0、應變總量不斷增長的情況下，觀測沉降量與計算沉降量之間的正差值范圍是123～407 m，這表示上新世-第四紀期間出現了加速沉降。渤海灣盆地的濟陽凹陷和渤中凹陷均已觀測到加速沉降現象，減速沉降則是首次被發現。

在C剖面的所有井與剖面A的Ma1到Ma4井中，當前（0 Ma時）的觀測沉降量與計算沉降量之間均為正差值，其余剖面的井均為負差值。該現象具有兩種不同的分布類型，正差異的特點是藍線（觀測沉降值）與粉線（預測沉降值）出現交叉點；負差異則以不相交為特征。藍、粉線相交的地質意義是，后期的加速沉降補充了早期減速沉降造成的厚度不足，剩余量為正差值表明后裂谷期總沉降高于均勻伸展模型預測的熱沉降量。兩線不相交則表明，后期加速沉降量不足以彌補早期減速沉降的缺失，產生了50～150 m的負差異。但按照當前藍線、粉線的延長趨勢，未來加速沉降可能能夠補償減速沉降造成的沉降量缺失。

4.2 裂后期正斷層活動與加速沉降

根據均勻伸展模型，正斷層活動在后裂谷期停止，但渤海灣盆地的渤中凹陷和冀中凹陷均發現了裂后期的正斷層活動[12,13]。因此可以證明，正斷層的重啟造成了渤海灣盆地的裂后期加速沉降。

在上新世-第四紀期間，臨清坳陷也出現了裂后期的正斷層活動。如圖4所示，基底斷層快速復活的特點是上盤的沉積層比下盤更厚。斷層的活動速率表明斷層活動相對強烈，這一時期的活動速率甚至大于Es4、Es2等裂谷期斷層活動相對平靜階段。因此，后裂谷期的正斷層活動與加速沉降密切相關。后裂谷期正斷層活動的機制較為復雜，人們提出了幾種構造模型來解釋這種現象，如郯廬斷層的右旋運動，新裂谷活動期和華北克拉通破壞的重新啟動等，對該問題的成因還有爭議。

平衡剖面和斷層活動率表明，中新世期缺少正斷層活動。同時，正斷層活動的結果只能導致加速沉降發生，與中新世的減速沉降因果關系較小。因此，中新世的減速沉降應歸因于其他構造因素。

圖 4 后裂谷期平衡地質剖面

4.3 古熱流與減速沉降

巖石圈減薄后的熱收縮是造成裂后沉降的主要機制，由于熱流減少，熱收縮沉降的速度隨時間呈指數下降。因此，熱流是裂后沉降的直接控制因素，加速和減速沉降的出現可能會反映熱流的轉變。為了探究沉降變化與古熱流之間的關系，本文設計了基于熱流模型預測值和觀測值的綜合分析。以平均應變總量=1.75的剖面A為例，使用改進的模擬熱流方程對剖面A的預測平均熱流和沉降進行估算。

由于后裂谷期沉降是同裂谷期沉降的函數，可以據此調整模型獲得熱流的觀測值：給從應變率反演中獲得的值增加一個Δ值（如0.05），計算出的沉降量也相應的增加。當結果與觀測的裂后沉降相吻合時，可認為此時的熱流與沉降的觀測值相對應。該值以式5的左邊表示。

其中式5等號左側是與觀測的沉降相對應的熱流。六口井的數據處理完成后，可采用線性指數擬合的方法得到與沉降觀測值相對應的熱流曲線（圖5a）。臨清坳陷的古熱流模型得出的數值為61～94 mW/m2，與本研究估計的熱流值一致。

a 比例地質年代尺度為20 Ma，比例熱流值以裂谷發育前的初始熱流為1

b 示模型計算的理論熱流值，紅線表示根據實際沉降量反演的實際熱流值

根據剖面A的模擬結果，在中新世期間，熱流呈線性下降，隨后在上新世-第四紀期間表現為曲線下降的趨勢。在整個裂后階段，熱流的觀測值比預測值要大（圖5b）。在7.7 Ma前后出現了兩條曲線的最大差異（約13.67）。在加速沉降期（18.5～5.3 Ma），觀測和預測的熱流差在3.01～13.67 m之間，而該值在減速沉降期約為10.44～12.67 mm2。

在中新世，熱收縮損失的熱量將由其他熱源繼續提供，盡管不足以阻止熱流的減少，但可以減緩熱流減少的速度，抑制巖石圈的熱收縮。因此，熱流減速過慢導致熱收縮不充分，沉降量并不明顯。在上新世早期，熱流冷卻模型從線性過渡到曲線，其他熱源消失，熱收縮以McKenzie的負指數模型重新啟動。熱流值約為78.5 m2，比預測值高16%（圖5b），其產生的熱收縮能量更強，導致明化鎮組出現了超充分沉降。因此，加速、減速沉降的出現與熱流變化有明顯的相關性。

5 結論

本文針對渤海灣盆地臨清坳陷的后裂谷階段，利用包含14口真實或合成井的3個地震截面，采用回剝法和改進的有限元應變速率反演技術進行數值模擬，獲取了拉伸因子、應變速率、計算沉降量和古熱流等的數據。后裂谷期觀測沉降量與計算沉降量的對比揭示了中新世的減速沉降，沉降差值為?96～343 m，和上新世-第四紀的加速沉降，沉降差值為+123～407 m。通過平衡剖面恢復和區域后裂谷期熱流值計算，認為中新世后裂谷期減速沉降主要歸因于熱源的持續流入，抑制了巖石圈的熱收縮；中新世加速沉降則與后裂谷期的正斷層活動密切相關。

[1] McKenzie D. Some remarks on the development of sedimentary basins [J]. Earth Planetary Science Letters, 1978,40(1):25-32

[2] McKenzie D, Jackson J, Priestley K. Thermal structure of oceanic and continental lithosphere [J]. Earth Planetary Science Letters, 2005,233(3-4):337-349

[3] McKenzie D, Watts A, Parsons B,. Planform of mantle convection beneath the Pacific Ocean [J]. Nature, 1980,288:442-446

[4] Madon M, Watts A. Gravity anomalies, subsidence history and the tectonic evolution of the Malay and Penyu Basins ( offshore Peninsular Malaysia ) [J]. Basin Research, 1998,10:375-392

[5] 漆家福.渤海灣新生代盆地的兩種構造系統及其成因解釋[J].中國地質,2004,31(1):15-22

[6] Qi J, Yang Q. Cenozoic structural deformation and dynamic processes of the Bohai Bay Basin province, China [J]. Marine and Petroleum Geology, 2010,27(4):757-771

[7] Qiu NS, Xu W, Zuo YH,. Meso-Cenozoic thermal regime in the Bohai Bay Basin, eastern North China Craton [J]. International Geology Review, 2015,57(3):271-289

[8] Qiu NS, Zuo YH, Chang J,. Geothermal evidence of Meso-Cenozoic lithosphere thinning in the Jiyang sub-basin, Bohai Bay Basin, eastern North China Craton [J]. Gondwana Research, 2014,26:1079-1092

[9] 任鳳樓,柳忠泉,邱連貴,等.渤海灣盆地新生代各坳陷沉降的時空差異性[J].地質科學,2008,43(3):546-557,575

[10] Ren J, Tamaki K, Li S,. Late Mesozoic and Cenozoic rifting and its dynamic setting in Eastern China and adjacent areas [J]. Tectonophysics, 2002,344(3-4):175-205

[11] 史卜慶,吳智平,王紀祥,等.渤海灣盆地東營運動的特征及成因分析[J].石油實驗地質,1999,21(3):196-200

[12] White N. Recovery of strain rate variation from inversion of subsidence data [J]. Nature, 1993,366(6454):499-452

[13] White N. An inverse method for determining lithospheric strain rate variation on geological timescales [J]. Earth Planetary Science Letters, 1994,122(3-4):351-371

[14] Silva J. Influence diagnostics and estimation algorithms for powell's SCLS [J]. Journal Bus. Economic Statistics, 2001,19(1):55-62

[15] 肖國強,楊吉龍,趙長榮,等.天津濱海地區G2孔磁性地層年代及其構造指示[J].地質通報,2014(10):1642-1650

[16] Song B, Chen L, Zhang J,A matlab program for 1d strain rate inversion [J]. Computer Geoscience, 2010,36(1):16-23

[17] 萬天豐.論構造事件的節律性[J].地學前緣,1997,4(3-4):257-263

[18] Wheeler P, White N. Quest for dynamic topography: observations from Southeast Asia [J]. Geology, 2000,28(11):963-966

Numerical Simulation and Abnormal Sedimentation of Linqing Depression on Bohai Bay Basin

WANG Zi-yang1, QIU Liang2*, WANG Kun1, WANG Yan-ju1, SU Zi-jun1, XIANG Yu-gen1

1.266590,2.100083,

In order to explore the mechanism of the Linqing depression in the southern Bohai Bay Basin, this paper uses the back-stripping method and the improved finite element strain rate inversion technique. Based on the data of 14 actual and synthetic Wells in 3 profiles, the data of tensile factor, strain rate and calculated settlement are obtained, and the abnormal settlement phenomenon based on the well difference between observed settlement and calculated settlement is revealed: The deceleration of negative difference in Miocene period, and the acceleration of positive difference in Pliocene to Quaternary period. The restoration of the equilibrium profile and the calculation of regional heat flow indicate that the mechanism of the slow subsidence in Miocene is mainly the continuous inflow of heat source, which inhibits the thermal contraction of the lithosphere. The accelerated subsidence since Miocene is closely related to the normal fault activity in the post-rift period.

Lingqing depression; ground subsidence; numerical simulation

P642.26

A

1000-2324(2022)06-0956-07

2022-05-02

2022-07-21

國家自然科學基金委青年基金(41702207)

王子揚(2002-),男,本科生,主要從事工程地質及地面沉降機理探究. E-mail:2652382461@qq.com

Author for correpondence. E-mail:qiuliang2011@gmail.com

10.3969/j.issn.1000-2324.2022.06.023