FIRST TIME DISCOVERY OF GUANGFULIN REMAINS AND LIVING PALEOENVIRONMENT OF ANCESTORS IN SHANGHAI REGION
-
摘要: 上海松江广富林遗址1999年重新发掘,在发掘中首次发现了类似豫东王油坊类型的文化——广富林遗存文化。本次采集了广富林遗址05SGIT1240探方样品,进行孢粉、藻类研究。根据孢粉、藻类成分特征,划分出了3个孢粉组合带,恢复了当时的古植被和古环境,并对广富林遗存先人生存环境、农耕发展作了探析,为环太湖地区文化发展及文化变迁提供了十分珍贵的新资料。Abstract: Guangfulin relics was re-explored in 1999 in Songjiang of Shanghai. Guangfulin remains culture was found for the first time there,and there are some similarities to Wangyoufang type cultural historical remains from eastern Henan Province. Samples were taken from the west wall in unit 05SGIT1240 of Guangfulin relics for sporopollen and algae research.According to the characteristics of the sporopollen and algae assemblages, three sporopollen zones were distinguished,and paleovegetation and paleoenvironment were reconstructed. And then, we analyzed the living environment and development of ancestors' cultivation. The discovery of Guangfulin remains provides precious new data for culture development and processes around the Taihu Lake region.
-
Keywords:
- Guangfulin remains /
- first time discovery /
- ancestor /
- paleoenvironment
-
天然气水合物是由水分子与气体分子组成的笼状结构化合物[1]。在自然界中，水合物中的气体分子主要是甲烷，因此是潜在的能源，同时也成为潜在的环境灾害因素。勘探结果表明，天然气水合物主要存在于海洋沉积物以及陆地冻土地区，其全球储量巨大，约为2×1016 m3[2]。
低温高压是天然气水合物稳定存在的必要条件，水合物的开采原理就是破坏水合物相平衡状态，使其分解为水和气体[3]。目前，常规的开采方法主要有3种：降压法[4]、注热法[5]以及注入抑制剂[6]。现场试采、实验以及数值模拟结果表明，降压法是最具潜力的开采方法，而其他开采方法适合作为辅助手段来提高水合物的开采效率[7]。
开采过程会破坏天然气水合物的相平衡状态，使储层划分为不同的区域：分解区域、未分解区域，在两者之间存在一个特殊的过渡区域，即分解前缘。以水合物分解前缘为界的分解区域与未分解区域之间的流体性质、储层地质条件以及地层温度压力规律等都不相同[8]。在开采过程中水合物分解前缘移动规律与水合物开采动态有着紧密关联，分解前缘能直接反映了水合物开采特征。另外，由于水合物在储层中起到了一定的胶结支撑作用，随着分解前缘的移动，地层稳定性降低，可能引发地质灾害[9]。因此，在水合物开采过程中，研究水合物分解前缘移动规律具有重要指示意义。
目前，对于天然气水合物分解前缘移动问题的理论研究多基于Stefan边界理论，即将水合物分解过程类比冰消融过程，是一个伴随相变的传热过程。Makogon首次借鉴Stefan问题，计算得到了降压分解水合物过程中压力分布的自相似解[10]。Verigin利用Stefan移动边界问题，建立了一维半无限大水合物藏降压开采模型，模型考虑分解前缘两侧的气相流动以及分解前缘处气体质量守恒，根据模型与Stefan问题的相似性对模型进行线性化自相似求解[11]。Ji等在Verigin模型基础上，建立水合物藏降压开采数学模型，假设水相静止，考虑了气相流动以及温度的变化，模型认为对流传热的作用比热传导强，分解区和水合物区的能量守恒方程中考虑了热对流以及节流和气体绝热效应，分解前缘处没有考虑水合物分解吸热的作用，将模型方程线性化处理后，自相似求解了分解前缘随时间的移动[12]。喻西崇等借鉴Ji等提出的数学模型，利用自相似原理推导出分解前缘移动表达式与温度、压力分布表达式[13]。唐良广，李刚等将水合物分解过程看作移动边界问题，建立了水合物层温度分布的一维传热模型，模型考虑了分解区和水合物区的热传导以及分解前缘处的能量守恒，根据自相似求解得到不同时刻水合物藏温度分布以及分解前缘的位置[14]。张旭辉等建立二维热传导模型，研究水合物储层有热水管垂直穿过时水合物最大分解范围，采用分离变量法对模型进行求解，结果表明分解前缘的最大移动距离随温度的增大而增大[15]。刘乐乐建立水平一维降压-加热数学模型，将模型有限差分离散后数值计算得到分解相变阵面的位置与时间的平方根呈正比[16]。李明川等建立注热移动界面的三相一维传质模型，数值差分计算得到分解前缘移动速度前期较高，后逐渐降低[17]。另外，随着数值模拟器的成熟，Long和Tjok利用HydrateResSim模拟水合物藏降压开采分解前缘移动，结果表明分解前缘的平均速度随绝对渗透率的增大而增大[18]。郑如意等数值模拟研究模型边界条件、渗透率、初始水合物饱和度、总热导率、井筒加热温度和井底压力等对水合物分解前缘移动的影响，发现分解前缘移动速度与渗透率、总热导率、井筒加热温度和边界供热成正比，相反，增大水合物初始饱和度和井底压力会降低分解前缘移动速度。此外，分解前缘移动规律也会随着参数的变化而变化[19]。
总的来说，前人采用数学建模以及数值模拟等多种方法对水合物分解前缘的移动规律进行了大量研究。与模型数值解相比，解析计算模型求解方便，计算过程也有助于深入了解某些物理变化的重要性。但大部分解析计算依据自相似原理，以分解前缘的移动与时间的平方根呈线性关系为前提假设，由此得到的分解前缘移动只在平均速度上存在差异，而分解前缘位置随时间变化规律已经确定，这将影响分解前缘移动规律探究。此外，利用现有数值模拟器探究水合物分解前缘移动规律时，要综合网格压力、相平衡压力以及水合物饱和度来判断分解前缘移动位置，这个过程相对复杂繁琐，并且单一储层模型的探究不具有普适性。
据此，本文建立了水合物降压分解一维三相数学模型，区别于传统自相似求解假设前提（即分解前缘的移动与时间的平方根呈线性关系），利用量级分析与偏微分方程无维化转化方法，解析计算探究了分解前缘随时间的移动规律，并由分解前缘移动计算相关产气量，对水合物开采动态进行了简单快速评估。
1. 概念模型
模型为Class3水合物藏。储层顶底板渗透率低，压力传递较慢，将顶部层和底部层认为是定温定压边界，允许发生热量和流水交换。
水合物分解是一个吸热的相变过程，会导致地层温度降低，此时，周围环境就会向水合物层传递热量。这种热补偿除了维持水合物分解所需热量外，还可以弥补吸热反应导致的地层温度下降。对于具有一定厚度的水合物层，可忽略顶底层围岩传入热量的影响，水合物分解相变前缘所吸收的热量主要来自于单位体积储层内能和未分解区热传导[20]。
随着时间的推移，水合物分解由近井区域向外扩散。假设水合物分解不是发生在整个储层内，而是发生在一定的狭窄区域内，可以将该区域视为一个表面，即所谓的水合物分解前缘[12]。它将储层分为两个部分，即水合物分解区与水合物未分解区。水合物分解区域自由气体在压力梯度的驱动下向井内流动，而分解前缘则向相反的方向移动。
根据前述假设，水合物分解过程可简化为分解前缘随时间向外移动的过程，而分解前缘就是天然气水合物发生分解的临界面，其厚度忽略不计，分解前缘处的地层压力即为该地层温度下的水合物相平衡压力[21]。在水合物层发生降压分解时，储层尺度内可采用水合物平衡分解模型，分解前缘界面（S（t）位置）把水合物层分为两个区域（图1）：已分解气水区（r0＜x＜S（t））与未分解水合物区（S（t）＜x＜∞）。
![]()
图 1 天然气水合物开采分解前缘迁移概念模型Figure 1. Conceptual model of migration of decomposition front of gas hydrate根据概念模型的简化，本文进一步假设为：
（1）模型考虑三相（水合物、甲烷、水）两组分（甲烷、水），不考虑甲烷和水合物的溶解；
（2）由于气体与水之间的压力差不大，忽略毛细管压力的影响；
（3）在水合物开采中，扩散作用贡献小于对流作用，忽略气水扩散对水合物分解的影响；
（4）水合物分解所需能量主要包括分解前缘所在区域的单位体积储层内能以及未分解区域传导热。
2. 数学模型
在数学模型中，将水合物储层看作x方向的一维流体场。描述多孔介质中水合物分解过程的主要方程包括水合物分解区域和分解前缘处的质量守恒方程、能量守恒方程。其中，分解前缘移动相关方程参考Stefan模型对冰水自由边界的描述[22]，其他方程则类似于Tsypkin[23-24]、Ahmadi[25]等提出的模型方程。
2.1 模型主要方程
水合物分解区域质量守恒方程：
$$\begin{split} &\dfrac{{\partial \left( {\phi {{{S}}_{\rm{i}}}{{{\rho }}_{\rm{i}}}} \right)}}{{\partial {\rm{t}}}}{\rm{ + }}\dfrac{{\partial \left( {{{{\rho }}_{\rm{i}}}{{{V}}_{\rm{i}}}} \right)}}{{\partial {\rm{x}}}}{\rm{ = 0}}\\ &{V_i} = - \dfrac{{K{K_{{\rm {ri}}}}}}{{{\mu _{\rm i}}}}\dfrac{{\partial {P_{\rm i}}}}{{\partial x}}{\rm{(i = w,}}\;{\rm{g)}} \end{split}$$ (1) 式中：
$ {\rm{\phi }} $ 是水合物层孔隙度，下标w，g分别代表水相与气相，S为水合物饱和度，K为水合物层绝对渗透率，$ {K}_{\rm {rg}} $ 与$ {K}_{\rm {rw}} $ 为气相与水相的相对渗透率，作为水相饱和度的函数。当分解区与未分解区的压力梯度都比较小时，水合物分解区域与未分解区域的能量守恒方程简化为（体现水合物分解区域与未分解区域温度变化关系）：
$$\frac{{\partial {{T}}}}{{\partial {{t}}}} = \frac{{{\lambda _{{\rm{1}},{\rm{2}}}}}}{{\rho {c_{{\rm{1}},{\rm{2}}}}}}\Delta T = {a_{{\rm{1}},{\rm{2}}}}\Delta T$$ (2) 式中：下标1，2表示分解区与未分解区，λ是导热系数，c是比热容，
$ \Delta $ 为拉普拉斯算子。水合物分解前缘处质量守恒方程为：
$$\begin{split} &\phi \left( {{\rho _{\rm h}}{S_{\rm h}}\varepsilon - {\rho _{\rm g}}{S_{\rm g}}} \right)\dfrac{{{\rm d}S\left( t \right)}}{{{\rm d}t}} = {\rho _{\rm g}}{V_{\rm g}}\\ &\phi {\rho _{\rm h}}{S_{\rm h}}(1 - \varepsilon ) = {\rho _{\rm w}}{S_{\rm w}} \end{split}$$ (3) 式中：
$ {\rm{\varepsilon }} $ 是单位体积水合物中气体所占的体积分数，S（t）为分解前缘位置。水合物分解前缘处能量守恒方程为（体现水合物分解前缘上的温度变化）：
$$\phi {\rho _{\rm h}}{S_{\rm h}}{q_{\rm h}}\frac{{{\rm d}S\left( t \right)}}{{{\rm d}t}} = {\lambda _1}{\rm {grad}}{T_1} - {\lambda _2}{\rm {grad}}{T_2}$$ (4) 式中：
$ {{q}}_{\rm{h}} $ 为单位质量水合物分解所需热量。2.2 模型辅助方程
分解前缘压力为水合物分解为气水的相平衡压力[26]：
$$\begin{split} &{\log _{10}}{P_{\rm D}} = a\left( {{T_{\rm D}} - {T_{\rm i}}} \right) + b{\left( {{T_D} - {T_0}} \right)^2} + c\\ &{{a}} = 0.034\;2/{{K}},\;{{b}} = 0.000\;5/{{K}},\\ &{{c}} = 6.480\;4,\;{{{T}}_0} = 273.15{{K}} \end{split}$$ (5) 式中：
$ {T}_{\rm D} $ 与$ {P}_{\rm D} $ 分别代表相平衡温度与压力，$ {T}_{i} $ 为参考温度，a、b，c是和水合物成分有关的经验常数，由Makogon平衡压力—温度数据得到。计算热传导系数方程：
$$ \begin{split} & {\mathrm{\lambda }}_{1}={\phi }{S}_{{\rm h}}{\mathrm{\lambda }}_{{\rm h}}+{\phi }{S}_{{\rm w}}{\mathrm{\lambda }}_{{\rm w}}+(1-{\phi }){\mathrm{\lambda }}_{{\rm s}}\\ &{\mathrm{\lambda }}_{2}={\phi }{(1-S}_{{\rm w}}{)\mathrm{\lambda }}_{{\rm g}}+{\phi }{S}_{{\rm w}}{\mathrm{\lambda }}_{{\rm w}}+(1-{\phi }){\mathrm{\lambda }}_{{\rm s}}\\ & {\mathrm{\rho }\mathrm{c}}_{1}={\phi }{{\rho }_{{\rm h}}S}_{{\rm h}}{\mathrm{c}}_{{\rm h}}+{\phi }{S}_{{\rm w}}{\rho }_{{\rm w}}{\mathrm{c}}_{{\rm w}}+(1-{\phi }){\rho }_{{\rm s}}{\mathrm{c}}_{{\rm s}}\\ &{\mathrm{\rho }\mathrm{c}}_{2}={\phi }{{\rho }_{{\rm g}}S}_{{\rm g}}{\mathrm{c}}_{{\rm g}}+{\phi }{S}_{{\rm w}}{\rho }_{{\rm w}}{\mathrm{c}}_{{\rm w}}+(1-{\phi }){\rho }_{{\rm s}}{\mathrm{c}}_{{\rm s}} \end{split} $$ (6) 2.3 模型初边值条件
模型初始条件与边界条件如表1所示。
表 1 初始条件与边界条件Table 1. Initial conditions and boundary conditions初边值条件 值 初始地层压力(${P}\left({x},0\right))$ $ {P}_{0} $ 初始地层温度(${T}\left({x},0\right))$ $ {T}_{0} $ 井底压力(${P}\left({r}_{0},t\right)$) ${P}_{{\rm w}}$ 外边界压力(${P}\left({\infty },{t}\right)$) $ {P}_{0} $ 外边界温度(${T}\left({\infty },{t}\right)$) $ {T}_{0} $ 2.4 求解结果
我们将分解区气相质量守恒方程（1）与分解前缘气相质量守恒方程（3）作量级分析（具体步骤见附录A）。量级分析的思想是，如果方程是基于无量纲和归一化变量的表现形式，方程不同项的系数能用来度量这些项的重要程度[27]。无量纲化后的方程（1）中第二项系数远大于第一项系数1，可忽略方程（1）中第一项对方程的影响。另外，由于气水相运动黏度相差100倍，在相同压力梯度下，水合物分解后，达到传输平衡，水相饱和度认为是不随时间变化。因此，在这个模型中水合物分解区流体流动简化为拟定常流动，得到分解区压力传导关系：
$$ {P}=\Bigg({P}_{{\rm w}}^{2}+\frac{{P}_{{\rm D}}^{2}-{P}_{{\rm w}}^{2}}{S\left(t\right)}S(t)\Bigg)^{1/2} $$ (7) 根据分解区压力传导方程（7）与分解前缘质量守恒方程（3）得到分解前缘随时间移动规律：
$$ {S}\left({t}\right)=\sqrt{\dfrac{\dfrac{{KK}_{{\rm{rg}}}}{{u}_{{\rm g}}}\dfrac{{P}_{{\rm D}}^{2}-{P}_{{\rm w}}^{2}}{{P}_{{\rm D}}}t}{\phi \Bigg[\dfrac{{\rho }_{{\rm h}}{S}_{{\rm h}}\varepsilon }{{\rho }_{{\rm g}}}-\left({1-S}_{{\rm w}}\right)\Bigg]}} $$ (8) 由式（8）可知分解前缘的移动与流体相渗透率，水合物分解相平衡压力与井底压力之间的差值有关，与时间的平方根呈线性关系。
将分解区与未分解区能量守恒方程作无量纲转换后代入积分得到（具体步骤见附录A）：
$$ \begin{split} &0 {\text{＜}} {{x}} {\text{＜}} {{S}}\left( {{t}} \right):{{T}} = {T_{\rm w}} + ({T_{\rm D}} - {T_{\rm w}})\dfrac{{{\rm{erf}}\left( {\dfrac{x}{{2\sqrt {{a_2}t} }}} \right)}}{{{\rm{erf}}\left( {\dfrac{{S\left( t \right)}}{{2\sqrt {{a_2}t} }}} \right)}}\\ &{{S}}\left( {{t}} \right) {\text{＜}} x {\text{＜}} \infty :{{T}} = {T_0} + ({T_{\rm D}} - {T_0})\dfrac{{{\rm{erfc}}\left( {\dfrac{x}{{2\sqrt {{a_1}t} }}} \right)}}{{{\rm{erfc}}\left( {\dfrac{{S\left( t \right)}}{{2\sqrt {{a_1}t} }}} \right)}} \end{split} $$ (9) 将式（9），（8）代入式（4）并与式（5）联立得到一个用于求解分解前缘相平衡压的超越方程组：
$$ \begin{split} & {\phi }{{\rho }_{{\rm h}}S}_{{\rm h}}{q}_{{\rm h}}\dfrac{{\rm d}S\left(t\right)}{{\rm d}t}-{\mathrm{\lambda }}_{1}\left({T}_{0}-{T}_{D}\right)\dfrac{\mathrm{exp}\left(\dfrac{-{x}^{2}}{4{a}_{1}t}\right)}{\mathrm{erfc}\left(\dfrac{S\left(t\right)}{2\sqrt{{a}_{1}t}}\right)\sqrt{{\text{π}} {a}_{1}t}}+\\ & {\mathrm{\lambda }}_{2}\left({T}_{{\rm D}}-{T}_{{\rm w}}\right)\dfrac{\mathrm{exp}\left(\dfrac{-{x}^{2}}{4{a}_{2}t}\right)}{\mathrm{erf}\left(\dfrac{S\left(t\right)}{2\sqrt{{a}_{2}t}}\right)\sqrt{{\text{π}} {a}_{2}t}}=0 \end{split} \!\!\!\!\!\!\!\!\!$$ (10) 单位横截面积生产井产气速率为：
$$ \begin{split} &{q}_{{\rm g}}=\dfrac{{K}_{{\rm g}}}{{\mu }_{{\rm g}}}\dfrac{\partial {P}(0,t)}{\partial x}=\dfrac{{K}_{{\rm g}}}{{2\mu }_{{\rm g}}}\dfrac{{P}_{{\rm D}}^{2}-{P}_{{\rm w}}^{2}}{{P}_{{\rm w}}S\left(t\right)} \\ & {q}_{{\rm w}}=\dfrac{{K}_{{\rm w}}{\mu }_{{\rm g}}}{{\mu }_{{\rm w}}{K}_{{\rm g}}}{q}_{{\rm g}} \end{split} $$ (11) 由式（11）可以看到，产气速率与分解前移动距离成反比，说明随着天然气水合物的分解，产气量逐渐减少。
当分解前缘移动到S（t）时，单位横截面积生产井总产气体积为：
$$ {{{Q}}_{\rm g}}\left( {\rm{t}} \right) = \int _0^t{q_{\rm g}}{\rm d}t $$ (12) 实际生产井总产气体积为：
$$ {{Q}}_{{\rm g}}\left(\mathrm{t}\right)={h}2\mathrm{{\text{π}} }{r}_{{\rm w}}\int_{0}^{t}{q}_{{\rm g}}{\rm d}t $$ (13) 式中：h为井射孔有效长度，
$ {r}_{{\rm w}} $ 为井孔半径。3. 讨论
3.1 模型验证
Yousif通过水合物砂岩样品降压分解实验，探究了分解前缘移动现象[28]。以Yousif实验为依据，将实验中所设置的相关参数代入模型中，并通过对超越方程（10）以及式（5），（8）联立求解，得到分解前缘移动位置，实验结果与模型计算结果进行对比，如图2所示，实验数据与模型计算结果近似，从而验证了模型结果的可靠性。
![]()
图 2 实验与模型下的分解前缘移动Figure 2. Model match of experimental results for location of hydrate dissociation front3.2 示例分析
以南海神狐海域天然气水合物藏[29-31]实际参数为例，模型计算所需的基本参数如表2所示。不考虑水合物分解过程中冰的生成对模拟结果的影响，井底压力为3 MPa。
表 2 相关物性参数Table 2. The correlated parameters used for calculation参数 数值 水合物密度${\rho }_{{\rm h}}$（kg/m3） 910 岩石密度${\rho }_{{\rm h}}$（kg/m3） 2 650 水合物层厚度h（m） 40 水合物饱和度${S}_{{\rm h}}$ 0.3 水合物层水饱和度${S}_{{\rm w}}$ 0.7 孔隙度$ {\rm{\phi }} $ 0.3 初始地层压力$ {P}_{0} $（MPa） 14.0 初始地层温度$ {T}_{0} $（K） 287.15 井底压力${P}_{{\rm w}}$（MPa） 3.0 单位质量水合物中气体体积分数$ {\rm{\varepsilon }} $ 0.129 气体运动黏滞系数${u}_{{\rm g}}$（Pas） 1.e-5 水相运动黏滞系数${u}_{{\rm w}}$（Pas） 1.e-3 水合物层绝对渗透率K（md） 2.9 井筒地层绝对渗透率K（md） 150 水合物分解吸收热${ {q} }_{ \rm{h} }$（kJ/kg） 430 水合物导热系数${\mathrm{\lambda } }_{{\rm h}}$（W/mK） 2.11 岩石导热系数${\mathrm{\lambda } }_{{\rm s}}$（W/mK） 2.0 水导热系数${\mathrm{\lambda } }_{{\rm w}}$（W/mK） 0.58 水合物比热容${\mathrm{c} }_{{\rm h}}$（KJ/kg K） 2.22 水比热容${\mathrm{c} }_{{\rm w}}$（KJ/kg K） 4.2 岩石比热容${\mathrm{c} }_{{\rm s}}$（KJ/kg K） 1.0 通过表2所示模型参数，进行水合物储层降压开采模型计算，得到200 d内分解前缘随时间变化规律见图3。在天然气水合物降压开采过程中，当地层压力低于水合物相平衡压力时，水合物开始分解，并出现水合物分解前缘。在开采60、120、200 d后，模型计算分解前缘随时间移动距离分别约为33.35、47.17、60.90 m。从图3a中，可以看到随着开采时间的推移，水合物分解前缘移动曲线斜率变小，说明分解前缘移动速度降低。
![]()
图 3 水合物开采特征a. 分解前缘随时间移动规律，b. 分解前缘移动速率随时间变化，c. 产气速率随时间变化，d. 总产气量随时间变化Figure 3. Characteristics of hydrate productiona. The moving of decomposition front with time，b. the velocity of decomposition front varies with time，c. the rate of gas production varies with time， d. the volume of gas production varies with time根据水合物分解前缘移动速率变化（图3b），看出水合物分解前缘移动速率在生产前期达到最大值，随开采时间推移，移动速率变慢最后将趋于平稳。出现这一现象的原因在于，当井底降压开始时，水合物平衡状态被打破，此时压降开始传递，井压与地层压力差作为水合物分解的主要驱动力。随着在地层中压力传递以及水合物分解过程能量消耗，导致地层温度下降，水合物相平衡压力与井压之间的压差减小，分解过程随之变慢；在开采后期，储层能量不足，水合物分解主要依靠储层热量的传导，分解前缘移动速率保持较低的平稳状态，在这种情况下分解主要受储层热物理性质的影响。此时，我们应该考虑储层注热等技术进一步促进水合物分解。
图3c给出了甲烷产气速率（单位高度和单位宽度的产量）随时间变化，我们可以看到，生产井（单位高度和单位宽度）在200 d内,最大产气速率约为250 m3/d，后期产气速率接近20 m3/d。水合物分解气体体积变化与分解前缘移动有关，即井口产气速率受到水合物分解前缘移动的影响，因此，产气速率随时间变化趋势与分解前缘移动速率一致。同时，如图3d所示，在200 d内甲烷总产气量（单位高度和单位宽度的产量）为18 000 m3。
为了进一步分析模型计算结果，将模型结果与2017年南海神狐海域第一次水合物试采情况进行比对。根据中国地质调查局的报告，连续产气60 d后试采结束，累产气量3 ×10 5 m3，平均日产气5 000 m3[32]。另外，天然气水合物储层厚度为40 m，假设井射孔贯穿整个水合物层，井孔半径为0.15 m。我们将60 d试采数据进行曲线拟合，得到图4所示模型计算总产气量与试采总产气量变化。从图4可以看到，在开采10 d内的两者产气量比较吻合，随着水合物开采过程的进一步推进，模型计算总产气值高于与实际开采结果，两者差异增大，但整体相对误差在可接受范围内。模型计算结果和实际试采结果之间存在较大差异的原因主要在于3个方面，首先本文模型中认为井射孔贯穿整个水合物层，水合物降压分解出的甲烷气体能快速有效的向井口方向移动，减小了孔隙中气体积聚所导致的地层压力增大的影响，从而有利于分解甲烷气的产出、后续压降传递和水合物的进一步分解，同时，模型忽略了实际开采中气体的溶解与扩散等作用的消耗；其次本模型没有考虑水合物二次生成、冰的形成等对水合物分解过程以及甲烷气流动的影响，在实际开采中，当储层局部温度低于冰点以下时，会有冰的生成以及水合物的再次形成，固相物质的出现会降低地层孔隙度和地层绝对渗透率，从而影响压降的传递，水合物分解减慢，并阻碍流体向井口方向流动，导致产气量减小，最后，实际的天然气水合物藏是一个三维空间，水合物分解过程是发生在三维空间内的物理化学变化，本文仅是一维简单模型，不可避免与实际情况存在差距。因此，本文模型计算所得天然气水合物降压开采总产气量高于实际开采值，是对总产气量变化的乐观预测，能简单快速地为实际开采提供大方向的参考。
3.3 敏感性分析
为了分析水合物储层相关参数对水合物分解前缘移动距离的影响，本文采用单次单因子敏感性分析方法，除了变量参数外，其他参数均保持在前述的参考数值[33]。
3.3.1 储层初始温度
地层初始温度是影响天然气水合物开采的重要储层参数[34]。探究不同地层温度下单井降压开采天然气水合物分解前缘移动变化规律。为了不改变储层初始压力，将温度变化控制在满足水合物相平衡条件内。
在地层初始压力与生产井压不变的情况下，不同地层温度下（286.15～288.15 K）水合物分解前缘移动变化如图5所示。可以看到，储层温度的改变显著影响水合物分解前缘的移动，当储层温度增大时，分解前缘移动距离增大，储层温度变化1 K时，200 d分解前缘移动距离在参考情况（287.15 K）基础上变化33%。出现这一现象的原因是，天然气水合物分解是一个吸热的过程，储层初始温度越高，储层所能提供水合物分解的能量越多，有利于水合物分解；此外，随着储层温度的升高，相应的水合物相平衡压力增大，水合物分解速度更快。因此，当水合物储层温度很低时，多考虑提高储层温度作为辅助手段来促进水合物分解，提高产气量。
![]()
图 5 不同储层初始温度下水合物分解前缘移动距离Figure 5. The location of hydrate dissociation front at different initial reservoir temperatures3.3.2 储层绝对渗透率
储层渗透率是反映流体运移能力的重要水力学参数[35]。最新的调查数据表明，南海储层类型主要为黏土质粉砂—低渗透粉砂，地层绝对渗透率平均为2～5 md[36]。据此，敏感性分析的地层渗透率范围为2～5 md。
对于不同地层渗透率，图6显示了分解前缘移动随时间变化。正如预期的那样，随着地层渗透率的减小，分解前缘移动速率变慢，移动距离减小。地层渗透率变化对水合物分解有明显的影响，较高的地层渗透率可以提高流体运移速率，有利于地层整体压降传递，为天然气水合物分解提供更大的驱动力，水合物分解前缘移动距离随之增大，产气量随之增加。当储层渗透率提高1 md（以2.9 md情况为基准），分解前缘移动距离较参考情况增大17%。
![]()
图 6 不同储层绝对渗透率下水合物分解前缘移动距离Figure 6. The location of hydrate dissociation front at different absolute permeability3.3.3 储层孔隙度
另一个重要的地层参数是孔隙度，而渗透率又是与孔隙度有关函数，在研究储层孔隙度变化对水合物分解前缘移动影响时，根据Kozeny-Carman方程[37]（孔隙度的三次方与渗透率之间成正比），在改变地层孔隙度的同时，地层渗透率也随之改变。
随着地层孔隙度的增大，分解前缘移动速率减小，移动距离也随之减小（图7）。在孔隙度为0.38时，此时水合物相平衡压力为3.05 MPa，井口和井口前端之间的压差很小，这意味着岩石的热量不足以用于水合物分解，此时只能通过水合物未分解区域的热量流入来获得额外的能量。在这种情况下，分解前缘移动由储层热物理参数决定。这是因为地层孔隙度很大时，含天然气水合物的沉积物单位体积潜热也较大，但由于天然气水合物含量高，比热低，热导率低，使得地层的热导率变小，因此,地层温度下降更快，从而不利于水合物分解。
![]()
图 7 不同储层孔隙度下水合物分解前缘移动距离Figure 7. The location of hydrate dissociation front at different reservoir porosity4. 结论
（1）通过参数量级分析，将天然气水合物分解渗流场简化为拟定常场，积分求解得到分解前缘随时间移动规律，同时将求解温度场变化方程进行无维化转换，推导出温度分布方程，并根据Yousif对含甲烷水合物砂岩样品降压分解实验进一步验证了模型结果的可靠性。
（2）通过南海神狐海域相关参数的实例分析，发现多孔介质中水合物分解前缘移动速率随时间减小，移动距离与时间的开平方呈线性关系；井口气体产量随时间减小，最后趋于稳定。另外，对比了南海水合物试采结果与模型计算的总产气量，发现模型计算值高于实际开采结果，且两者相对误差在25%范围内。
分析误差出现的原因，主要是建立的模型忽略了气体的溶解与扩散等作用，且没有考虑水合物二次生成、冰的形成等对水合物分解过程以及甲烷气流动的影响。因此，本文通过研究水合物分解前缘移动规律对产气速率、产气量进行了较为乐观的预测，为水合物开采潜力评价提供了一种新的简单快速的计算方法，同时对开采动态评估给出大方向的参考。
（3）通过对地层初始温度、绝对渗透率以及孔隙度敏感性分析发现，地层初始温度和渗透率与水合物分解前缘移动距离之间成正相关关系，初始地层温度对水合物分解过程影响显著，另外，地层孔隙度越大，分解前缘移动速率反而降低，移动距离减小，井口与分解前缘压差减小，此时分解前缘移动速率由储层热物理参数决定。
附录A. 模型求解
将分解区气相质量守恒方程（1）中基本物理量表示为特征值与无量纲数的形式：
$$\tag{A-1} \frac{{\phi }{{ P}}}{{{ T}}{T}_{0}}\frac{\partial }{\partial \bar{t}}\left({S}_{{\rm g}}\frac{\bar{p}}{\bar{T}}\right)-\frac{K}{{\mu }_{{\rm g}}}\frac{{\mathrm{{\rm P}}}^{2}}{{T}_{0}{L}^{2}}\frac{\partial }{\partial \bar{x}}\left({K}_{{\rm r}{\rm g}}\frac{\bar{p}}{\bar{T}}\frac{\partial \bar{p}}{\partial \bar{x}}\right)=0 $$ 式中：
$ {T}_{0} $ 、P、T，L分别为温度、压强、时间，长度的特征值，$ \bar{T} $ 、$ \bar{p} $ 、$ \bar{t} $ ，$ \bar{x} $ 分别为温度、压强、时间，长度的无量纲数。对于分解前缘气相质量守恒方程（3）中基本物理量表示为特征值与无量纲数的形式：
$$\tag{A-2} \begin{split} & {\rm{\phi }}\left({\rho }_{{\rm h}}{S}_{{\rm h}}\varepsilon -{\rho }_{{\rm g}}{S}_{{\rm g}}\right)\dfrac{L}{T}\dfrac{{\rm d}\bar{s}}{{\rm d}\bar{t}}=\dfrac{K{K}_{{\rm r}{\rm g}}{\rho }_{{\rm g}}}{{\mu }_{g}}\dfrac{{{ P}}}{L}\dfrac{\partial \bar{p}}{\partial \bar{x}}\\ & \dfrac{{\rho }_{{\rm h}}{S}_{{\rm h}}\varepsilon }{{\rho }_{{\rm g}}}\dfrac{1}{1-{S}_{{\rm w}}}-1\sim\dfrac{K}{{\mu }_{{\rm g}}\phi }\dfrac{{T}{{ P}}}{{L}^{2}} \end{split} $$ 由于式（A-2）左右项相似，左边项远远大于1，右边项也远大于1，则式（A-1）第一项系数远小于方程第二项系数。
在这个模型中分解区的流动就可以看作拟定常流动：
$$\tag{A-3} \frac{K}{{\mu }_{i}}\frac{\partial }{\partial x}\left({K}_{ri}{\rho }_{i}\frac{\partial p}{\partial x}\right)=0,i={\rm w},{\rm g} $$ 代入初边值条件于式（9）中，得到分解区压力传导关系：
$$\tag{A-4} {P}=({P}_{{\rm w}}^{2}+\frac{{P}_{{\rm D}}^{2}-{P}_{{\rm w}}^{2}}{S\left({\rm t}\right)}S({\rm t}){)}^{1/2} $$ 根据分解区压力传导方程（A-4）与分解前缘质量守恒方程（3）得到分解前缘随时间移动规律：
$$\tag{A-5} {S}\left({t}\right)=\sqrt{\dfrac{\dfrac{{KK}_{\rm {rg}}}{{u}_{{\rm g}}}\dfrac{{P}_{{\rm D}}^{2}-{P}_{{\rm w}}^{2}}{{P}_{{\rm D}}}t}{\phi \left[\dfrac{{\rho }_{{\rm h}}{S}_{{\rm h}}\varepsilon }{{\rho }_{{\rm g}}}-\left({1-S}_{{\rm w}}\right)\right]}} $$ 将分解区与未分解区能量守恒方程作无量纲转换：
$$\tag{A-6} \begin{split} &{t}=\mathrm{\tau }\bar{t}\;\; {{x}} ={S}\left({t}\right)\bar{x} \;\;\mathrm{T}={T}_{{\rm D}}+({T}_{0}-{T}_{{\rm D}})\bar{T} \\ & \mathrm{\tau }={S\left(t\right)}^{2}\rho c/\lambda ={S\left(t\right)}^{2}/a \end{split} $$ 将式（A-6）代入式（2）后，我们进一步将无维化方程作不变性伸缩变换：
$$\tag{A-7} {{t}}^{*}=\mathrm{\beta }\bar{t} \;\;{\mathrm{x}}^{*}=\mathrm{\delta }\bar{x} \;\; {T}^{*}=\mathrm{\gamma }\bar{T} $$ 将无维化方程作不变性伸缩变换后得到γ=1，
$ {\delta }^{2}=\beta $ 。据此，得到一个新的隐式关系：$$\tag{A-8} {F}=\left(\bar{x}{\bar{t}}^{-\frac{1}{2}},\bar{T}\right) $$ 根据式（A-8），式（3）无维以后的解为
$ \bar{T}= \theta \left(\eta \right)=\mathrm{\theta }\left(\bar{x}{\bar{t}}^{-\frac{1}{2}}\right) $ ，代入能量守恒方程后积分得到：$$\tag{A-9} \begin{split} &0 < {{x}} < {{S}}\left( {{t}} \right):{{T}} = {T_{\rm w}} + ({T_{\rm D}} - {T_{\rm w}})\dfrac{{{\rm{erf}}\left( {\dfrac{x}{{2\sqrt {{a_2}t} }}} \right)}}{{{\rm{erf}}\left( {\dfrac{{S\left( t \right)}}{{2\sqrt {{a_2}t} }}} \right)}}\\ &{{S}}\left( {{t}} \right) < x < \infty :{{T}} = {T_0} + ({T_{\rm D}} - {T_0})\dfrac{{{\rm{erfc}}\left( {\dfrac{x}{{2\sqrt {{a_1}t} }}} \right)}}{{{\rm{erfc}}\left( {\dfrac{{S\left( t \right)}}{{2\sqrt {{a_1}t} }}} \right)}} \end{split} \!\!\!\!\!\!\!\!\!\!$$ 将式（A-9），（A-5）代入式（4）并与式（5）联立得到一个用于求解分解前缘相平衡压的超越方程组：
$$\tag{A-10} \begin{split} & {\phi }{{\rho }_{{\rm h}}S}_{{\rm h}}{q}_{{\rm h}}\dfrac{{\rm d}S\left(t\right)}{{\rm d}t}-{\mathrm{\lambda }}_{1}\left({T}_{0}-{T}_{{\rm D}}\right)\dfrac{\mathrm{exp}\left(\dfrac{-{x}^{2}}{4{a}_{1}t}\right)}{\mathrm{erfc}\left(\dfrac{S\left(t\right)}{2\sqrt{{a}_{1}t}}\right)\sqrt{{\text{π}} {a}_{1}t}}+\\ & \quad \quad {\mathrm{\lambda }}_{2}\left({T}_{{\rm D}}-{T}_{{\rm w}}\right)\dfrac{\mathrm{exp}\left(\dfrac{-{x}^{2}}{4{a}_{2}t}\right)}{\mathrm{erf}\left(\dfrac{S\left(t\right)}{2\sqrt{{a}_{2}t}}\right)\sqrt{{\text{π}} {a}_{2}t}}=0 \\[-32pt] \end{split}$$ -
[1] 上海市文物保管委员会. 上海市松江县广富林新石器时代遗址试掘[J]. 考古,1962(9):46.[Archaeological Section Shanghai Municipal Commission of Preservation Ancient Mounments. Trial excavation of Guangfulin Relics in Songjiang county, Shanghai[J].Archaeology,1962(9 ):46.]
[2] ZHANG Yu-lan, ZHANG Min-bin, SONG Jian. Development of ancestors'cultivation revealed in phytolith assemblages from Guangfulin relics[J]. Chinese Science Bulletin, 2003, 48(3):287-291.
[3] 宋建. 王油坊类型与广富林遗存[M]//河南省考古研究所.华夏文明的形成与发展.郑州:大象出版社,2003:183-189.[SONG Jian.Wangyoufang type and Guangfulin remains[M]//Archaeology Institute of Henan Province.Form and development of civilization of an ancient name for China.Zhengzhou:Daxiang Press,2003:183 -189.]
[4] 栾丰实. 王油坊类型初探[M]//栾丰实.海岱地区考古研究.济南:山东大学出版社,1997.[LUAN Feng-shi. Preliminary exploration of Wangyoufang type[M]//LUAN Feng-shi.Archaeological research of Haidai region.Jinan:Shandong University Press, 1977.] [5] 南京博物院考古研究所,杨州博物馆,兴化博物馆. 江苏兴化戴家舍南荡遗址[J].文物,1995(4):16-30.[Archaeology Institute of Nanjing Museum, Yangzhou Museum, Xinghua Museum. Nandang relics of Daijiashan of Xinghua County in Jiangsu Province[J]. Cultural Relics, 1995 (4):16-30.]
[6] 广富林遗址考古队. 广富林遗存的发现与思考[N].中国文物报,2000年9月13日.[Archaeological Team of Guangfulin Relics. Development and Ponderation of Guangfulin Remains[N].Cultrual Relics of China, 13 September, 2000.] [7] 张玉兰,宋建,吕炳全. 广富林遗址考古新发现及先人生活环境探析[J]. 同济大学学报, 2002,30(12):1454-1457. [ZHANG Yu-lan, SONG Jian, LÜ Bing-quan. New discovery in archaeology and exploration and living environment of ancestors from Guangfulin Relics in Shanghai[J]. Journal of Tongji University, 2002,30(12):1454-1457.]
[8] 中国植被编辑委员会. 中国植被[M]. 北京:科学出版社,1983.[Editorial Committee of China Vegetation.Vegetation of China[M]. Beijing:Science Press,1983.] [9] Lin J, Zhang S, Qiu J, et al. Quaternary marine transgressions and paleoclimate in the Yangtze River delta region[J]. Quaternary Research,1989,32(32):296-306.
[10] 张玉兰,张敏斌,宋建. 从广富林遗址中的植硅体组合特征看先民农耕发展[J].科学通报, 2003,48(1):96-99. [ZHANG Yu-lan, ZHANG Min-bin, SONG Jian. Development of ancestors'cultivation revealed in phytolith assemblages from Guangfulin relics[J].Chinese Science Bulletin, 2003, 48(3):287-291.]
下载:
计量
- 文章访问数: 2105
- HTML全文浏览量: 206
- PDF下载量: 12
