弯曲轴对称多孔通道中注入气流的演化

《Physical Review Fluids》：Evolution of a gas current injected into a curved axisymmetric porous channel

　　我们研究了气体注入轴对称水饱和多孔通道的过程，这些通道具有高斯型和抛物型剖面，作为穹顶状背斜中地下储气的理想化模型。利用每个通道的细长特性，我们采用复合渐近近似导出了气液界面的演化方程，该近似能适应较大的通道坡度，并具有描述弱弯曲通道中铺展过程的简化小坡度形式。在高气体迁移率极限下，与平面平直通道不同，浮力在不同几何结构中通过不同机制影响动力学。对于稳定注入高斯通道的气体，由于轴对称性导致的气体速度衰减，浮力可以持续影响流动。在抛物型通道中，通道坡度的增加确保了浮力最终在取决于注入速率和流体性质的时间尺度上影响流动。抛物型通道流动的渐近分析揭示了五个时间域，每个时间域具有多个空间区域和不同的铺展速率，反映了注入与浮力之间时空竞争的发展。最初，一层薄气膜沿上边界铺展；随后，通道坡度和气膜的伸长产生静水压力梯度，该梯度增强直至浮力阻止上接触线移动并使气膜变厚。在气膜下方，液体随后排出，直至界面在浮力作用下变平。降阶模型的解析解捕捉了每个时间域中的界面演化和接触线运动，并与完整数值模拟结果进行了验证。这些结果对地下氢气和二氧化碳储存具有重要意义，其中水平界面垂直推进既提高了安全性又提高了储存效率。I. 引言向更清洁、可再生的能源转型需要大规模储能解决方案。天然地下构造因其巨大的容量、广泛的可用性以及储存不同形式能量的多功能性，提供了一种有前景的选择。能量可以以多种形式储存在地下：通过压缩空气储存机械能[1]，通过加热流体储存热能[2]，以及以气态氢形式储存化学能[3]。对于压缩空气和氢气的情况，多余的可再生电力用于压缩空气或通过电解产生氢气，储存的能量随后在需求超过供应时回收[4,5]。拟议的地下储存地点包括枯竭油气藏、含水层和盐穴[5]。这些地层是多孔的，并充满盐水。虽然水平或缓倾斜几何结构最为常见，但弯曲的背斜为储存氢气等浮力气体提供了特别有利的条件[6]。背斜是由沉积层的向上弯曲形成的，其上由称为盖层的不渗透地层限定，这些盖层作为浮力流体的天然圈闭。背斜与邻近含水层连接的位置称为溢出点；迁移超过该位置的气体可能无法回收。因此，详细了解控制该几何结构中流体运动的流体动力学对于防止气体迁移超过溢出点至关重要。注入氢气在此类结构中的流体动力学构成了本研究的主要动机。铺展流[7]已在水平受限多孔通道中得到广泛研究，特别是当流体不可压缩且运动主要由注入或浮力驱动、受粘性应力抵抗时。在此背景下，Huppert和Woods[8]、Pegler等人[9]以及Zheng等人[10]分析了二维平面域中的注入过程，表明流动结构和铺展速率敏感地依赖于流体之间的粘度比。Nordbotten和Celia[11]以及Guo等人[12]证明，轴对称构型中的演化进一步受一个量化浮力与粘性力相对重要性的参数控制。其他人则纳入了额外的物理效应。Zheng和Neufeld[13]研究了毛细管力的影响，表明对于大邦德数或在单分散多孔介质中，可以恢复尖锐界面极限。Farahani等人[14]通过Darcy-Forchheimer方程建模探索了惯性修正，揭示了在较高注入速率下铺展动力学的变化。Castellucci等人[15]表明，气体可压缩性的影响可在长含水层的整个铺展过程中持续存在，从而减缓迁移并降低气体压力。涉及垂直流动的混溶性重力流也由Szulczewski和Juanes[16]进行了研究。这些研究共同提供了对控制此类流动的关键物理机制在广泛注入情景下的详细理解。然而，天然含水层很少是完全平坦的，这些机制与几何变化的相互作用产生了更丰富、更复杂的流动行为。这种效应可能对流体演化以及二氧化碳和氢气等气体的圈闭和储存效率产生重要影响。MacMinn等人[17]研究了在受背景地下水流动影响的倾斜含水层中二氧化碳的注入和随后的迁移，并纳入了毛细管圈闭效应。他们表明，在注入后，最大储存效率（以保留的二氧化碳体积相对于含水层占据体积衡量）发生在温和的顺坡地下水流动条件下，这是由于毛细管圈闭增强所致。该模型后来在MacMinn等人[18]中扩展以考虑溶解圈闭，发现溶解圈闭通过减缓流体演化和大幅提高储存效率而具有显著影响。Gunn和Woods[19]研究了向带有位于注水井上坡或下坡的泄漏断层的倾斜含水层中的连续注入。他们的结果强调了含水层的几何结构如何决定气体流是否充满通道高度，这将影响储存效率。受不渗透盖层下方二氧化碳流演化的启发，若干研究考察了地形对无约束或深层流动铺展的影响。在这种情况下，滞留流体的运动可忽略至一阶，界面演化由非线性扩散方程控制。Pegler等人[20]分析了由幂律剖面表示的平滑地形变化的影响，而Di等人[21]考虑了更现实的、起伏的曲线地形。先前关于气体注入背斜的研究包括Hagemann等人[22]的计算研究和Mortimer等人[23]的理论分析。Hagemann等人[22]证明流动行为强烈依赖于注入速率。在低注入速率下，重力通过使界面变平来稳定驱替。在高注入速率下，粘性力占主导，导致指进不稳定性和快速横向铺展，使流体能够延伸超过溢出点。Mortimer等人[23]研究了浮力主导状态下的流动，重点关注几何结构和注入条件如何影响总体储存容量。他们考虑了向双层背斜的下层注入，其中各层由具有特定毛细管进入压力的低渗透泥岩隔层分开。他们确定了决定储存效率的临界注入速率。低于该速率时，流体达到足够深度，使得泥岩两侧的静水压力差超过毛细管进入压力。流体随后开始充填上层，从而最大化储存。对于高于此阈值的注入速率，流体铺展超过溢出点并流入邻近含水层。在本工作中，我们在这些研究的基础上，考虑向由不渗透弯曲壁限定的多孔通道中的注入，作为穹顶状背斜的简化模型。我们专注于单层系统，并探索广泛的流动状态，以考察浮力、注入速率、粘性力和几何结构如何控制铺展速率和界面演化。多孔介质的非均质性和各向异性虽然在应用中无疑重要，但在此不予考虑。对向高斯型和抛物型通道中稳定注入的考虑揭示了各种铺展动力学。我们利用含水层的细长特性导出了界面位置的单一演化方程，要求气液界面坡度相对于局部含水层方向较小，或浮力足够强以强制水平界面。为了适应可能在大距离上产生大坡度的弱弯曲含水层，我们使用内禀坐标公式化该演化方程，发展了一种精确纳入浮力效应的复合渐近近似。我们用气体注入的模拟补充复合模型简化版本的渐近分析，该分析提供了动力学主导特征的明确预测。II. 模型公式化 A. 控制方程我们考虑气体向具有弯曲边界的轴对称多孔通道中的注入，该通道最初完全被水饱和（图1）。通道具有固定宽度$H$，假设该宽度远小于通道的特征径向长度尺度$L$。其几何结构由中心线曲率指定，边界位于距中心线法向距离$H$处。中心线曲线表示为\begin{equation} z^* = L f(s^*/L), \tag{1} \end{equation}其中$f$是无量纲形状函数，$L$是通道的特征垂向变化，$s^*$表示沿中心线的弧长。星号表示有量纲变量。在恒定方位角$\theta$的平面中，中心线的曲率由下式给出\begin{equation} \kappa^* = \frac{\mathrm{d}\theta}{\mathrm{d}s^*}, \tag{2} \end{equation}其中$\theta$为水平面的倾角。下面，我们考虑两个不同的极限：弱弯曲状态，其中$L \gg H$（符号“$\gg$”表示“远大于”），使得曲率标度为$\kappa^* \sim H/L^2$（本研究的主要焦点）；以及有限曲率状态，其中$L \sim H$，曲率标度为$\kappa^* \sim 1/H$。图1. 本文考虑的两种薄通道几何结构的示意图，由其中心线（虚线曲线）的曲率定义：(a) 抛物型；(b) 高斯型。在(a)中，显示了圆柱参考系的基向量以及用于公式化模型的曲线坐标系$(s, n)$（指向页内的方位角基向量$\mathbf{e}_\theta$被省略）。还标出了中心线相对于水平的局部倾角$\theta$。注入发生在$s^* = 0$处的点源，浅蓝色区域表示气体充填域。界面与通道上、下边界的交点分别发生在$s^* = s_u^*(t^*)$和$s^* = s_l^*(t^*)$处。我们在与通道中心线对齐的局部轴对称曲线坐标系$(s^*, n^*)$中公式化模型[图1(a)]，其中$n^*$是距中心线的法向距离。中心线与水平面的角度记为$\theta$，对于我们考虑的几何结构，满足$\theta(0) = 0$且$\mathrm{d}\theta/\mathrm{d}s^* > 0$。该曲线坐标系是通过从具有坐标$(r^*, \theta, z^*)$的圆柱坐标描述域进行变换而获得的。由于通道具有均匀宽度$H$，法向坐标满足$n^* \in [-H/2, H/2]$。中心线的曲率通过下式与水平倾角相关：\begin{equation} \kappa^* = \frac{\mathrm{d}\theta}{\mathrm{d}s^*}. \tag{3} \end{equation}通道中一点的位置向量，以笛卡尔基表示但由圆柱坐标$(r^*, \theta, z^*)$参数化，为\begin{equation} \mathbf{r}^* = r^* \mathbf{e}_r(\theta) + z^* \mathbf{e}_z. \tag{4} \end{equation}圆柱坐标中的坐标基向量[图1(a)]定义为\begin{equation} \mathbf{e}_r = (\cos\theta, 0, \sin\theta), \quad \mathbf{e}_\theta = (-\sin\theta, 0, \cos\theta), \quad \mathbf{e}_z = (0, 0, 1). \tag{5} \end{equation}位置向量以圆柱基表示但由曲线坐标$(s^*, n^*)$参数化，为\begin{equation} \mathbf{r}^* = \mathbf{r}_c^*(s^*) + n^* \mathbf{n}(s^*), \tag{6} \end{equation}其中$\mathbf{r}_c^*(s^*)$是中心线的位置向量，$\mathbf{n}(s^*)$是法向单位向量。坐标变换的雅可比元素可写为\begin{equation} \frac{\partial (r^*, z^*)}{\partial (s^*, n^*)} = (1 + \kappa^* n^*) \begin{pmatrix} \cos\theta & -\sin\theta \\ \sin\theta & \cos\theta \end{pmatrix}. \tag{7} \end{equation}对法向导数积分得到\begin{align} r^*(s^*, n^*) &= r_c^*(s^*) - n^* \sin\theta, \tag{8a}\\ z^*(s^*, n^*) &= z_c^*(s^*) + n^* \cos\theta, \tag{8b} \end{align}（对于$1 + \kappa^* n^* > 0$），其描述了径向和垂向坐标沿中心线法向方向的变化。曲线坐标系中的切向和法向单位向量定义为\begin{equation} \mathbf{s} = \frac{\partial \mathbf{r}^*/\partial s^*}{|\partial \mathbf{r}^*/\partial s^*|} = \cos\theta \mathbf{e}_r + \sin\theta \mathbf{e}_z, \quad \mathbf{n} = -\sin\theta \mathbf{e}_r + \cos\theta \mathbf{e}_z, \tag{9} \end{equation}分别如图1(a)所示。体积元为$(1 + \kappa^* n^*) \mathrm{d}s^* \mathrm{d}n^* \mathrm{d}\theta$。假设流体不可压缩且不混溶，由位于$n^* = h^*(s^*, t^*)$处的尖锐界面分开。与下边界和上边界的接触点分别记为$s^* = s_l^*(t^*)$和$s^* = s_u^*(t^*)$，使得$h^*(s_l^*, t^*) = -H/2$且$h^*(s_u^*, t^*) = H/2$。界面投影到中心线上的长度为$\ell^* = s_u^* - s_l^*$。我们将$h^*$的定义扩展到整个域，写作\begin{equation} h^*(s^*, t^*) = \begin{cases} H/2, & s^* < s_u^*, \\ h^*, & s_l^* \le s^* \le s_u^*, \\ -H/2, & s^* > s_l^*. \end{cases} \tag{10} \end{equation}我们使用下标$g$和$l$分别区分气相和液相（水）的变量。各相中的流体速度$\mathbf{u}_j^*$（$j = g, l$）满足\begin{equation} \mathbf{u}_j^* = -\frac{k^*}{\mu_j^*} (\nabla^* p_j^* - \rho_j^* \mathbf{g}^*), \quad j = g, l, \tag{11} \end{equation}这是在曲线坐标系中表达的达西定律，其中$\mathbf{g}^* = -g^* \mathbf{e}_z$是重力加速度，$k^*$是介质的固有渗透率。密度$\rho_g^*$和$\rho_l^*$以及粘度$\mu_g^*$和$\mu_l^*$在各相内假设为均匀。气体通过位于原点（取为通道中心线的顶点）的点源注入。气相和液相的连续性方程为\begin{align} \phi^* \frac{\partial S_g}{\partial t^*} + \nabla^* \cdot (S_g \mathbf{u}_g^*) &= Q^* \delta(s^*) \delta(n^*), \tag{12a}\\ \phi^* \frac{\partial S_l}{\partial t^*} + \nabla^* \cdot (S_l \mathbf{u}_l^*) &= 0, \tag{12b} \end{align}其中$\delta$是表示点源的狄拉克δ函数，$Q^*$是源强度。在两种流体的界面处，我们施加运动学条件（13），该条件将界面的速度与流体速度联系起来，适用于 。式（13）中的下标 和 表示偏导数，而下标 和 标记相。我们忽略任何毛细效应，并假设压力在界面上连续（14）。式（10）中 的定义意味着式（13）确保对于 没有穿过下边界的穿透，对于 没有穿过上边界的穿透。要求液体填充区域的下边界和气体填充区域的上边界无穿透，要求速度满足（15）。在注入之前，我们假设初始气泡完全位于下边界之上，因此下接触线尚未形成。初始构型假定处于流体静力平衡，因此气液界面初始为水平。我们推迟指定初始条件的具体函数形式，因为这取决于式（1）中通道形状函数 的特定选择。气体体积由下式给出（16）。

B. 无量纲化 我们用 、 和 重新缩放带星号的变量，其中 。式（7）中的坐标导数变为（17a）（17b）（17c）（17d）。将重新缩放的变量代入达西方程（11），并投影到 和 方向，得到（18a）（18b），其中 和 分别测量气相和液相中相对于浮力的流体静压力贡献，其余无量纲参数为（19）。这里 是粘度比， 类似于达西-瑞利数[24]，衡量浮力与注入引起的粘性力的相对强度；类似的无量纲群出现在轴对称受限流的其他研究中[11,12]。将这些重新缩放的变量代入质量守恒方程（12a）和（12b），得到（20）。在重新缩放的坐标中，界面处的边界条件（13）和（14）变为（21a）（21b），边界处的无穿透变为（22）。达西方程（18）、质量守恒方程（20）以及边界条件（21a）、（21b）和（22）构成了模型的完整无量纲表述。在式（18）–（22）中，下标 、 和 表示偏导数，而 和 标记相。

C. 细通道近似 在本节中，我们考虑当 时模型无量纲表述的简化。目前，我们假设通道具有强曲率，对应于式（1）中的 。达西方程中大量重力项的存在使极限 成为奇异的。因此，动力学可能表现出多种机制：一个早期时间机制，其中通道坡度和与水平的夹角 保持较小，所有项在切向达西方程中平衡；以及在 在通道中不保持较小的几何形状中，一个后期时间机制，其中浮力主导流动并使界面变平。受此结构启发，我们在简化表述中保留重力项，并推导出一个演化方程，该方程在两种机制中都具有正确的极限行为。取 ，我们在法向达西方程（18）中丢弃 和 ，因此压力在通道上呈流体静力分布；在 方向积分并应用界面压力连续性条件（21b），我们恢复压力场（23），其中某个标量场 表示中心线气体压力。然后，式（18a）和（18b）中的切向达西方程变为（24a）（24b）。式（17b）表明，通道的细长意味着，到一阶，径向坐标独立于 ，因此 ，其中 是通道中心线的径向坐标。在质量守恒方程（20）中忽略 量级的项，则对于 ，得到（25）在一阶。在这个极限中，点源成为施加在 处的边界条件。将式（25）在 上对气相积分，在 上对液相积分，应用运动学条件（13）和无穿透条件（15），得到气相和液相深度积分的输运方程（26a）（26b）。将式（26a）和（26b）相加并积分，恢复通量约束（27），该约束平衡通过通道的体积通量与源输送的体积通量。将式（23）、（24a）和（24b）代入式（27）中的积分，得到（28a）（28b）。将式（28a）和（28b）代入式（27）并重新排列，得到一个将压力梯度 与界面位置联系起来的表达式：（29）。将压力梯度的表达式（29）与深度积分速度（28b）和液体的输运方程（26b）结合，得到界面位置 的演化方程。通过引入 可以获得更简单的表述， 测量界面距下通道边界的法向距离[图1(a)]。我们假设气体密度远小于液体密度 ；因此 。用 表示，演化方程以及相应的边界和初始条件可以写为（30a）（30b）（30c）（30d）（30e）。在式（30a）中，与 成比例的项表示由垂直坐标随弧长的变化引起的沿通道中心线的流体静压力梯度；大的 因子反映了当通道坡度 较大时 随 的强烈变化。与 成比例的项捕捉了与法向方向上流体静压力差相关的浮力贡献。平流项考虑了施加的注入。在接触线附近对式（30a）进行局部分析，得到运动学条件（30b）和（30c），确保没有通量穿过移动的接触线。在界面到达通道顶点 处的下边界之前，应用边界条件（30d）。一旦界面离开原点 且 ，则改为应用无通量边界条件（30c）。式（30e）是表示平直界面的初始条件；这里 是一个常数，定义在 处界面上方距通道中心线的高度。由于气体不可压缩，气体总质量等同于其总体积；因此，气体体积守恒确保（31）。我们用 表示通道中气体的初始体积。为了闭合控制方程（30），必须指定通道中心线的几何形状。对于给定的中心线曲率 ，径向坐标 和局部水平倾角 可以作为中心线弧长的函数计算，使用（32）。一旦这些几何量确定，就可以求解式（30a）–（30e）以获得界面从其规定的初始条件的演化。在这项工作中，我们考虑两种通道几何形状：一种弱弯曲的抛物线通道，它在长距离上产生大坡度；以及一种弱弯曲的高斯通道，它在长距离上渐近变平。每种情况对应的几何量在下面第II E节中给出。在指定这些之前，我们推导出式（30a）–（30e）的一个更简单的版本，该版本适用于小中心线曲率和浅通道坡度。

D. 小坡度极限 当通道坡度较小——要么在整个通道中，要么在通道顶点附近——可以获得控制方程（30a）–（30e）的简化形式。我们引入缩放 和 ，在此缩放下列（32）表明中心线的径向坐标近似等于其弧长， 。在这个小坡度极限下，控制方程（30a）–（30e）则（在一阶）简化为（33a）（33b）（33c）（33d）（33e）。通道坡度通过式（33a）中的 规定；与 成比例的项表示浮力效应。与式（30a）一样，与 成比例的项来自浮力的法向分量，而依赖于 的项则考虑浮力的切向分量。平流项表示注入的影响。在极限 下，式（33）简化为Guo等人[12]及其他研究的平坦轴对称通道的控制方程。用于求解式（30a）–（30e）和（33a）–（33e）的数值方法的细节在附录A中提供。

E. 中心线的几何形状

1. 抛物线通道 对于图1(a)中绘制的抛物线通道，式（1）中的形状函数定义为（34），这对应于式（2）中的缩放 。将式（34）代入式（2）得到通道中心线的曲率：（35）。因此，曲率为 ，即使当 。平面曲线 的弧长由下式给出（36）。将式（34）代入式（36）得到抛物线通道的弧长，（37a）其逆函数定义中心线的径向坐标 。中心线角度 可以通过数值积分获得（37b）。我们将式（32）与式（37a）和（37b）结合，求解式（30）以得到通过抛物线通道的流动。尽管曲率在整个域内保持较小，但局部坡度随弧长单调增长，当 时接近 。因此，小坡度近似式（33a）–（33e）仅在 时有效。在这个渐近区域，我们引入小坡度重新缩放 和 ，在此缩放下一阶变量指定的几何量简化为（38）。将 代入式（33a）–（33e）则完成抛物线通道小坡度演化模型的指定。

2. 高斯通道 我们考虑的第二个几何形状是一个高斯通道，其中心线处处弱弯曲，如图1(b)所示。通道几何形状通过式（1）中的形状函数规定，我们取为 。在小坡度重新缩放下，这在一阶简化为 。通过这种选择，中心线的一阶曲率可以显式写为（39）。曲率在原点附近为负，并在 处改变符号，对应于从局部凹几何向渐近平坦通道的过渡。与抛物线通道相比，高斯通道的坡度增长直到 ，之后指数衰减并在远场渐近趋近于零。中心线相对于水平线的角度近似为（40）。将式（40）代入式（33a）–（33e）完成高斯形状通道中扩散控制方程的指定。

F. 模型公式的总结 总结来说，我们使用复合近似推导了一个演化方程[式（30a）]，用于细长弯曲通道中的气液界面，用通道中心线局部的坐标表示。在此近似下，压力在中心线法向上是流体静力的，而重力强迫以奇异的 项进入切向达西方程。这种结构反映了所选的缩放，其中垂直坐标的变化发生在与弧长相同的长度尺度上。取式（30a）–（30e）的小坡度极限，得到式（33a）–（33e）给出的简化系统，这构成了本研究分析的主要模型。该近似的准确性取决于通道几何形状。对于抛物线通道，如第II E 1节所示，近似要求 ，因为大坡度在注入点下游很远处发展。然而，在弱弯曲通道中，式（33a）–（33e）即使在严格渐近区域之外，也能继续捕捉更一般模型式（30a）–（30e）的基本物理机制。这在图2中有所说明，稍后更详细地讨论。因此，我们采用小坡度系统式（33a）–（33e）作为阐明弯曲通道中主导扩散动力学以及解释复合模型式（30a）–（30e）所描述的更广泛行为的简化框架。相比之下，对于浅坡度的高斯通道，该近似在整个域内仍然有效，因为坡度保持较小。

图2。在抛物线通道中，对于 ， ，和 ，（a）小坡度模型式（33a）–（33e）和（b）复合模型式（30a）–（30e）的数值解[式（38）]。上下接触线的演化由实黑色曲线表示。颜色图表示界面高度 ；对于每个时间 ，颜色的垂直变化表示 在界面上的空间变化，从 到 。在(a)中识别出五个渐近区域（I–V）。时间 和 也被展示。为了帮助解释颜色图，图3中展示了代表五个区域中每一个的膜厚度和界面轮廓。

图3。抛物线通道中扩散的五个渐近区域的说明。每一行对应一个区域，每个区域的时间尺度由垂直箭头指示。第1列显示小 极限下 的空间区域示意图，指示渐近区域的长度尺度和 的特征量级。第2栏展示了方程(33a)–(33e)在柱坐标系中当 和 时的解，黑色曲线标出了抛物线形通道的边界。在本文剩余部分，我们专门考虑极限 ，即气体（如氢气）驱替黏度高得多的液体的特征情形，并为抛物线形和高斯形通道几何构建降阶渐近模型。抛物线形通道的渐近分析详见附录B，高斯形通道的渐近分析详见附录C。在第II G节和第II H节中，我们概述了这些方法的关键步骤。在第III节中，我们将完整模型方程(33a)–(33e)的数值模拟与降阶模型的预测结果进行比较，并突出主导的物理平衡关系。

G. 抛物线形通道：小 极限

对于抛物线形通道，当方程(33a)–(33e)和(38)中 且 时，在不同的时空区域会出现多个不同的极限。为了表征这些极限，我们将解划分为五个不同的时间区域，标记为I–V。这些区域在图2(a)中用垂直虚线标出，对应 和 时小斜率方程(33a)–(33e)的一个代表性解；图2(b)显示了使用完整几何变量（方程(35)、(37a)和(37b)）计算的 时完整复合模型[方程(30a)–(30e)]的对应解。这些图展示了上下接触线在 之前的演化，接触线之间的颜色图在每个时刻指示了界面高度 沿通道的空间变化。对于所选参数值，小斜率模型和完整复合模型的解本质上无法区分。

在区域I中，沿上边界形成一层长而薄的气体薄膜，在 接近1处显示为红色区域。在区域II中，上接触线 的推进速度减慢，且 在薄膜内减小，表明气层变厚。到区域III时，该区域已显著增厚，气体和液体填充层的深度相当。在区域IV中， 在大部分域内都很小，反映了沿通道下边界排液的薄液膜的形成。最后，在区域V中，下接触线追上上接触线，之后两者一起推进；界面随后以 中的陡峭前沿传播，对应物理空间中的平界面。

附录B中的降阶模型揭示了每个区域内驱动动力学的物理机制，从而解释了图2中观察到的主要特征，这将在第III节中展示。图3展示了与区域I–V相关的时间尺度及其层级顺序。在区域I–IV内， 的解表现出多个渐近空间区域。我们将域划分为用字母下标标记的子区域，从原点开始按字典序排列；例如， 表示区域I中最靠近原点的区域。图3的第一列展示了空间区域的示意图、其特征长度尺度以及 的量级。第二列展示了方程(33a)–(33e)在 和 时映射到柱坐标系的代表性解，每个面板对应五个时间区域中的不同时刻。选择这些参数值是为了说明在极限 和 下的典型界面形状，同时避免当 时出现的长程扩展，以便关键特征保持可见。因此，每一行展示了 的空间结构（在内禀坐标中）如何映射到柱坐标系中的界面形状。在每个时空区域内，我们推导了 和 的降阶方程和渐近近似（附录B）。表I汇总了所有区域I–V中接触线的近似结果，这些结果将与方程(33a)–(33e)的数值解在第III节中进行比较。

表I. 抛物线形通道中区域I–V内接触线位置的前导阶渐近近似及相关时间尺度。近似的推导见附录B。

时间 I II III IV V

H. 高斯形通道：小 极限

对于浅斜率高斯形通道，我们再次关注小黏度比极限 ，但考虑更广泛的 范围。我们的渐近框架（概述于附录C）确定了当 时的一个特殊极限，其中浮力和注入达到平衡；该极限作为参数空间中的组织区域。动力学分为早期和晚期时间区域。在早期，当 时，界面保持在原点附近并受边界曲率影响。一旦下接触线形成（ ）且 时，解发展出内外结构：一个靠近 的内区域和一个包含 的更广阔的外区域。在大多数情况下，内区域被抑制或缓慢演化，因此在晚期仍能感受到下边界的曲率。相比之下，外区域由快速扩展的薄膜组成，向远场推进，通道在远场渐近变平。

在外区域，由方程(33a)–(33e)给出的小斜率模型允许自相似解，将控制偏微分方程简化为常微分方程（ODE）。在附录C 2 b中，我们证明了极限 和 分别恢复Guo等人[12]分析的注入主导区域以及Lyle等人[25]、Guo等人[12]等人研究的浮力主导扩展区域。所得ODE在不同 下的解在第III节中展示，数值积分使用附录C 2 b中描述的打靶法进行。

III. 结果

在第III A节中，我们展示抛物线形通道的结果，重点放在参数空间 和 的区域。在第III B节中，我们展示小斜率高斯形通道的结果，覆盖参数空间 的区域。对于所有模拟，初始界面高度设为 。

A. 抛物线形通道

我们首先展示小斜率模型（方程(33a)–(33e)）和复合模型（方程(30a)–(30e)）在 和 下的代表性解。图4(a)显示了两种模型在 和 下界面演化的比较。对于 ，在绘制的时间间隔内，两个解在图形上无法区分。相比之下，对于 ，在 时出现小差异，误差随时间增加而增大，随着弧长 的增长。这种行为与小斜率近似有效性要求 一致。这些比较证实，当 足够小时，方程(33a)–(33e)能捕捉复合模型的动力学。特别是，对于小 ，小斜率模型可靠地再现了区域I–V之间的动力学转变；然而，在区域V中，一旦界面延伸到通道斜率不再足够小的程度，误差不可避免地累积，如 所见。

图4. 抛物线形通道中 和 时的扩展数值解。(a) 界面在 时刻的演化。这些轮廓分别对应区域I–V。实线红色曲线显示小斜率方程(33a)–(33e)的解；虚线曲线显示完整复合方程(30a)–(30e)在 （黑色）和 （蓝色）下的对应解。(b) 复合方程(30a)–(30e)在 下上下接触线的演化。虚线红色和虚线蓝色曲线分别显示区域II和III的渐近预测（表I）。

图4(a)中的界面演化遵循图2和图3所示的模式。最初，沿上边界形成一层薄气体膜（区域I），类似于平坦几何中的小 扩展[9,12]；然而，与平坦几何不同，运动并非仅由注入主导。随着薄膜伸长，沿其长度方向建立静水压力梯度，使浮力增加薄膜厚度并减慢 的推进，导致鼻部隆起[区域II；见图4(a)中的 ]。在区域III中，浮力开始驱动界面在其大部分长度上缓慢且近似空间均匀的重新分布，逐渐使气体下方的液层变薄。注入的气体体积通量主要由区域 和 的排液容纳（图3），因此到达上接触线附近区域 的通量可忽略不计。此过程持续到注入体积几乎完全由浮力驱动的排液容纳，此时上接触线实际上静止， 在 处平台化，如图4(b)所示。同时，下接触线在注入和浮力的共同作用下继续推进。 的平台化由区域II近似在 时刻的 和区域III解在 时预测（见表I）。两种近似在图4(b)中在其各自有效范围内与完整数值解吻合良好。一旦界面下方的液层排干，仅留下厚度 的薄膜沿下边界，扩展过渡到区域IV。

在区域IV中，足够的气体通量被输送到区域 （图3），使上接触线恢复运动[见图4(a)中的 ]。在该区域中，一个陡峭的传播前沿在 中（对应物理空间中的水平界面）连接薄膜和上边界，并以由注入体积决定的速率推进（见表I）。图4(b)显示区域III的 近似（见表I）捕捉了从有效静止到重新传播的转变。一旦薄膜中的液体完全排干，下接触线追上传播前沿，此后浮力在整个长度上使界面变平（区域V）。在该区域中，界面高度 近似线性。扩展速率由注入气体体积决定（参见Pegler等人[20]中描述的“水平顶部”区域），而浮力仅在一个与 成比例的短界面长度尺度上起作用，引起 的逐渐倾斜，从而在物理空间中产生有效的变平。因此，两个接触线以相同的前导阶位置 传播，并带有 修正（表I），该修正考虑了浮力引起的倾斜[图4(b)]。

为了评估不同区域内接触线近似的准确性并验证所提出的标度关系，在图5中，我们将方程(33a)–(33e)和方程(38)数值计算的接触线演化与表I中总结的相应渐近预测进行比较，固定黏度比 ，并变化 的值。在图5(a)中，除 的情况外，在早期时间两个接触线的传播遵循区域I解，此时浮力次于注入。接触线随后经历区域II和III；对这些转变的更仔细检查见图5(b)和5(c)。对于所有考虑的 值，在足够晚的时间 ，如图5(a)所示，两个接触线接近区域V特征的 解。在 的情况下，区域I和II之间的区别不太明显。

图5. 抛物线形通道中接触线的演化，通过数值求解方程(33a)–(33e)获得， 和 在递减的 值下。在所有图中，实线和点划线分别对应于 和 的数值解。虚线表示不同渐近区域的近似，每个近似用颜色标记和匹配。在(a)中，区域I的接触线近似（表I）用虚线红色显示，区域IV和V预测的两个接触线的渐近线[方程(B34)]用虚线黑色显示。面板(b)将(a)中 的数据重新绘制在区域 和 的坐标中，虚线粉色曲线表示该区域中 的指数近似[方程(B12)]。面板(c)将(a)中 的数据重新绘制在区域 和 的坐标中， 和 。虚线粉色曲线表示方程(B19)中的近似，而黑色点表示 ，即方程(B37)的解，标记 追上移动前沿的时刻。此后，下接触线根据方程(B43a)与前沿一起传播，用虚线黑色显示。

图5(a)中 的数值解在图5(b)的区域 坐标中重新绘制，在区域IV出现之前，曲线在 时塌缩到渐近解方程(B12)上。图5(a)中 的数值数据在图5(c)的区域 和区域 坐标中重新绘制。在初始传播期间， 的曲线再次塌缩到方程(B19)的预测上，直到下接触线截获前进的前沿，标志着区域V的开始。由式(B37)预测的区域IV与V之间的过渡时间由黑点标出。在两种情况下，数值解与渐近预测式(B12)和(B19)之间的高度一致性表明，这些近似准确捕捉了图5(a)中展示的从区域I动力学到区域V扩展特征的转变。因此，当 时，区域I–V的接触线渐近解共同捕捉了整个演化过程中的主导行为。当 时，数值解在图5(a)中与区域I和II的解明显偏离，并在图5(c)中与区域III和IV的解明显偏离。这表明渐近近似在 时未能完全捕捉区域I–IV中的动力学。这种差异产生的原因是，与 成比例的局部浮力项开始影响空间区域 。图6比较了当 和 时，方程(33a)–(33e)的数值解与附录B中推导的跨越区域I–V的渐近近似。在图6(a)中，对应于区域I，早期动力学由近似式(B5)很好地描述，该式准确捕捉了界面轮廓。随着 增大至 ，浮力效应变得越来越显著，导致薄膜逐渐增厚，上部接触线的传播速度降低。因此，当解接近区域II的过渡时，近似式(B5)高估了界面的横向范围，并与数值解产生偏离。通过区域II的演化如图6(b)所示。以 处的数值轮廓初始化双曲解式(B9)，该近似捕捉了界面的主导行为。然而，随着时间推移，在前缘附近出现差异，标志着区域II描述失效和区域III的开始。图6. 当 和 时，方程(33a)–(33e)的数值解。面板(a)–(e)显示了跨越区域I–V的代表性时间；彩色曲线表示数值解，黑色曲线表示附录B中相应的渐近预测。(a) 区域I， ；解析近似式(B5)。(b) 区域II， ；双曲近似式(B9)和(B11)。(c) 区域III， ；虚线曲线显示预测的指数衰减式(B21)。(d) 区域IV， ；结合式(B27)、(B30a)、(B30b)和(B36)的复合近似（彩色曲线）与模拟结果（黑色曲线）几乎无法区分。(e) 区域V， ；近似式(B42)。(f) 投影到中心线上的界面长度演化， 。红色和蓝色虚线分别显示由 和 的相应表达式导出的区域II和III近似（见表I）。图6(c)展示了区域III，其中界面从邻近上边界的薄气膜转变为沿下边界的薄液膜；在区域 （图3）中，足够远离接触线处，界面的主体以式(B21)预测的指数速率下降。区域IV的复合近似，由区域 中的解构建[式(B27)、(B30a)、(B30b)和(B36)]，如图6(d)所示，一旦解进入其有效范围 时，该近似在整个界面上提供了良好的一致性。图6(e)显示了区域V中的演化，此时界面变得有效平坦，近似式(B42)准确捕捉了随后的界面传播。最后，图6(f)显示了 的演化，以及由区域II和III中 和 的解获得的近似（见表I），两者均与数值解显示出良好的一致性。B. 高斯通道在极限 下，高斯形状通道中的扩展控制方程(33a)–(33e)和(40)呈现出内外结构（见附录C）。在此框架下，外部区域自相似演化，并由单一参数 控制，该参数衡量浮力相对于注入的相对强度。如附录C所示，对于 ，将外部区域的相似变量 和 引入控制方程(33a)–(33e)和(40)，连同全局体积约束式(31)，得到边值问题(41a) (41b) (41c) (41d) 这些方程描述了在注入和浮力平衡下外部区域的自相似演化。图7(a)显示了气膜厚度 随 增大的方程(41a)–(41d)的解，说明了薄膜如何随着浮力影响的增强而变薄。图7(b)显示了相应的特征值 ，它通过 确定上部接触线位置，作为 的函数。 随 的单调增加表明更强的浮力加速了上部接触线的扩展。对于小 ，扩展由注入主导， 。在极限 下，特征值渐近趋于 。图7. (a) 描述高斯通道中外部区域动力学的方程(41a)–(41d)的解，其中 且相似变量 ，显示为 递增的值。解是使用附录C中概述的打靶法获得的。(b) 相似解的相应特征值 作为 的函数。蓝色和红色虚线分别表示 和 极限下的渐近近似。下面，我们将完整系统方程(33a)–(33e)和(40)的数值解与使用附录C中概述的打靶法获得的方程(41a)–(41d)的解进行比较。我们还考虑了次极限 下的解析解和极限 下的降阶无参数模型（见附录C）。在极限 下，注入主导浮力，方程(41)允许解析解式(C12)，该解此前由Guo等人[12]报告。所得解的內外结构与抛物线通道中的区域I相似。在此区域中，接触线按以下方式演化(42a) (42b) 图8展示了 时方程(33a)–(33e)和(40)的代表性解，这意味着浮力效应较弱。界面的初始演化如图8(a)所示，以圆柱坐标绘制。如抛物线通道区域I的结果所预期，沿上边界快速形成薄气膜。接触线在短暂瞬态后迅速接近渐近预测式(42a)和(42b)[图8(b)]。演化的后期阶段[图8(c)]显示薄膜伸长。将轮廓重新绘制在式(C6)的相似变量中，表明其坍缩到解析解式(C12)上[图8(d)]。图8. 高斯通道中 且 时方程(33a)–(33e)与(40)的数值解。(a) 界面以圆柱坐标在时间 处显示，时间递增由浅至深的颜色梯度表示。通道边界以黑色显示。(b) 上部（黑色实线）和下部（黑色点划线）接触线的演化。为比较，显示了 的后期渐近预测式(42a)（品红色虚线）和 的后期近似式(42b)（绿色虚线）。(c) 时间 处的界面轮廓 。(d) (c)中的轮廓重新绘制在式(C6a)的自相似变量 和 中。红色虚线曲线显示解析相似解式(C12)。将 增加到1（图9）展示了浮力对流动的影响。在下部接触线形成之前，浮力使界面变平，与渐近解式(C2)一致。在此早期阶段，上部接触线 以注入体积设定的速率推进，由近似式(C3)准确捕捉[图9(b)]。一旦 在原点形成，下部接触线迅速移动到附近位置并被钉住，同时薄气膜沿上边界扩展[图9(a)和9(b)]。图9(c)显示了后期内外结构出现时的 。在邻近 的内部区域，界面保持局部水平，接触线保持钉住，如式(C4)所预测。在外部区域，注入和浮力之间的平衡驱动薄膜更快扩展[比较图9(b)与图8(b)]。方程(41a)–(41d)的相似解，通过附录C中的打靶法数值积分式(C10a)–(C10c)获得，预测上部接触线传播为 (43)，与数值解密切一致[图9(b)]。图9(d)显示了重新绘制在相似坐标中的界面轮廓，除最早时间（ ）外，所有时间均坍缩到外部相似解上。内部解式(C4)也显示在其有效范围内表现出良好的一致性。然而，图9(d)显示外部解式(41c)的对数内部极限（回忆Lyle等人[25]报告的边界条件）与内部解式(C4)的常数外部极限不兼容，暗示存在一个被动的中间空间区域，我们不再进一步研究。图9. 高斯通道中 且 时方程(33a)–(33e)与(40)的数值解。(a) 界面以圆柱坐标在时间 处显示，递增时间由浅至深的颜色梯度指示。通道边界以黑色显示。(b) 上部（黑色实线）和下部（黑色点划线）接触线的演化。为比较，显示了 的渐近预测式(C3)（青色虚线）和后期近似式(43)（品红色虚线）。(c) 时间 处的界面轮廓 。(d) (c)中的轮廓重新绘制在方程(41a)–(41d)的外部自相似变量 和 中。红色虚线曲线显示方程(41a)–(41d)的外部相似解，通过附录C中概述的打靶法获得，而绿色虚线曲线显示内部解式(C4)。进一步增加浮力影响，取 ，改变了外部区域动力学（图10）。最初，界面演化与 的情况保持相似，上部接触线 的早期运动再次由式(C3)捕捉[图10(b)]。在后期，邻近 的内部区域与 的情况相似，尽管下部接触线被钉住在更靠近原点的位置[图10(b)和10(c)]。在此参数区域，注入在外部区域次于浮力，薄膜厚度和轴向长度尺度变得与粘度比无关，具有 和 ，如图7(b)中 的大 渐近所示。在这些标度下，式(41a)简化为多孔介质方程(C13a)，该方程描述了无约束多孔层中浮力驱动的扩展[25]。虽然该简化模型的一个版本此前由Guo等人[12]获得，但据我们所知，内外结构的细节尚未被研究。使用附录C中描述的方案对式(C13a)进行数值积分得到 (44)，与完整数值解紧密匹配[图10(b)]。相应的相似解也显示除最早时间（ ）外，在后期区域建立之前，具有极好的一致性[图10(d)]。图10. 高斯通道中 且 时方程(33a)–(33e)与(40)的数值解。(a) 界面以圆柱坐标在时间 处绘制。(b) 上部（黑色实线）和下部（黑色点划线）接触线的演化。为比较，显示了 的早期渐近预测式(C3)（青色虚线）和后期近似式(44)（品红色虚线）。(c) 时间 处的界面轮廓 。(d) (c)中的轮廓重新绘制在式(C13a)–(C13c)的自相似变量 和 中。红色虚线曲线显示方程(C13a)–(C13c)的相似解，使用打靶法计算。IV. 讨论我们已经展示了含水层形状在确定注入气体电流扩展速率方面的核心作用。我们考虑的两个例子（抛物形和Gaussian通道）说明了不同的结果。A. 抛物形通道抛物形通道被视为弱弯曲的，沿整个通道具有无量纲曲率 。在这种几何结构中，局部斜率缓慢但单调地随弧长增加。虽然小斜率方程(33a)–(33e)没有捕捉到后期发展的大斜率，但与完整复合模型的比较（图2和4）表明，降阶模型保留了所有基本物理特性，并在大斜率产生之前准确再现了完整系统的动力学。即使在后期，当斜率变得显著时，界面继续根据预期的物理机制演化，尽管由于几何项如 和 的近似，界面位置和形状可能出现差异[见图4(a)中的 。在建立了小斜率模型在弱弯曲抛物形通道中扩展的有效性之后，我们采用该公式作为研究动力学的主要框架。在抛物形通道中，我们在极限 和 下识别了五个时间区域，对应于高度可动的气相和相对于注入诱导粘性力较弱的浮力。各区域中接触线和相关时间尺度的渐近近似总结在表II中，以有量纲变量表示，并可进行如下物理解释。状态I的时间尺度由注入体积通量沿特征弧长 扩展所需的时间决定，在此期间界面和接触线主要在注入作用下演化。在该状态下，接触线以正比于 的速率推进，其前置系数由注入通量和流体黏度决定。上接触线比下接触线传播得更快，这是由于上边界附近形成了薄气膜，其中流动性更强的气体平衡了黏性更强的液体的通量。然而，轴对称几何结构会减弱速度，而通道坡度增大则增强切向浮力分量，使浮力变得重要。表II汇总了在抛物型通道中扩散所识别的动力学状态，以有量纲参数和变量表示，适用于 和 。第2列列出了定义各状态的特征时间（其中 表示量级相当），第3列和第4列分别给出了下接触线和上接触线的相应近似。系数 和 定义于式（B41）。时间 I II III IV V 状态II对应于气相中浮力效应的起始。浮力在气体中引起量级为 的速度，从而给出浮力开始影响气体流动的时间尺度 。在此阶段，浮力导致气膜增厚，上接触线呈指数减速，最终趋于平台，其有量纲位置由表II中状态II的 近似系数给出。在状态III中，浮力效应在液相中变得显著。浮力在液体中引起量级为 的速度，从而给出浮力影响液体的时间尺度 。物理上，液体响应比气体慢，因为 。超过此时间后，浮力驱动液体从气体充填区域下方排出，导致下接触线在注入和浮力共同作用下呈指数加速。上接触线在注入作用下继续缓慢推进，其运动受可用体积通量限制，因为注入流体的大部分通过下方液体的排出得以容纳。这种通量再分配反映在状态III中 近似的指数前置系数中（见表II）。状态IV的时间尺度由状态III的时间尺度乘以因子 得出。该因子反映了液体从气体充填区域下方几乎完全排出所需的额外时间，仅留下沿下边界的薄残余层，厚度为 。状态V对应于远大于状态IV的时间，此时液体已从气体充填区域下方完全排出，下接触线已追上上接触线，界面近似水平。B. 高斯通道考虑的第二种几何结构是高斯通道，其曲率和坡度均保持较小，满足 和 。通道坡度随弧长呈指数衰减，因此在远场中几何结构变得有效平坦。因此，在后期，边界几何仅影响 附近的内部区域，而外部动力学简化为平坦轴对称多孔通道的动力学。在附录C中，我们在小 极限下进行了渐近分析，识别出一个新的特殊极限，其中 ，由式（41a）–（41d）给出，对应于浮力与注入之间的平衡。该极限同时恢复了Guo等人[12]和Nordbotten与Celia[11]先前报道的注入主导行为[式（C12）]和浮力主导行为[式（C13a）–（C13c）]。这两个区域之间的边界 在 参数空间内平行于Guo等人[12]图7中给出的两个数值驱动的有限参数交叉边界。这些边界限定了一个有限区域，在该区域内注入和浮力同时影响流动。在极限 下，即使通道是平坦的，浮力也能影响轴对称几何结构外部区域的动力学，这与平面通道形成对比，Castellucci等人[15]发现平面通道中浮力效应始终弱于注入。这种差异是由轴对称通道中注入速度的几何衰减引起的：随着半径增大，注入通量分布在不断增大的横截面积上（或者在缓坡近似下，随 增大的弧长上），导致径向速度衰减。这一机制可以通过尺度论证来理解。浮力通过诱导量级为 的径向压力梯度来平滑界面坡度的变化。根据高斯通道的外部尺度[式（41a）–（41d）]，我们发现 和 ，因此界面梯度的有量纲尺度满足 ，从而给出浮力诱导的压力梯度量级为 。源项产生量级为 的水平速度。考虑外部区域在径向长度尺度 上延伸，该速度尺度变为 。由于液体比气体黏度大得多，它决定了维持该速度所需的压力梯度，量级为 。将该压力梯度与浮力诱导的压力梯度平衡，得到浮力影响外部区域的条件，即 。C. 实际意义表III给出了模型应用于 和 注入典型背斜构造时有量纲和无量纲参数的数量级估计。正如预期， -盐水和 -盐水的黏度比均处于小 状态。 的合理取值范围跨越多个数量级，反映了不同储层中渗透率、通道高度和流体性质的变化，以及可能采用的不同注入速率。表III列出了有量纲和无量纲模型参数。氢气的流体性质取自Castellucci等人[15]， 取自Bickle等人[26]。代表性背斜的通道宽度和径向范围遵循Mortimer等人[23]和Harati等人[27]。由于现场数据有限[27,28]，氢气注入速率反映了数值模拟中常用的值，而 速率基于当前场地注入速率1–2 的估计[26]。渗透率值对应适合储存的典型含水层范围[15]。参数 氢气 气体密度 () 420–610 气体黏度 () 3–5 密度差 () 390–680 注入速率 () 盐水密度 () 盐水黏度 () 通道宽度 () 10–50 径向长度尺度 () 渗透率 () 无量纲参数： 0.3–0.5 对 的估计表明，在高斯通道中，流动处于注入主导状态，后期扩散由式（C12）描述。相比之下，与 储存相关的更高黏度和通常更低的注入速率导致更宽的参数范围，浮力效应可能在某些储存场地的动力学中发挥作用。浮力主导状态导致更薄的外部气膜，因此与注入主导状态相比，储存效率较低。以往研究表明，对于足够缓慢的背斜注入，浮力在早期使界面趋于平坦[23]。相比之下，在较高注入速率下，长气指沿上边界发展并快速传播[22]，这与抛物型通道状态I和高斯通道中状态 观察到的行为一致。对于抛物型通道，我们的结果表明，即使在较高注入速率下，只要背斜在足够大的长度尺度上延伸，浮力最终也会恢复平坦界面，尽管时间较晚。此外，在曲率更大（由 量化）和渗透率更高 的通道中，平坦化的起始更早（表II）。对于地下储氢而言，平坦界面是有利的，因为它最大化了通道体积的利用率。此外，弯曲背斜具有溢出点，超过该点气体将损失至邻近含水层[3]。值得注意的是，我们发现在渐近状态 和 下，上接触线在 处趋于平台，并在 时保持有效钉扎。如果平台位置位于溢出点上游，则浮力可能通过阻止 的前进而抑制溢出。随后的注入主要通过下方液体层的排出得以容纳，而非上接触线的进一步传播，从而也最大化了可用通道体积的利用率。该机制提出了一种实用策略，即在保持对当前演化控制的同时提高储存能力，即使在相对较高的注入速率下也是如此。使用表III中列出的参数范围，我们可以估计抛物型通道中状态I–V相关时间尺度的数量级。考虑以速率 将氢气注入近似抛物型背斜，其渗透率为 ，特征宽度为 ，垂向坐标在径向长度尺度 上变化 。利用氢气与盐水之间的密度差以及表III中列出的氢气黏度，我们得到 。状态II的时间尺度约为一天，而状态III和IV预计分别在大约115天和1.5年后出现。在更长时间后，状态V出现，界面将被浮力有效平坦化。由于氢气-盐水体系的小黏度比（ ），当 足够小时，状态IV–V在实际运行时间尺度上不太可能实现。然而，在曲率更大的背斜中，该限制可能放宽，此时动力学由式（30a）–（30e）控制，或者对于 注入，由于其黏度更大。此外，在实际应用中，流动是否达到这些晚期状态还将取决于注入给定储存场地的气体总体积。这可以通过比较注入给定目标体积 所需的时间与不同状态的时间尺度来评估。尽管本节主要关注 和 储存，但该建模框架也可能为压缩空气等其他浮力气体的储存提供见解[29]。然而，气体可压缩性和热力学效应预计将在气泡演化和整体储能效率中发挥更突出的作用[30]。可压缩性效应可以采用类似于Castellucci等人[15]开发的方法纳入当前模型。D. 结语总而言之，我们扩展了现有的受限多孔层中黏性重力流理论模型，以考虑弯曲通道几何结构。通过在通道中心线局部坐标中建立控制方程，并利用通道的细长特性，我们构建了一个复合近似，得到气液界面的演化方程（30a），精确保留了浮力效应。该框架描述了注入气体流经具有有限曲率的细长通道的演化，并能适应大通道坡度。模型以一般几何函数表示，因此在完整复合方程（30a）–（30e）中指定曲率 、倾角 和径向坐标 即可从给定初始条件模拟广泛几何结构下的界面演化。利用复合模型，我们导出了适用于弱弯曲小坡度通道的简化演化方程（33a）。受穹顶状背斜中氢气与 储存的启发，我们使用小坡度模型研究了抛物型和高斯型背斜中的气体注入。为了隔离几何效应，在推导式（30a）–（30e）和（33a）–（33e）给出的模型时做了若干简化假设。特别地，我们假设界面保持轴对称。然而，众所周知，黏性流体被黏度较低的气体驱替（即 ）容易发生黏性指进不稳定性[31]。Hinton和Jyoti[32]最近的工作考察了平坦轴对称通道中气液界面的稳定性，其状态为 和 ，采用允许界面方位角变化 的三维公式。他们表明，指进可能在早期出现，但一旦浮力将流体分层，指进不稳定性即被抑制。一旦分层发生，注入气体沿上边界形成径向广泛扩展的薄膜，如抛物型通道状态I和高斯通道中 所观察到的，总压降分布在更大的径向范围上。局部压力梯度的减小降低了其对微小界面扰动的敏感性，从而抑制了导致黏性指进的的正反馈机制。在该自相似区域内，作者通过线性稳定性分析证实，界面对轴对称和非轴对称扰动都是稳定的。在某些条件下，其他物理机制可能变得重要。对于足够小的孔隙尺寸或大的界面张力，毛细效应可能产生毛细 fringe 和跨越界面的压力跃变，从而使尖锐界面假设失效[13]。此类效应还可能促进毛细捕集。如果溶解时间尺度远小于储层特征时间尺度，则溶解度可能通过降低当前迁移速率显著影响演化过程。这种效应在小- 区域尤为显著，因为流体间增强的界面面积促进了传质增加，正如MacMinn等人[18]所证明的那样。同样，如果沿流向的压力降足够大，则气体可压缩性可能减慢扩散速率[15]。此外，流体性质依赖于压力和温度，因此在不同含水层之间可能有所差异。温度的空间变化可能引入额外的动力学效应，包括热驱动对流，或者在存在浓度梯度时，双扩散对流[33]。在本研究中，我们关注弱弯曲通道中的重力流，由小斜率近似方程(33a)–(33e)描述。然而，完整的复合模型对于曲率更强的通道仍然有效，为研究高曲率细长几何提供了框架。在更一般的情况下，方程(30a)和(30b)中与 成正比的项（在小斜率近似中被忽略）可能变得显著并影响界面演化。该项修正了浮力的切向分量，可能影响动力学行为。一个自然的扩展是考虑包含溢出点的几何构型，超过该点通道斜率变为正。在这种情况下，注入的气体可能分离，需要模拟多个不连续气泡。我们将此类扩展留给未来的工作。