CN101903805A - 沉积盆地中的建模 - Google Patents

Abstract

本发明的实施例用于产生从压实和流体流动方面描述盆地的盆地模型。可以同时求解定义压实和流体流动的方程。本发明的实施例使用定义一组未知量的多个方程，所述方程在盆地上是相容的(有解的)。所述方程可以定义总压力、流体静压、厚度和有效应力。

Description

沉积盆地中的建模

相关申请的交叉引用

本申请要求于2007年12月21日递交的、名称为“MODELING INSEDIMENTARY BASINS”、代理人卷号为2207EM385的美国临时专利申请61/008,801的权益，该临时专利申请通过引用整体并入本文。

技术领域

本申请一般涉及计算机建模，且更具体地涉及对沉积盆地中的压力建模。

背景技术

在地质勘探中，获得与存在于地表以下的各种岩层和结构有关的信息是令人期望的。此类信息可包括地层、密度、孔隙度、组成等信息。之后，此类信息用于建模地表下的盆地以预测碳氢化合物储藏并且在碳氢化合物的开采方面提供帮助。

盆地分析是对作为地球动力学实体的沉积盆地的综合研究。之所以研究沉积盆地是因为这种盆地包含随着时间的推移发生在地表之上和之下的沉积过程的记录。在盆地的几何体系中，盆地包含地壳构造演变和地层史以及岩石圈如何变形的迹象。因此，盆地是地质信息的主要储存库。此外，过去和现在的沉积盆地是几乎所有世界上商用碳氢化合物矿床的来源。

盆地仿真或模拟对沉积盆地的构成和演变建模。仿真致力于各种物理和化学现象，这些现象控制在移动正在下沉的盆地的构架的过程中碳氢化合物矿床的形成，例如热传递、密实度、水流量、碳氢化合物产生以及流体的多相移动。盆地建模提供对流体流动和孔压模式的重要了解。注意压力评估对勘探估计和计划都很重要，因为压力能够接近一些压实不足区(under-compacted area)中的静岩。

在盆地的典型历史中，沉积物在层顶部的沉积随着时间的推移逐渐累积形成另一层。随着越来越多的层添加到顶部表面，地下层被顶表面层的重量压实。地下层的孔隙度也随压实而改变。因此，随着时间的推移，孔隙度一直在变化。在盆地形成的过程中，可能会在一沉积层的顶部上形成有机物质层。随着时间的推移，有机层被其他沉积层覆盖。该有机物质层被称为源岩。源岩暴露在热量和压力下并且有机物质被转换成碳氢化合物矿床。随后的压力导致碳氢化合物物质从源岩排出并且迁移到圈闭(entrapment)位置。因此，对于盆地建模来说，至关重要的是了解碳氢化合物在源岩中形成所处的条件(例如，温度和压力)以及碳氢化合物在其迁移期间被暴露或已经暴露的条件。准确的建模将实现对盆地的更为成功的勘探。

主要条件之一是压力，其可以利用达西定律来定义，达西定律的内容是，液体会从较高的压力区移至较低的压力区并且移动的速率与压降成正比。对于多孔介质中与实验室压实规律和应力应变行为有关联的单相流体流动，非平衡压实以及由此产生的水流可利用达西定律表示。一个例子可以在作者为P.A.Allen和J.R.Allen、1990年在马萨诸塞州剑桥的Blackwell Scientific Publications上的“Basin Analysis：Principles and Applications”中找到。这种耦合过程的数字建模是复杂的，并且历史上已在三个领域实施：以计算应力应变行为为主目标的地质力学建模、在多孔介质中的流体流动建模和断裂力学。请注意，对于涉及这些过程中两个或三个过程的建模，建模总是假设这些过程是分开的。换句话说，每个过程独立于其他过程而建模。因此，这种方案在这些过程存在强联系的情况下是不可接受的(例如，在高沉积速率的情况下)，此时，由于应力改变而造成的孔隙度和渗透率的快速变化会导致压实不足、关于流体静力分布的高过压的形成以及可能情况下的固体介质断裂。这种分开方案的例子中可以在作者为I.L′Heureux和A.D.Flowler的″A simple model of flow patterns inoverpressured sedimentary basins with heat transport and fracturing″，Journal of Geophysical Research，Vol.105，No.BlO，23741-23752页，2000中找到。

发明内容

本说明书针对通过评估盆地中发生的现象准确建模地质盆地中状态的系统和方法的实施例。这种建模可包括描述随地质时间而演变的沉积盆地中的压实过程和流体流动。

在建模压实过程和流体流动的同时，考虑包括多孔固相的沉积系统，多孔固相的空隙体积充满称为孔隙流体的液体。由于重力作用和固相与液相之间的密度差异，固相通过减少其孔隙度在其自身的重量(以及其他层的重量)下压实，从而导致孔隙流体从固相基质排出。

本发明的实施例使用连续力学方法表达质量和动量守恒的方程。本发明的实施例假设一维的垂直压实(压紧)以简化压实现象。这允许本发明的实施例同时求解流体流动和压实的方程。

本发明的实施例利用达西定律支配的一维垂直压实和三维孔隙流体运动导出一个非线性方程系统。一个方程是利用关于流体静力(hydrostatic)负荷的超压表示的扩散方程。另一个方程涉及固岩厚度及其孔隙度。另一个方程利用力平衡定义有效应力。又一方程是将总的垂直应力和孔隙压力与孔隙度关联的本构律。该方程假定岩石基质的弹塑性行为，换句话说，岩石的压实状态是不可逆转的，并且呈现出滞后性。

本发明的实施例使用将流体密度与压力、多孔岩石的渗透率及其孔隙度关联的本构律。尽管本发明的实施例可以使用任何现有的关系，但是流体密度ρ_a对压力p的依赖关系被假定为其中

是大气压ρ_atm下的流体密度，且渗透率对孔隙度K的以下普通依赖关系假定为

其中K₀，n，m是某些常数。

本发明的实施例运行以更为有效、更少的时间并且利用较少的计算资源的方式产生盆地模型。本发明的实施例允许同时解决压实和流体流动而不是利用重复迭代。本发明的实施例甚至在地质盆地经历快速变化时(例如沉积物以高速率沉积时)也产生准确的结果。

在一个一般方面，例如在计算机上建模物理区域的方法包括：接收定义物理区域的至少一个物理特性的数据；选择第一现象和第二现象，其中第一和第二现象在物理区域上联系以进行建模；定义描述第一和第二现象的一组方程，其中这些方程在物理区域上是有解的(相容的)；通过对第一现象、第二现象和该组方程中的至少一种施加至少一种假设来简化该组方程；和求解该组方程以同时利用数据描述这两种现象。

此方面的实施可包括以下特征的一个或更多个。例如，物理区域可以是地表下的地质盆地并且这两种现象可以是流体流动和盆地中的物质压实，其中流体位于盆地内。流体可以是石油、天然气、水、液体、气体和具有放射性同位素的液体中的至少一种。所述物质可以是沉积物。所述至少一种假设可以包括以下中的至少一种：沉积物累积率是已知的；压实仅发生在垂直方向；和/或压实是相对不可逆转的。该方法可包括在物理区域的模型上提供网格，其中网格包括多个单元。求解可以针对网格的每个单元进行。在建模期间，可以在垂直方向发展网格的每个单元以对物质随时间的累积建模。在建模过程中，至少一个单元可能变得掩埋在模型中，因为在该一个单元上方发展其他单元。每个单元可以是平行六面体单元。定义单元的水平平面的x方向和y方向与地层的时间线(或时代线)对准。

流体可以是可压缩流体，并且方程组可包括第一方程、第二方程、第三方程和第四方程，其中第一方程定义每个单元的过压，第二方程定义每个单元的单元厚度，第三方程定义每个单元的物质负荷，第四方程定义每个单元的流体静压。流体可以是不可压缩的流体，并且该组方程可包括第一方程、第二方程和第三方程，第一方程定义每个单元的过压，第二方程定义每个单元的单元厚度，第三方程定义每个单元的物质负荷。向单元施加至少一种转换；其中这种转换是沉积、下举(downlift)、抬升和侵蚀中的一种。可以向邻接区域边缘的单元施加至少一个边界条件。物理区域可以是地表下的地质盆地，并且该模型涉及地表下的石油，并且求解有助于从盆地开采石油。数据可能根据来自传感器的信息而获得，所述传感器测量物理区域的至少一种物理特性。

该方法可包括基于该组被求解的方程产生地表下的地质盆地的盆地模型。可以基于盆地模型预测物理区域内碳氢化合物的位置。可以安排生产基础设施以在物理区域内基于预测的碳氢化合物的位置开采碳氢化合物。物理区域生产碳氢化合物的潜力可以基于盆地模型来安排。

在另一一般方面，公开了具有其上记载了计算机程序逻辑的计算机可读介质的计算机程序产品，用于在计算机上建模地表下的地质盆地，其包括：定义描述流体流动和沉积物压实的一组方程的代码，其中这些方程在盆地上是有解的，并且其中代码通过对流体流动、沉积物压实和该组方程中的至少一个施加至少一种假设而被简化，以及用于求解该组方程以同时描述盆地中的流体流动和沉积物压实的代码。

此方面的实施可包括以下特征的一个或更多个。例如，计算机程序逻辑可包括用于在盆地的模型上提供网格的代码，其中网格包括多个单元。该组方程可包括描述第一方程的代码、描述第二方程的代码和描述第三方程的代码，第一方程定义每个单元的过压，第二方程定义每个单元的单元厚度，第三方程定义每个单元的物质负荷。

该计算机程序产品可包括用于向单元施加至少一种转换的代码，其中转换是沉积、下举、抬升和侵蚀中的一种。至少一种假设可包括第一假设、第二假设和第三假设，第一假设假设沉积物累积率是已知的，第二假设是压实仅发生在垂直方向，第三假设是压实是相对不可逆转的。流体可以是石油。

在另一一般方面中，用于在计算机上建模地表下地质盆地的方法包括接收定义盆地的至少一种物理特性的数据；定义描述盆地中流体流动和沉积物压实的一组方程，其中所述方程在物理区域上是有解的；通过施加压实仅发生在垂直方向的假设简化该组方程；以及求解该组方程以利用数据同时描述两种现象。

这方面的实施可包括以下特征的一个或更多个。该模型可涉及地表下的石油。该方法可进一步包括根据来自传感器的信息获得数据，所述传感器测量物理区域的至少一种物理特性。被求解的方程可帮助从盆地中开采石油。物理区域可以是地表下的地质盆地并且两种现象可以是盆地中的流体流动和物质压实，其中流体位于盆地中。可以基于该组被求解的方程产生地表下的地质盆地的盆地模型。可以基于盆地模型预测物理区域内碳氢化合物的位置。可以安排生产基础设施(例如，泵、压缩机和/或各种地表和地表下装置和设备)以在物理区域内基于预测的碳氢化合物的位置开采碳氢化合物。物理区域生产碳氢化合物的潜力可以基于盆地模型来评估。

以上内容相当广泛概述了本发明的特征和技术优势，以便更好地理解接下来对本发明的更详细的描述。本发明的附加特征和优势将在下文描述，其构成本发明权利要求的主题。本领域技术人员应当理解公开的概念和具体实施例可易于用作修改或设计用于实现本发明的相同目的的其他结构的基础。本领域技术人员还应当认识到这种等价结构并不偏离所附权利要求记载的本发明的精神和范围。结合附图和以下描述，将更好地理解被视为本发明特点的关于其操作组织和方法的新颖性特征以及本发明的另外目标和优势。但是，要清楚地理解，提供每幅附图的目的仅仅是为了图解说明和描述而不是要限定本发明。

附图说明

为了更完整地理解本发明，请结合附图参照以下描述，其中：

图1描述根据本发明的各实施例显示域内单元随时间积累而压实的模型示例；

图2描述根据本发明的各实施例通过沉积形成模型单元的示例；

图3是根据本发明的各实施例位于一层多层域内的单元的示例；

图4描述根据本发明的各实施例从一个域的一个单元移动到该域的另一个单元的流量示例；

图5描述根据本发明的各实施例用于建模物理区域的示例性方法；

图6描述适于是用本发明的计算机系统的方框图。

具体实施方式

本发明的实施例用于建模地表下的油田。本文描述的各实施例示例可参照此类油田。但是，实施例可用于建模涉及其他物质和/或过程的其他域。例如，各实施例可用于建模液体污染物在地下盆地中的分布、放射性物质自地下存储设备的迁移、或者其他液体、水、天然气或其他气体的迁移。

此类仿真中使用的数据可通过各种技术使用测量盆地的各种特性的传感器获得，各种技术例如是地层分析、地震反演或地质学家对这些内容的地质诠释。

接下来根据本发明的各实施例描述压实、松散和流体流动建模。这些模型优选考虑介质的力学平衡版。在说明书中，考虑了一组假设，其导致压实域中一般流体流动模型的公式表示。此外，说明书定义可施加到模型的简化假设，其减少了所需的计算。

质量、动量和本构关系的平衡

根据本发明的各实施例，可考虑沉积物和流体的物质平衡、力平衡以及流变本构关系以提供适当的盆地模型。该模型使用一般假设并使用特定考虑来简化建模过程。

地质盆地可表示为一组堆叠在一起的具有不同厚度的层。在盆地中的一些位置，层的厚度退化到零，形成尖灭(pinch-out)。为了在后面简化描述，盆地将被拓扑地视为一个平行六面体区或多个平行六面体区域，被称为单元。请注意，可替代使用根据2007年12月14日递交的题为″MODELING SUBSURFACE PROCESSES ONUNSTRUCTURED GRID″的美国专利申请61/007,761(代理机构卷号为2007EM361)中限定的实施例形成的棱柱形网格。

以下方程可表示单个平行六面体区域：

{(x，y，z；t)：0≤x≤X，0≤y≤Y，Z_top(x，y；t)≤z≤Z_bot(x，y；t)}，其中X和Y构成水平平面，Z_top(x，y；t)是沉积物的上层，而Z_bot(x，y；t)是沉积物的下层或基岩。

图1描述在计算域或区域104上压实过程的示例。在时刻t₁，域104具有顶部表面101和基底层103。感兴趣区域示为分区102。该区域可包括源岩。请注意，如图1所示，Z_top可以在地表、地表下的表面或海底。区域104以沉积速率q_s累积额外的沉积物，并且在时刻t₂，原始顶层101为地下层101′，并且该区域具有新的顶部层105。额外沉积物的重量已经导致感兴趣区102越来越深且越来越实，如区102’所示。底层103’也已经距离地表面移动得更深。包含在区域102’内的液体将经历压力增加，这会导致液体自区域102’排出。

请注意假设顶部表面101中的变化是已知的，即规定函数Z_top(x，y；t)。基岩的深度Z_bot(x，y；t)可以在每个点(x，y)以及每个时刻t处计算。计算域以曲线Z_top(x，y；t)和Z_bot(x，y；t)为边界，并且由于沉积物的沉积或侵蚀的原因随着时间扩大或收缩。沉积速率q_s可以是未知的，但是为了描述本发明的实施例，它是时间和空间的已知函数。

【0039】即便各实施例不受域维数的限制，但是以下段落假设压实展示在一维内(例如，垂直的)并且可以是非线性的。在以下内容中，z(t)将表示物质点在时刻t相对于海平面z＝0的坐标。请注意，同一物质点在时刻t’将具有不同的坐标z(t′)。坐标z＜0的负值表示海平面以上的高度。

压实的模型可以被看作是土壤加固的过程。沉积物充当可压缩的多孔基质。由于孔大小压实的原因，在时刻t₁占据体积Ω(t₁)的多孔岩石的元素将在时刻t₂占据体积Ω(t₂)并且具有相同的岩石基质密度和相同的质量，参见图1中的区102和102’。岩石质量守恒方程将具有以下形式：

$\frac{\∂(1-\mathrm{\φ}){\mathrm{\ρ}}_s}{\∂t}+\▿\·\left((1-\mathrm{\φ}){\mathrm{\ρ}}_s{\mathrm{\ν}}_r\right)=0,---\left(1.1\right)$

其中ρ_s是固岩质量密度，φ是孔隙度，而ν_r是岩石粒子速度。假设岩石是惰性的并且针对每种类型的沉积物具有恒定的岩石基质密度。在方程(1.1)右侧中的零意味着不考虑任何的固体物质体积源。沉积物的沉积可考虑作为边界条件。请注意，应当在考虑孔隙度对压力和有效应力的依赖关系之后单独描述侵蚀。

当考虑一维垂直压实时，岩石粒子速度是具有非零分量的向量，从而v_r＝(0，0，u_s)^T并且方程(1.1)变为

$\frac{\∂}{\∂t}\left({\mathrm{\ρ}}_s\mathrm{\φ}\right)-\frac{\∂}{\∂z}\left((1-\mathrm{\φ}){\mathrm{\ρ}}_s{u}_s\right)=0---\left(1.2\right)$

方程(1.2)的边界条件是通过岩石基质的沉积速率设置的。多孔岩石每次以已知的沉积速率q_s(t)≥0和已知的孔隙度φ₀(t)沉积。在短时间段Δt内，以下数量的岩石被添加到域

ΔM_rock＝a·Δt·q_s(t)，(1.3)

其中a是某点(x，y，Z_top(x，y；t))周围的小表面面积。

图2描述在图1的表层101上的沉积行为。图2显示当沉积物的部分沉积在域的顶表面时在t₁和t₂时刻的盆地顶层，其中t₂＝t₁+Δt。请注意部分A和B最初位于顶表面上，z坐标z(t₁)＝Z_top(t₁)，并且在沉积之后被掩埋并具有新的z坐标z(t₂)＞Z_top(t₂)。

由于岩石基质的密度是已知的，因此岩石的沉积量应该等于

ΔM_rock＝a·(z(t₂)-z(t₁))-(Z_top(t₂)-Z_top(t₁))·(1-φ₀(t))·ρ_s(z(t₁)).(1.4)

对于无穷小的时间段Δt，比较方程(1.3)和(1.4)并将极限视为Δt趋于0，则可获得域的顶部边界处关于物质点速度的以下表达式

${u_s|}_{Z_{\mathrm{top}}\left(t\right)}=\frac{\∂Z_{\mathrm{top}}\left(t\right)}{\∂t}+\frac{q_s\left(t\right)}{\mathrm{\ρ}\left(Z_{\mathrm{top}}\left(t\right)\right)\·(1-{\mathrm{\φ}}_0\left(t\right))}.---\left(1.5\right)$

由于函数Z_top(t)是已知的，因此其导数也是已知的，并且因此只要提供沉积速率q_s，则方程右侧是完全确定的。

在侵蚀的情况下，即自顶表面去除岩石，q_s＜0，岩石应当具有在压实过程中获得的孔隙度。在此情况下，方程(1.5)变为如下

${u_s|}_{Z_{\mathrm{top}}\left(t\right)}=\frac{\∂Z_{\mathrm{top}}\left(t\right)}{\∂t}+\frac{q_s\left(t\right)}{\mathrm{\ρ}\left(Z_{\mathrm{top}}\left(t\right)\right)\·(1-\mathrm{\φ}\left(Z_{\mathrm{top}}\left(t\right)\right))}.---\left(1.6\right)$

内部侵蚀的情形(例如从地下层去除岩石物质)应该以类似的方式处理。为了目前描述的目的，假定去除速率是已知的。

因此，边界条件变为相对孔隙度函数是非线性的。

对于xy-平面内点(x，y)周围的小面积ds考虑岩柱(column ofrock)

C(x，y；t)＝{ds×(Z_top(x，y；t)，Z_bot(x，y；t))}.

在任何固定时刻t，岩石在该岩柱中的总质量将由积分给定

$M_s(x,y;t)=\underset{C(x,y;t)}{\∫}(1-\mathrm{\φ}(x,y,z;t)){\mathrm{\ρ}}_s(x,y,z;t)\mathrm{dxdydz}.$

用面积ds细分该表达式的两部分并将极限视为ds趋于0提供

$M\left(t\right)=\underset{Z_{\mathrm{top}}\left(t\right)}{\overset{Z_{\mathrm{bot}}\left(t\right)}{\∫}}(1-\mathrm{\φ}(\mathrm{\ξ};t)){\mathrm{\ρ}}_s(\mathrm{\ξ};t)\mathrm{d\ξ}.---\left(1.7\right)$

该表达式对于任意点(x，y)成立，所以为了简化的目的可以忽略对(x，y)的依赖。

采用M(t)关于时间的实质导数(物质导数)并且使用方程(1.5)和(1.6)将导出以下方程

$\frac{D}{\mathrm{Dt}}M\left(t\right)=q_s\left(t\right).---\left(1.8\right)$

现在将方程(1.8)在时间间隔上求积分并代入方程(1.7)则可获得沉积物质量平衡的积分形式

$\underset{Z_{\mathrm{top}}\left(t\right)}{\overset{Z_{\mathrm{bot}}\left(t\right)}{\∫}}(1-\mathrm{\φ}(\mathrm{\ξ};t)){\mathrm{\ρ}}_s(\mathrm{\ξ};t)\mathrm{d\ξ}=\underset{0}{\overset{t}{\∫}}{q}_s\left(\mathrm{\τ}\right)\mathrm{d\τ}.---\left(1.9\right)$

该方案允许在之后的时刻t＞t₀确定在某时刻t₀＞0沉积的物质点的位置。考虑在时刻t₀位于域的顶表面的物质点，即具有垂直位置z(t₀)≡Z_top(t₀)。如果沉积速率为非零，则该点将被掩埋并且在时刻t＞t₀时，其位置将为z(t)＞Z_top(t)。考虑从Z_top(t)至z(t)的岩柱的质量平衡，可能获得以下等式：

$\underset{Z_{\mathrm{top}}\left(t\right)}{\overset{Z\left(t\right)}{\∫}}(1-\mathrm{\φ}(\mathrm{\ξ};t)){\mathrm{\ρ}}_s(\mathrm{\ξ};t)\mathrm{d\ξ}=\underset{t_0}{\overset{t}{\∫}}{q}_s\left(\mathrm{\τ}\right)\mathrm{d\τ}.---\left(1.10\right)$

利用方程(1.10)可以导出质量平衡的更通用形式。考虑在位置z(t₀)≥Z_top(t₀)某时刻t₀≥0的物质点。这些方程联系压实和流体流动。那么在之后的时刻t₁≥t₀，该点将具有由以下方程给出的位置z(t₁)：

$\underset{Z_{\mathrm{top}}\left(t_1\right)}{\overset{Z\left(t_1\right)}{\∫}}(1-\mathrm{\φ}(\mathrm{\ξ};t_1)){\mathrm{\ρ}}_s(\mathrm{\ξ};t_1)\mathrm{d\ξ}=\underset{Z_{\mathrm{top}}\left(t_0\right)}{\overset{z\left(t_0\right)}{\∫}}(1-\mathrm{\φ}(\mathrm{\ξ};t_0)){\mathrm{\ρ}}_s(\mathrm{\ξ};t_0)\mathrm{d\ξ}+\underset{t_0}{\overset{t_1}{\∫}}{q}_s\left(\mathrm{\τ}\right)\mathrm{d\τ}.---\left(1.11\right)$

为了简明，考虑单相流体流动的情形。对于用于确定盆地的沉积/侵蚀史以及前向压实过程的流体的物质平衡方程则具有以下形式：

$\frac{{\∂\mathrm{\ρ}}_a\mathrm{\φ}}{\∂t}+\▿\·\left({\mathrm{\ρ}}_a{\mathrm{\φ\ν}}_r\right)-\▿\·{\mathrm{\ρ}}_a\frac{K}{{\mathrm{\μ}}_a}(\▿p-{\mathrm{\ρ}}_{a}g{\▿}_z)=0,---\left(1.2\right)$

其中ρ_a是流体密度，μ_a是流体粘性，而K是渗透率。假设这些变量是已知的函数。

在引入压力势Φ＝p-ρ_agz之后，方程(1.12)可重写成

$\frac{{\∂\mathrm{\ρ}}_a\mathrm{\φ}}{\∂t}+\▿\·\left({\mathrm{\ρ}}_a{\mathrm{\φ\ν}}_r\right)-\▿\·{\mathrm{\ρ}}_a\frac{K}{{\mathrm{\μ}}_a}(\▿\mathrm{\Φ})=0.---\left(1.13\right)$

在底部边界或盆地基底103，可以假设无流动状态。在垂直边界上，如盆地顶表面101，可能具有无流动状态或流动状态边界。为了简单起见，将假定垂直边界具有无流动状态，但是，本发明的实施例可以具有流动状态。

以下示例的另一个假设是感兴趣的域低于海(地下水位)平面。这接着导致海平面下的岩石充满水的假设。换句话说，沉积的沉积物的孔体积包含水。沉积速率表示为q_a(x，y；t)。对于小时间段Δt期间点(x，y，Z_top(x，y；t))周围的小面积ds，以下水量将添加到盆地(请注意为了简明省略了(x，y))

$\mathrm{\ΔM}=\mathrm{ds}\·\mathrm{\Δt}\·q_a\left(t\right)=\mathrm{ds}\·\mathrm{\Δt}\·({u_s|}_{Z_{\mathrm{top}}\left(t\right)}-\frac{{\∂Z}_{\mathrm{top}}}{\∂t}\·{\mathrm{\φ}}_0\left(t\right)\·{\mathrm{\ρ}}_a\left(Z_{\mathrm{top}}\left(t\right)\right)).$

应用方程(1.5)产生

$q_a\left(t\right)=\frac{{\mathrm{\ρ}}_a\left(Z_{\mathrm{top}}\left(t\right)\right){\mathrm{\φ}}_0\left(t\right)}{{\mathrm{\ρ}}_a\left(Z_{\mathrm{top}}\left(t\right)\right)(1-{\mathrm{\φ}}_0\left(t\right))}{q}_s\left(t\right).---\left(1.14\right)$

在侵蚀的情形下，方程(1.14)变为：

$q_a\left(t\right)=\frac{{\mathrm{\ρ}}_a\left(Z_{\mathrm{top}}\left(t\right)\right){\mathrm{\φ}}_0(Z_{\mathrm{top}}\left(t\right);t)}{{\mathrm{\ρ}}_a\left(Z_{\mathrm{top}}\left(t\right)\right)(1-{\mathrm{\φ}}_0(Z_{\mathrm{top}}\left(t\right);t))}{q}_s\left(t\right).---\left(1.15\right)$

流体质量平衡在与移动沉积物连接的平行六面体单元(例如单元102)上的积分形式的推导开始于

C(t)＝{(x，y，z)：x₀≤x≤x₁，y₀≤y≤y₁，z₀(t)≤z≤z₁(t)}.

在任意固定的时刻t，该单元中流体的总质量由以下积分给出

$M_a\left(C\left(t\right)\right)=\underset{C\left(t\right)}{\∫}{\mathrm{\ρ}}_a\mathrm{\φd\Ω}.$

对于任意单元，其边框和物质点一起移动

$\frac{{\∂z}_0}{\∂t}={u_s|}_{z_0}$ 和 $\frac{{\∂z}_1}{\∂t}={u_s|}_{z_1}.---\left(1.16\right)$

将方程(1.13)和(1.16)组合产生

$\frac{D}{\mathrm{Dt}}{M}_a\left(C\left(t\right)\right)-\underset{C\left(t\right)}{\∫}\▿\·({\mathrm{\ρ}}_a\frac{K}{{\mathrm{\μ}}_a}\▿\mathrm{\Φ})\mathrm{d\Ω}=\underset{C\left(t\right)}{\∫}{q}_a\mathrm{d\Ω}.---\left(1.17\right)$

对于邻近顶部边界101的任意单元，例如如果单元102的上表面包括表面101的一部分，则使用以下方程

$C_{\mathrm{top}}\left(t\right)=\{(x,y,z):x_0\≤x\≤x_1,y_0\≤y\≤y_1,Z_{\mathrm{top}}\left(t\right)\≤z\≤z_1\left(t\right)\},$

使用方程(1.5)和(1.6)代替方程(1.16)提供时间导数如下

$\frac{\∂Z_{\mathrm{top}}}{\∂t}={u_s|}_{Z_{\mathrm{top}}\left(t\right)}-\frac{q_s\left(t\right)}{{\mathrm{\ρ}}_s(1-\stackrel{\‾}{\mathrm{\φ}})}{|}_{Z_{\mathrm{top}}}$ 和 $\frac{{\∂z}_1}{\∂t}={u_s|}_{z_1},---\left(1.18\right)$

其中，对于沉积

或对于侵蚀

对于这种单元，方程(1.17)应当修改如下

$\frac{D}{\mathrm{Dt}}{M}_a\left(C\left(t\right)\right)-\underset{C\left(t\right)}{\∫}\▿\·({\mathrm{\ρ}}_a\frac{K}{{\mathrm{\μ}}_a}\▿\mathrm{\Φ})\mathrm{d\Ω}=\underset{C\left(t\right)}{\∫}{q}_a\mathrm{d\Ω}+\underset{y_0}{\overset{y_1}{\∫}}\underset{x_0}{\overset{x_1}{\∫}}\frac{{\mathrm{\ρ}}_a\stackrel{\‾}{\mathrm{\φ}}{q}_s\left(t\right)}{{\mathrm{\ρ}}_s(1-\stackrel{\‾}{\mathrm{\φ}})}{|}_{Z_{\mathrm{top}}}\mathrm{ds},---\left(1.19\right)$

其中最后的积分表示由于沉积或侵蚀过程而分别自系统添加或移除的流体质量。

对于多孔介质中的流体流动，总动量方程可以写成

其中

是应力张量。容积密度ρ是利用容积分数加权的成分的密度之和，具体如下：

ρ＝ρ_s(1-φ)+ρ_aφ.           (1.21)

应力张量可以以

的形式考虑，其中引入负号是为了与岩石力学用法一致。那么方程(1.20)可以以另一形式表示

$\frac{\∂\mathrm{\σ}}{\∂z}=(({\mathrm{\ρ}}_{s}1-\mathrm{\φ})+{\mathrm{\ρ}}_a\mathrm{\φ})g.---\left(1.22\right)$

有效应力σ_E和静岩负荷L可以分别表示为应力σ与流体孔压p以及σ和流体静压p_h之间的差

σ_E＝σ-p和L＝σ-p_h.   (1.23)

利用压力势的定义，有效应力的另一形式为

σ_E＝L-Φ.

因此，力平衡方程(1.22)可以用σ_E和L表示如下

$\frac{{\∂\mathrm{\σ}}_E}{\∂z}=\frac{\∂(L-\mathrm{\Φ})}{\∂z}.---\left(1.24\right)$

对于可压缩流体，任意点的流体静压可以由以下方程给出

$p_h(z;t)=p\left(Z_{\mathrm{top}}\left(t\right)\right)+g\underset{Z_{\mathrm{top}}\left(t\right)}{\overset{z}{\∫}}{\mathrm{\ρ}}_a\left(p\left(\mathrm{\ξ}\right)\right)\mathrm{d\ξ}.---\left(1.25\right)$

将方程(1.22)和(1.25)结合，可将静岩负荷写成

$L(z;t)=g\underset{Z_{\mathrm{top}}\left(t\right)}{\overset{z}{\∫}}({\mathrm{\ρ}}_s-{\mathrm{\ρ}}_a\left(p\left(\mathrm{\ξ}\right)\right))(1-\mathrm{\φ})\mathrm{d\ξ}.---\left(1.26\right)$

基于Athy in L.Athy在″Density，porosity，and compaction ofsedimentary rocks″，BuI.Am.Assoc.Geol.，14(1930)，pp.1-24中对沉积盆地的实验观察，提出孔隙度φ和深度z之间存在直接关系。以其最简单的形式，这种关系可表示为

φ＝φ₀e^-bz.              (1.27)

观察发现孔隙度是有效应力σ_E的函数φ＝φ(σ_E)，并且正是借助有效应力σ_E对常压下沉积物的深度的依赖性能够推断出诸如方程(1.27)中阐明的关系。例如，参见P.Allen和J.Allen的文章″Basin Analysis，Principles and Applications″，Blackwell Scientific Publications，马萨诸塞州剑桥，1990，其通过引用整体并入本文。因此，尽管针对常压下岩石的孔隙度-深度关系看似很强，但是φ与z之间的关系的推论只是有益的工具。换句话说，孔隙度和负荷在每个点都是有关联的。在本发明的各实施例中，孔隙度被视为有效应力的函数。请注意本发明的其他实施例可使用其他类型的流变学。而且，本构孔隙度-有效应力关系可以以双指数的形式假设为

$\mathrm{\φ}={\mathrm{\φ}}_c+{\mathrm{\φ}}_1{e}^{-b_1{\mathrm{\σ}}_E}+{\mathrm{\φ}}_2{e}^{-b_2{\mathrm{\σ}}_E},---\left(1.28\right)$

其中φ_c是孔隙度的下限(不可约)孔隙度，而φ_c+φ₁+φ₂是地表条件的沉积物的孔隙度。

通常，沉积物随着时间的推移被掩藏又不随着时间的推移被发掘，因此模型中的应力随着时间趋于增长。但是，在考虑侵蚀的模型中，有效应力σ_E可能降低。在此情形下，孔隙度允许根据以下表达式出现微小的增长

$\mathrm{\φ}={\mathrm{\φ}}_c+({\mathrm{\φ}}_0-{\mathrm{\φ}}_c)e^{-b{\mathrm{\σ}}_E-b_{\mathrm{ul}}({\mathrm{\σ}}_E^{\mathrm{new}}-{\mathrm{\σ}}_E)},---\left(1.29\right)$

其中，是同一物质点处新的、降低的有效应力，而b_ul是去荷可压缩性(unloading compressibility)。

为了考虑压实不可逆转的性质以及在有效应力由于侵蚀的原因降低时允许小的松散，孔隙度被假设为两个变量的时间依赖函数，两个变量也就是任意给定时刻t的有效应力和模型所有之前寿命内达到的应力的历史最大值，并且可以表示为

$\mathrm{\φ}\left(z\left(t\right)\right)=\mathrm{\φ}({\mathrm{\σ}}_E\left(z\left(t\right)\right),{\mathrm{\σ}}_E^{\max }\left(z\left(t\right)\right)),$

其中，

而z(t)是物质点在时间t的z坐标，函数σ_E(z)由方程(1.23)定义。

如果在任意给定的时间有效应力变得小于其历史最大值，则应用方程(1.29)计算孔隙度。否则使用方程(1.28)。

全耦合压力模型

基于以上部分描述的质量、动量和本构关系的平衡，通过以下方程组描述了压实域中的单相流体流动。请注意，在每个特定单元处存在四个要求解的未知量，分别为孔隙度φ(z(t))、压力势Φ(z；t)、静岩负荷L(z；t)和流体静压p_h(z；t)。

方程组2.1：

$\frac{{\∂\mathrm{\ρ}}_a\mathrm{\φ}}{\∂t}+\▿\·\left({\mathrm{\ρ}}_a\mathrm{\φ}{u}_s\right)-\▿\·{\mathrm{\ρ}}_a\frac{K}{{\mathrm{\μ}}_a}\▿\mathrm{\Φ}=q_a,$

$\frac{\∂}{\∂t}\left({\mathrm{\ρ}}_s\mathrm{\φ}\right)-\frac{\∂}{\∂z}\left((1-\mathrm{\φ}){\mathrm{\ρ}}_s{u}_s\right)=0,$

$\frac{\∂L}{\∂z}=({\mathrm{\ρ}}_s-{\mathrm{\ρ}}_a)(1-\mathrm{\φ})g,$

σ_E＝L-Φ，

${\mathrm{\σ}}_E^{\max }\left(z\left(t\right)\right)=\underset{\mathrm{\τ}\≤t}{\sup }\left\{{\mathrm{\σ}}_E^{\max }\left(z\left(\mathrm{\τ}\right)\right)\right\},$

$\mathrm{\φ}\left(z\left(t\right)\right)=\mathrm{\φ}({\mathrm{\σ}}_E\left(z\left(t\right)\right),{\mathrm{\σ}}_E^{\max }\left(z\left(t\right)\right)),$

$p_h(z;t)=p\left(Z_{\mathrm{top}}\left(t\right)\right)+g\underset{Z_{\mathrm{top}}\left(t\right)}{\overset{z}{\∫}}{\mathrm{\ρ}}_a(p_h+\mathrm{\Φ}))\mathrm{d\ξ},$

${u_s|}_{Z_{\mathrm{top}}\left(t\right)}=\frac{\∂Z_{\mathrm{top}}\left(t\right)}{\∂t}+\frac{q_s\left(t\right)}{{\mathrm{\ρ}}_s(1-\mathrm{\φ})}{|}_{Z_{\mathrm{top}}},$

p(Z_top(t))＝p_atm+ρ_seag·max{0，Z_top}，

(x，y，z(t))∈{(x，y，z；t)：0≤x≤X，0≤y≤Y，Z_top(x，y；t)≤z≤Z_bot(x，y；t)}，

其中，Z_top(x，y；t)是盆地顶表面101(或海底)，Z_bot(x，y；t)是图1中的盆地基底103。因此，只要规定了沉积速率q_s(x，y；t)，方程系统，方程组2.1，是完全确定的。

在遵循盆地地层学的曲线坐标系统内考虑方程系统(以上方程组2.1所定义的)。换句话说，x和y方向沿地层时间线延伸并且因此弯曲而遵循盆地区域的倾斜。这种约束保持坐标系统的轴沿着最大渗透率的方向(渗透率椭圆体的主轴)，该方向在断裂盆地地层中通常沿着地层排列。

z方向被处理为好像与x垂直，但z方向实际上沿着垂直方向延伸。此方位是正向下，原点位于盆地顶表面或海平面。除了考虑水平产状沉积物时之外，坐标系统并非真正正交的事实将引入计算上的误差。在典型的盆地地层的倾斜处，该误差相当小，尤其是当与坐标系统为正交但是相对渗透率椭圆体的轴是倾斜的时引入的误差相比。

本发明的实施例假设渗透介质具有分层的结构并且每层具有相同的属性。换句话说，来自方程(1.28)的系数φ_c、φ₁、φ₂、b₁和b₂以及来自方程(1.21)的岩石密度ρ_s被假设为是分段常数。因此，如果与表面点(x，y)对应的每柱被视为分成n_z层，使得

$z_0\≡Z_{\mathrm{top}}\≤z_1\≤...\≤Z_{n_z-1}\≤z_{n_z}\≡z_{\mathrm{bot}},$

则在任何时刻，在每层中系数和

为常数l＝1，...，n_z。

在本发明的其他实施例中，假设盆地的发展是从过去的某时刻T_s＜0到现在时刻T_e＝0被建模。列举了从顶部到底部的多个层。每层的起始和停止沉积时刻分别表示为t_sl和t_el。这导致每个第l层在第l-1层之前沉积的假设，从而

$T_s=t_{{\mathrm{sn}}_z}<t_{{\mathrm{en}}_z}<t_{{\mathrm{sn}}_z-1}<...<t_{e2}<t_{s1}<t_{e1}<T_e.---\left(2.2\right)$

本发明的实施例使用拉格朗日方案导出离散化。因此，网格跟随移动的沉积物。根据本发明的实施例，计算网格通过以下的方式构造。首先，在xy平面内构造网格。之后，网格被垂直拉伸以形成列。为了简明的目的，假设网格为矩形。但是，xy-网格可以是不均匀的，并且x和y方向的网孔尺寸可以是任意的。因此，在xy平面内构造矩形网格，使得

$x_0=0<x_1<...<x_{n_{x-1}}<x_{n_x}=X,y_0=0<y_1<...<y_{n_{y-1}}<y_{n_y}=Y,$

并且n_x×n_y列通过以下表达式限定

Col_i，j(t)＝{(x，y，z；t)：x_i-1≤x≤x_i，y_j-1≤y≤y_j，Z_top(t)≤z≤Z_bot(t)}.

根据本发明的实施例，可以既在整个列集合上又在这些列的子集上甚至在单个列上实施计算。

每列具有相同数量的层n_z并且其中一些层在域的一部分中具有零厚度，这表示已经在xy平面的该部分修剪过特定层。启动仿真的一种方式是将计算域设为零厚度，即Z_top(x，y；T_s)＝Z_bot(x，y；T_s)。其他方式可能具有非零厚度，从而已经存在一个或更多层。

优选将总的时间间隔[T_s，T_e]分割成M个较小的间隔Δt＝t_n-t_n-1，T_s＝t₀＜...＜t_M＝T_e，分割方式使得针对来自方程(2.2)的每个t_ej(或t_sj)存在使t_i＝t_ej的下标i。

当计算过程从一个时步移动到下一时步[t_n-1，t_n]时，本发明的实施例假设时步开始时(即在时刻t_n-1)的计算几何体系是已知的。推理地，在时刻t_n时单元的厚度是未知的并且应当是仿真的一部分。由于利用本发明的实施例，岩石运动发生在垂直方向，因此单元优选被视为厚度可以随时间变化的平行六面体。图3图示了计算单元301的一个示例。单元301位于由x_i-1和x_i限定的列内。如图3所示，由于层k可具有位于单元301之上的单元和位于单元301之下的单元，所以每层在一列内可具有一个以上单元。计算单元可以表示为

C_i，j，k(t)＝{(x，y，z；t)：x_i-1≤x≤x_i，y_j-1≤y≤y_j，z_i，j，k-1(t)≤z≤z_i，j，k(t)}，

其中k＝1，...，n_z。在以下论述中，为了简明的目的，单元可通过一个下标而不是三个一组的下标来标注。换句话说，单元301可以利用下标k标注为单元C_k而不使用三个一组的下标i，j，k得到的标记C_i，j，k。

根据本发明的实施例，在仿真开始时，每个单元起源于域的顶部。随着沉积物的沉积，单元随着时间而发展。之后，当完全沉积时，由于新单元沉积在该单元顶部，该单元接着被掩埋并且被压实。在不存在成岩作用的情况下，被完全沉积之后的任何单元保持恒定的岩石质量，除非该单元通过上层单元的侵蚀移动到其开始受到侵蚀的表面。

可以向任何计算单元应用不同类型的转换。一种类型是沉积，借此将单元沉积在域的顶表面。单元随着时间的推移而发展，岩石质量增加并且孔隙度可能变化。另一类型是下举(downlift)，借此由于新的单元沉积在单元的顶部而掩埋并压实该单元。该单元的岩石质量不改变，而该单元的孔隙度通常可能减小。另一类型是隆起或抬升(uplift)，借此由于海底抬升或上面单元的侵蚀，单元在列中被提升。单元的岩石质量不变化，而单元的孔隙度可能稍微增长。另一类型是侵蚀，其中单元受到侵蚀。当单元局部或全部被侵蚀时，单元的岩石质量减小，而孔隙度可能稍微增大。

由于任何单元C_i，j，k(t_n)的孔隙度会随时间而改变，因此单元的厚度也会随时间而改变，如下式表示的

$\mathrm{\&Delta;}{z}_{i,j,k}^n=z_{i,j,k}^n-z_{i,j,k-1}^n,$

因此，在时间t_n的计算网格是明确未知的并且应当是仿真的一部分。

方程组2.1的第一个方程采用关于流体静力负荷的超压Φ书写。超压用作主变量并且在整个计算单元中被视为常数，因此其值与单元中心有关。之后，总的孔压力将表示成流体静压和超压的和，如p_i，j，k＝p_{h；i，j，k}+Φ_i，j，k表示的。

本发明的实施例利用有限体积方法将孔隙度离散化，其中孔隙度的离散值是整个单元的平均孔隙度，如下式表示的

${\mathrm{\&phi;}}_{i,j,k}=\frac{1}{V_{i,j,k}}\underset{C_{i,j,k\left(t\right)}}{\&Integral;}\mathrm{\&phi;}(x,y,z)\mathrm{d\&Omega;}=\frac{1}{{\mathrm{\&Delta;z}}_{i,j,k}}\underset{z_{i,j,k-1}}{\overset{z_{i,j,k}}{\&Integral;}}\mathrm{\&phi;}(x,y,z)\mathrm{dz},$

其中V_i，j，k是单元的体积。

设s_i，j，k表示单元固体厚度，其值为单元的总压缩的岩石体积除以单元的水平表面积，如下式表示的

$s_{i,j,k}=\frac{1}{{\mathrm{\&Delta;x}}_i{\mathrm{\&Delta;y}}_j}\underset{C_{i,j,k}\left(t\right)}{\&Integral;}(1-\mathrm{\&phi;}(x,y,z))\mathrm{d\&Omega;},$

其中Δx_i和Δy_j分别是单元在x和y方向的尺寸。在不存在成岩作用的情况下，根据方程组(2.1)的第二个方程，可以断定固体厚度的值可以表示成

$s_{i,j,k}=\underset{z_{i,j,k-1}}{\overset{z_{i,j,k}}{\&Integral;}}(1-\mathrm{\&phi;}\left(z\right))\mathrm{dz}=(1-{\mathrm{\&phi;}}_{i,j,k}){\mathrm{\&Delta;z}}_{i,j,k},---\left(2.3\right)$

并且在单元C_i，j，k被完全沉积之后该值不随时间而变化。如果单元的沉积历史的起始时间t_sk和终止时间t_ek以及沉积速率q_{s；i，j，k}是已知的，则在t_sk之后的任意时刻单元的固体厚度可以由下式确定

$s_{i,j,k}\left(t\right)=(t_{\mathrm{ek}}-t_{\mathrm{sk}})\frac{q_{s;i,j,k}}{{\mathrm{\&rho;}}_{s;i,j,k}}\min \{1,\frac{{t-t}_{\mathrm{sk}}}{t_{\mathrm{ek}}-t_{\mathrm{sk}}}\}.$

反之亦然，如果层k中单元的固体厚度s_i，j，k及其沉积的起始和终止时间t_sk和t_ek是已知的，则层的沉积速率可利用下式计算

$q_{s;i,j,k}=\frac{s_{i,j,k}{\mathrm{\&rho;}}_{s;i,j,k}}{t_{\mathrm{ek}}-t_{\mathrm{sk}}}.$

表达式(2.3)提供一种根据给定的单元的固体厚度和能渗透的厚度计算平均孔隙度的方法

${\mathrm{\&phi;}}_{i,j,k}=1-\frac{s_{i,j,k}}{{\mathrm{\&Delta;z}}_{i,j,k}}.---\left(2.4\right)$

由于引入固体厚度，单个单元上累积的静岩负荷可以用以下式子表示

${\mathrm{\&Delta;L}}_{i,j,k}=\underset{z_{i,j,k-1}}{\overset{z_{i,j,k}}{\&Integral;}}({\mathrm{\&rho;}}_s-{\mathrm{\&rho;}}_a)(1-\mathrm{\&phi;})\mathrm{dz}=g({\mathrm{\&rho;}}_{s;i,j,k}-{\mathrm{\&rho;}}_{a;i,j,k})s_{i,j,k}.---\left(2.5\right)$

在不可压缩流体的情形中，则仿真过程中流体密度不会改变，并且因此可以表示为：

${\mathrm{\&rho;}}_{a;i,j,k}={\mathrm{\&rho;}}_a^0.---\left(2.6\right)$

在可压缩流体的情形中，应当考虑流体密度对孔压力的依赖性，孔压力表示为静压头p_h和超压Ω的和。由于为了计算的目的，每个单元内的孔压被视为常数，由此水密度在单元内也被视为常数，并且由单元上的压力值定义。在此情形下，流体密度可以表示为：

${\mathrm{\&rho;}}_{a;i,j,k}={\mathrm{\&rho;}}_a^0\&CenterDot;e^{\mathrm{\&beta;}(p_{h;i,j,k}+{\mathrm{\&Phi;}}_{i,j,k}-p_{\mathrm{atm}})}.---\left(2.7\right)$

注意，单个单元上累积的流体静压将具有以下形式：

$\mathrm{\&Delta;}{p}_{h;i,j,k}=\underset{z_{i,j,k-1}}{\overset{z_{i,j,k}}{\&Integral;}}g{\mathrm{\&rho;}}_a\mathrm{dz}=g{\mathrm{\&rho;}}_{a;i,j,k}{\mathrm{\&Delta;z}}_{i,j,k},$

并且单元C_i，j，k中心的流体静压p_{h；i，j，k}的值可以计算如下：

$p_{h;i,j,1}=p_{h;i,j,\mathrm{surf}}+\frac{1}{2}\mathrm{\&Delta;}{p}_{h;i,j,1},---\left(2.8\right)$

$p_{h;i,j,k}=p_{h;i,j,k-1}+\frac{1}{2}({\mathrm{\&Delta;p}}_{h;i,j,k-1}+p_{h;i,j,k}),k=2,...,n_2,$

其中p_{h；i，j，surf}是列C_i，j，k的顶表面处流体静压的值。

如上所述，仿真时网格并非确切地已知并且应当是计算的一部分。单元厚度依赖于单元顶部的被掩埋的沉积物量和超压的值。方程组(2.1)的第三方程用于获得单元厚度的该组离散方程。将两部分都除以方程右侧以及从

到

进行积分提供

通常，右手侧的积分不能以解析的方式计算，替代地可以近似得到。由于孔隙度-有效应力关系的指数形式的原因，一点和两点近似可能无法提供良好的准确性。三点辛普森(Simpson)公式可为不是很厚(例如，≤1km)的计算单元的积分提供良好的近似。在特定情形，例如，厚的单元或高度变化的孔隙度关系，可能不得不用多点求积以在近似方程(2.9)中的积分。下面的讨论仅仅使用Simpson规则作为举例，因为可以使用其他近似。因此，使用方程(2.5)，方程(2.9)的近似值变成以下表达式(请注意，为了简明省略下标i和j)

其中

$L_{i,j,k}^{\mathrm{top}}\&equiv;L_{i,j,k-1/2}=L_{i,j,k}-\frac{1}{2}\mathrm{\&Delta;}{L}_{i,j,k},$

$L_{i,j,k}^{\mathrm{bot}}\&equiv;L_{i,j,k+1/2}=L_{i,j,k}+\frac{1}{2}\mathrm{\&Delta;}{L}_{i,j,k},---\left(2.11\right)$

并且L_i，j，k是单元C_i，j，k的中心处静岩负荷的值，其计算如下

$L_{i,j,1}=\frac{1}{2}{\mathrm{\&Delta;L}}_{i,j,1,}$

$L_{i,j,k}=L_{i,j,k-1}+\frac{1}{2}({\mathrm{\&Delta;L}}_{i,j,k-1}+{\mathrm{\&Delta;L}}_{i,j,k}),k=2,...,n_z.---\left(2.12\right)$

超压在单元边界

和

的值在以下段落中提供。

方程组(2.1)的第一个方程优选利用有限体积方法离散化，该有限体积方法可以通过以下方式应用。第一个方程在计算块(例如C^t)上并且在时步[t_n-1，t_n]内积分。请注意，每个计算块与物质坐标有关，因此时间上以某一速率v_r移动。应用散度定理(divergence theorem)和在时步内积分方程(1.17)得到

$\underset{t_{n-1}}{\overset{t_n}{\&Integral;}}\frac{D}{\mathrm{Dt}}{M}_a\left(C^t\right)\mathrm{dt}-\underset{t_{n-1}}{\overset{t_n}{\&Integral;}}\underset{{\&PartialD;C}^t}{\&Integral;}(n,{\mathrm{\&rho;}}_a\frac{K}{{\mathrm{\&mu;}}_a}\&dtri;\mathrm{\&Phi;})\mathrm{dsdt}=\underset{t_{n-1}}{\overset{t_n}{\&Integral;}}\underset{C^t}{\&Integral;}{q}_a\mathrm{d\&Omega;dt}.$

请注意，左手侧中的第一项可以明确地对时间积分以形成

$(M_a\left(C^{t_n}\right)-M_a\left(C^{t_{n-t}}\right))\underset{t_{n-1}}{\overset{t_n}{\&Integral;}}\underset{\&PartialD;C^t}{\&Integral;}(n,{\mathrm{\&rho;}}_a\frac{K}{{\mathrm{\&mu;}}_a}\&dtri;\mathrm{\&Phi;})\mathrm{dsdt}=\underset{t_{n-1}}{\overset{t_n}{\&Integral;}}\underset{C^t}{\&Integral;}{q}_a\mathrm{d\&Omega;dt}.---\left(2.13\right)$

由于单元C_i，j，k中的流体质量由下式提供

$M_a\left(C_{i,j,k}\right)=\underset{C_{i,j,k}}{\&Integral;}{\mathrm{\&rho;}}_a\mathrm{\&phi;}(x,y,z)\mathrm{dxdydz},$

其在时刻t_n-1被近似为

$M_a\left(C_{i,j,k}^{t_{n-1}}\right)\&ap;{\mathrm{\&Delta;x}}_i{\mathrm{\&Delta;y}}_j{\mathrm{\&Delta;z}}_{i,j,k}^{n-1}{\mathrm{\&rho;}}_{a;i,j,k}^{n-1}{\mathrm{\&phi;}}_{i,j,k}^{n-1}.---\left(2.14\right)$

而利用方程(2.4)其在时刻t_n可近似为

$M_a\left(C_{i,j,k}^{t_n}\right)\&ap;{\mathrm{\&Delta;x}}_i{\mathrm{\&Delta;y}}_j{\mathrm{\&rho;}}_{a;i,j,k}^n({\mathrm{\&Delta;z}}_{i,j,k}^n-s_{i,j,k}^n).---\left(2.15\right)$

如果没有至单元的流体的内部源，则函数q_a为零并且只有流体增加或移除通过沉积或侵蚀过程。利用方程(1.19)，方程(2.13)中右手侧变为

$\underset{t_{n-1}}{\overset{t_n}{\&Integral;}}\underset{C^t}{\&Integral;}{q}_a\mathrm{d\&Omega;dt}=\mathrm{\&Delta;x\&Delta;y}{\mathrm{\&rho;}}_{a;i,j,k}^{*}{\mathrm{\&phi;}}_{i,j,k}^{*}\frac{q_{s;i,j,k}\left(t_{n-1}\right)}{{\mathrm{\&rho;}}_{s;i,j,k}(1-{\mathrm{\&phi;}}_{i,j,k}^{*})},---\left(2.16\right)$

其中，标记^*表示该值分别从表面(输入数据)或从针对沉积或侵蚀的之前时步t_n-1取得。

请注意，每个计算块

是平行六面体的形式，所述平行六面体的面与坐标平面平行。因此，(2.13)左手侧中的表面积分项可以通过以下表达式近似得到

$\underset{t_{n-1}}{\overset{t_n}{\&Integral;}}\underset{\&PartialD;C_{i,j,k}^t}{\&Integral;}(n,{\mathrm{\&rho;}}_a\frac{K}{{\mathrm{\&mu;}}_a}\&dtri;\mathrm{\&Phi;})\mathrm{dsdt}\&ap;\mathrm{\&Delta;}{t}_n\underset{m=1}{\overset{6}{\mathrm{\&Sigma;}}}{\left(\underset{\&PartialD;C_{i,j,k;m}^{*}}{\&Integral;}(n_m,{\mathrm{\&rho;}}_a\frac{K}{{\mathrm{\&mu;}}_a}\&dtri;\mathrm{\&Phi;})\right)}^{*},---\left(2.17\right)$

其中数量(-)^*代表在时间上积分的某一近似值，而表示平行六面体的第m面。为了产生完全的隐公式，应当在t^*＝t_n时考虑所有参数并且方程(2.17)变为

$\underset{t_{n-1}}{\overset{t_n}{\&Integral;}}\underset{{\&PartialD;C}_{i,j,k}^t}{\&Integral;}(n,{\mathrm{\&rho;}}_a\frac{K}{{\mathrm{\&mu;}}_a}\&dtri;\mathrm{\&Phi;})\mathrm{dsdt}\&ap;{\mathrm{\&Delta;t}}_n\underset{m=1}{\overset{6}{\mathrm{\&Sigma;}}}{\left(\underset{\&PartialD;C_{i,j,k;m}^{t_n}}{\&Integral;}(n_m,{\mathrm{\&rho;}}_a\frac{K}{{\mathrm{\&mu;}}_a}\&dtri;\mathrm{\&Phi;})\right)}^{\left(n\right)}.---\left(2.18\right)$

来自方程2.18的面积分

的近似被定义如下。

考虑通量

的法向分量的面积分，其被表示成

${\left({\mathrm{Flux}|}_{{\&PartialD;C}_m^{*}}\right)}^{(*)}={\left(\underset{{\&PartialD;C}_m^{*}}{\&Integral;}(n_m,{\mathrm{\&rho;}}_a\frac{K}{{\mathrm{\&mu;}}_a}\&dtri;\mathrm{\&Phi;})\right)}^{*}---\left(2.19\right)$

方程(2.19)的近似值的例子如图4中所示，其表示从单元C_i，j，k403到单元C_i+1，j，k403的通量401。请注意，通量401在x方向并且自单元402的中心发源并移动到单元403的中心。单元402和403的x面的面积分别表示为S_x，i和S_x，i+1。请注意，立方体C_i具有六个侧面，其中的一个侧面S_x，i邻接立方体C_i+1，请参见列出方程2.18的那段。

Φ(r_i+1，j，k)和Φ(r_i，j，k)之间的差可通过积分表示为

$\mathrm{\&Phi;}\left(r_{i+1,j,k}\right)-\mathrm{\&Phi;}\left(r_{i,j,k}\right)=\underset{r_{i,j,k}}{\overset{r_{r+1,j,k}}{\&Integral;}}(\&dtri;\mathrm{\&Phi;},{\stackrel{\&RightArrow;}{\mathrm{\&tau;}}}_x)\mathrm{dl}.$

移动性可以表示为

并且指示w＝aK▽Φ。则

并且积分可以重写为

$\underset{r_{i,j,k}}{\overset{r_{i+1,j,k}}{\&Integral;}}(\&dtri;\mathrm{\&Phi;},{\stackrel{\&RightArrow;}{\mathrm{\&tau;}}}_x)\mathrm{dl}=\underset{r_{i,j,k}}{\overset{r_{i+1,j,k}}{\&Integral;}}\frac{1}{a}(K^{-1}w,{\stackrel{\&RightArrow;}{\mathrm{\&tau;}}}_x)\mathrm{dl}=\underset{r_{i,j,k}}{\overset{r_{i+1,j,k}}{\&Integral;}}\frac{1}{a}(w,K^{-1}{\stackrel{\&RightArrow;}{\mathrm{\&tau;}}}_x)\mathrm{dl}.$

由于渗透率张量K是在和层结构对齐的坐标系统中的对角线，那么向量

是K的特征向量，即因此上述差可以表示为

$\mathrm{\&Phi;}\left(r_{i+1,j,k}\right)-\mathrm{\&Phi;}\left(r_{i,j,k}\right)=\underset{r_{i,j,k}}{\overset{r_{i+1,j,k}}{\&Integral;}}\frac{1}{{\mathrm{ak}}_x}(w,{\stackrel{\&RightArrow;}{\mathrm{\&tau;}}}_x)\mathrm{dl}.$

通过相同的方式，每个单元中的通量可以被视为就好像引入了单元的公共面上在点r₀404的势能，其中

${\mathrm{\&Phi;}}^{(i,j,k)}\left(r_0\right)-\mathrm{\&Phi;}\left(r_{i,j,k}\right)=\underset{r_{i,j,k}}{\overset{r_0}{\&Integral;}}\frac{1}{{\mathrm{ak}}_x}(w,{\stackrel{\&RightArrow;}{\mathrm{\&tau;}}}_x)\mathrm{dl},$

$\mathrm{\&Phi;}\left(r_{i+1,j,k}\right)-{\mathrm{\&Phi;}}^{(i+1,j,k)}\left(r_0\right)=\underset{r_0}{\overset{r_{i+1,j,k}}{\&Integral;}}\frac{1}{{\mathrm{ak}}_x}(w,{\stackrel{\&RightArrow;}{\mathrm{\&tau;}}}_x)\mathrm{dl}.$

之后，这些积分可以通过以下表达式来近似得到

其中，

是通量分量沿着向量在单元r_i，j，k中心的值。由于系数a和k_x与其在计算块中心处的值有关，因此在下文中将它们作为a_α，j，k＝a(r_α，j，k)和k_{x，α，j，k}＝k_x(r_α，j，k)，α＝i，i+1.提及。

由于相邻单元的同一表面的势能Φ^(i，j，k)(r₀)和Φ^{(i+1，j，k)}(r₀)的值一致并且从单元C_i，j，k流出的通过面S_x，i的通量的值等于通过面S_x，i+1流入单元C_i+1，j，k的通量的值，即

Φ^(i，j，k)(r₀)＝Φ^{(i+1，j，k)}(r₀)≡Φ₀

和

因此，得到副势能(auxiliary potential)Φ₀的值是可能的：

${\mathrm{\&Phi;}}_0=\frac{\mathrm{\&Phi;}\left(r_{i,j,k}\right)\frac{a_{i,j,k}{k}_{x,i,j,k}{S}_{x,i}}{\mathrm{\&Delta;}{x}_i}+\mathrm{\&Phi;}\left(r_{i+1,j,k}\right)\frac{a_{i+1,j,k}{k}_{x,i+1,j,k}{S}_{x,i+1}}{\mathrm{\&Delta;}{x}_{i+1}}}{\frac{a_{i,j,k}{k}_{x,i,j,k}{S}_{x,i}}{\mathrm{\&Delta;}{x}_i}+\frac{a_{i+1,j,k}{k}_{x,i+1,j,k}{S}_{x,i+1}}{\mathrm{\&Delta;}{x}_{i+1}}}.---\left(2.21\right)$

由于从单元C_i，j，k到单元C_i+1，j，k的通量可以利用下式计算，

${\mathrm{Flux}}_{i,j,k}^{i+1,j,k}=\underset{S_i}{\&Integral;}(n_x,w)\mathrm{ds}\&ap;({\mathrm{\&Phi;}}_0-\mathrm{\&Phi;}\left(r_{i,j,k}\right))\frac{a_{i,j,k}{k}_{x,i,j,k}{S}_{x,i}}{{\mathrm{\&Delta;x}}_i/2},$

在从上式中去掉Φ₀的值之后因此其变为

${\mathrm{Flux}}_{i,j,k}^{i+1,j,k}=\frac{2(\mathrm{\&Phi;}\left(r_{i+1,j,k}\right)-\mathrm{\&Phi;}\left(r_{i,j,k}\right))}{\frac{{\mathrm{\&Delta;x}}_i}{a_{i,j,k}{k}_{x,i,j,k}{S}_{x,i}}+\frac{{\mathrm{\&Delta;x}}_{i+1}}{a_{i+1,j,k}{k}_{x,i+1,j,k}{S}_{x,i+1}}}.---\left(2.22\right)$

由于S_x，i＝Δy_jΔz_i，j，k是当前计算单元的面的面积，因此使用(Δx_i)²/V_i，j，k替代Δx_i/S_x，i表达式是可能的，其中V_i，j，k＝Δx_iΔy_jΔz_i，j，k是单元的体积。

利用移动性项的标准迎风(upwinding)技术

则可将(2.22)重新写成以下形式

${\mathrm{Flux}}_{i,j,k}^{i+1,j,k}=\frac{{2a}_{i+1,j,k}^{\mathrm{upw}}(\mathrm{\&Phi;}\left(r_{i+1,j,k}\right)-\mathrm{\&Phi;}\left(r_{i,j,k}\right))}{\frac{{\left({\mathrm{\&Delta;x}}_i\right)}^2}{k_{x,i,j,k}{V}_{i,j,k}}+\frac{{\left({\mathrm{\&Delta;x}}_{i+1}\right)}^2}{k_{x,i+1,j,k}{V}_{i+1,j,k}}}.---\left(2.23\right)$

通过计算单元C_i，j，k的所有其他面的通量以相同的方式获得。C_i，j，k的面的传输率系数表示为

其中(α，β，γ)的集合包括{(i±1，j，k)，(i，j±1，k)，(i，j，k±1)}，因为

${\mathrm{Tr}}_{i,j,k}^{i-1,j,k}=\frac{{2a}_{i-1,j,k}^{\mathrm{upw}}}{\frac{{\left({\mathrm{\&Delta;x}}_i\right)}^2}{k_{x,i,j,k}{V}_{i,j,k}}+\frac{{\left({\mathrm{\&Delta;x}}_{i-1}\right)}^2}{k_{x,i-1,j,k}{V}_{i-1,j,k}}},$ ${\mathrm{Tr}}_{i,j,k}^{i+1,j,k}=\frac{{2a}_{i+1,j,k}^{\mathrm{upw}}}{\frac{{\left({\mathrm{\&Delta;x}}_i\right)}^2}{k_{x,i,j,k}{V}_{i,j,k}}+\frac{{\left({\mathrm{\&Delta;x}}_{i+1}\right)}^2}{k_{x,i+1,j,k}{V}_{i+1,j,k}}},$

${\mathrm{Tr}}_{i,j,k}^{i,j-1,k}=\frac{{2a}_{i,j-1,k}^{\mathrm{upw}}}{\frac{{\left({\mathrm{\&Delta;y}}_j\right)}^2}{k_{y,i,j,k}{V}_{i,j,k}}+\frac{{\left({\mathrm{\&Delta;y}}_{j-1}\right)}^2}{k_{y,i,j-1,k}{V}_{i,j-1,k}}},$ ${\mathrm{Tr}}_{i,j,k}^{i,j+1,k}=\frac{{2a}_{i,j+1,k}^{\mathrm{upw}}}{\frac{{\left({\mathrm{\&Delta;y}}_j\right)}^2}{k_{y,i,j,k}{V}_{i,j,k}}+\frac{{\left({\mathrm{\&Delta;y}}_{j+1}\right)}^2}{k_{y,i,j+1,k}{V}_{i,j+1,k}}},$

${\mathrm{Tr}}_{i,j,k}^{i,j,k-1}=\frac{{2a}_{i,j,k-1}^{\mathrm{upw}}}{\frac{{\left({\mathrm{\&Delta;z}}_{i,j,k}\right)}^2}{k_{z,i,j,k}{V}_{i,j,k}}+\frac{{\left({\mathrm{\&Delta;z}}_{i,j,k-1}\right)}^2}{k_{z,i,j,k-1}{V}_{i,j,k-1}}},$ ${\mathrm{Tr}}_{i,j,k}^{i,j,k+1}=\frac{{2a}_{i,j,k+1}^{\mathrm{upw}}}{\frac{{\left({\mathrm{\&Delta;z}}_{i,j,k}\right)}^2}{k_{z,i,j,k}{V}_{i,j,k}}+\frac{{\left({\mathrm{\&Delta;z}}_{i,j,k+1}\right)}^2}{k_{z,i,j,k+1}{V}_{i,j,k+1}}}.$

之后，在完全隐式的方案中方程(2.23)的通量可以在当前时间级近似为

${\left({\mathrm{Flux}}_{i,j,k}^{\mathrm{\&alpha;},\mathrm{\&beta;},\mathrm{\&gamma;}}\right)}^{\left(n\right)}={\left({\mathrm{Tr}}_{i,j,k}^{\mathrm{\&alpha;},\mathrm{\&beta;},\mathrm{\&gamma;}}\right)}^{\left(n\right)}({\mathrm{\&Phi;}}_{\mathrm{\&alpha;},\mathrm{\&beta;},\mathrm{\&gamma;}}^{\left(n\right)}-{\mathrm{\&Phi;}}_{i,j,k}^{\left(n\right)}).---\left(2.24\right)$

请注意存在不同类型的边界条件，其可被强加于给定层的面的部分上，例如封闭的边界、指定的流入物或流出物(诺埃曼类型)和指定压力(狄利克雷类型)。以下段落将描述这些边界条件。在说明书中，假设计算块C_i，j，k的一个面属于域的边界。为了使符号更为统一，该面表示为

并且该面上的势能增长将表示为

请注意对于内部面

${\mathrm{\&Delta;\&Phi;}}_{i,j,k}^{\mathrm{\&alpha;},\mathrm{\&beta;},\mathrm{\&gamma;}}={\mathrm{\&Phi;}}_{\mathrm{\&alpha;},\mathrm{\&beta;},\mathrm{\&gamma;}}-{\mathrm{\&Phi;}}_{i,j,k}.$

对于封闭的边界，不存在通量，因此

${\left({\mathrm{Flux}}_{i,j,k}^{\mathrm{\&alpha;},\mathrm{\&beta;},\mathrm{\&gamma;}}\right)}^{(*)}=0,$

其中，(^*)表示强加该条件时的时间级。

对于指定的流入物或流出物条件，本发明的流入物边界条件实施例假设流体流入域。因此，表达式

$Q_{i,j,k}^{\mathrm{\&alpha;},\mathrm{\&beta;},\mathrm{\&gamma;}}\&equiv;{\left({\mathrm{Flux}}_{i,j,k}^{\mathrm{\&alpha;},\mathrm{\&beta;},\mathrm{\&gamma;}}\right)}^{(*)}={\left(\underset{\&PartialD;C_{i,j,k}^{\mathrm{\&alpha;},\mathrm{\&beta;},\mathrm{\&gamma;}}}{\&Integral;}(n_s,\mathrm{aK}\&dtri;\mathrm{\&Phi;})\right)}^{(*)}$

应当是负的，因为n_s是外法向向量并且▽Φ方向向内。因此，对于该类型的边界，其设为

${\left({\mathrm{Flux}}_{i,j,k}^{\mathrm{\&alpha;},\mathrm{\&beta;},\mathrm{\&gamma;}}\right)}^{(*)}=Q_{i,j,k}^{\mathrm{\&alpha;},\mathrm{\&beta;},\mathrm{\&gamma;}}.$

其中，

否则，对于流出边界条件，应当指定流体流出物的正值

对于指定的压力边界条件，本发明的实施例假设在域边界的毛细管压力是小的并且层的倾斜是可忽略的。因此，边界面可以被视为与轴线d垂直(其中d可以是x、y或z)。当在边界面提供压力时(在面的中间点r_b)，方程组(2.20)的第一个方程可以如下修改为

${\mathrm{\&Phi;}}^{(i,j,k)}\left(r_b\right)-\mathrm{\&Phi;}\left(r_{i,j,k}\right)=\frac{1}{a\left(r_{i,j,k}\right)k_d\left(r_{i,j,k}\right)}{(w,n_d)}_{i,j,k}\frac{{\mathrm{\&Delta;d}}_{i,j,k}}{2},$

其中Δd_i，j，k是Δx_i或Δy_j或Δz_i，j，k，这取决于面。因此，对应的通量由下式确定

${\mathrm{Flux}}_{i,j,k}^{\mathrm{\&alpha;},\mathrm{\&beta;},\mathrm{\&gamma;}}=(\mathrm{\&Phi;}\left(r_b\right)-\mathrm{\&Phi;}\left(r_{i,j,k}\right))\frac{a_{i,j,k}{k}_{d,i,j,k}{V}_{i,j,k}}{{\left({\mathrm{\&Delta;d}}_{i,j,k}\right)}^2/2},$

并且传输率由下式确定

${\mathrm{Tr}}_{i,j,k}^{\mathrm{\&alpha;},\mathrm{\&beta;},\mathrm{\&gamma;}}=\frac{{2a}_{i,j,k}{k}_{d,i,j,k}{V}_{i,j,k}}{{\left(\mathrm{\&Delta;}{d}_{i,j,k}\right)}^2}.$

对于垂直于x或y轴的边界面，边界项为

$\left({\mathrm{Flux}}_{i,j,k}^{i\&PlusMinus;1,j,k}\right)=\left({\mathrm{Tr}}_{i,j,k}^{i\&PlusMinus;1,j,k}\right)({\mathrm{\&Phi;}}_b-{\mathrm{\&Phi;}}_{i,j,k}),{\mathrm{Tr}}_{i,j,k}^{i\&PlusMinus;1,j,k}=\frac{{2a}_{i,j,k}{k}_{x,i,j,k}{V}_{i,j,k}}{{\left({\mathrm{\&Delta;x}}_i\right)}^2},---\left(2.25\right)$

或

$\left({\mathrm{Flux}}_{i,j,k}^{i,j\&PlusMinus;1,k}\right)=\left({\mathrm{Tr}}_{i,j,k}^{i,j\&PlusMinus;1,k}\right)({\mathrm{\&Phi;}}_b-{\mathrm{\&Phi;}}_{i,j,k}),{\mathrm{Tr}}_{i,j,k}^{i,j\&PlusMinus;1,k}=\frac{{2a}_{i,j,k}{k}_{y,i,j,k}{V}_{i,j,k}}{{\left({\mathrm{\&Delta;y}}_j\right)}^2}.---\left(2.26\right)$

对于垂直于z轴的边界面(z_b-z_i，j，k)＝±Δz_i，j，k/2，其中对于上面的面(z_b＜z_i，j，k)，符号为“-”，对于底面(z_b＞z_i，j，k)，符号为“+”。通常，在盆地建模仿真中，计算域的底部不存在流边界条件并且域顶部的压力被指定。对于顶面，通量由下式给出

$\left({\mathrm{Flux}}_{i,j,k}^{i,j,k-1}\right)=\left({\mathrm{Tr}}_{i,j,k}^{i,j,k-1}\right)({\mathrm{\&Phi;}}_b-{\mathrm{\&Phi;}}_{i,j,k}),{\mathrm{Tr}}_{i,j,k}^{i,j,k-1}=\frac{{2a}_{i,j,k}{k}_{y,i,j,k}{V}_{i,j,k}}{{\left(\mathrm{\&Delta;}{z}_k\right)}^2}.---\left(2.27\right)$

可以通过将项

添加到右手侧的向量中并将所有其余的添加到矩阵的对角项的方式将表达式(2.25)、(2.26)和(2.27)合并到矩阵方程(2.13)中。之后，通过将解Φ_i，j，k替换回方程(2.25)、(2.26)和(2.27)中获得速率。

对于方程(2.10)中单元厚度的计算，超压的近似值在单元的顶部边界和底部边界是有用的，也就是和

因此，由于最顶部单元的顶部边界处的压力边界条件，应当强加以下条件

${\mathrm{\&Phi;}}_{i,j,1}^{\mathrm{top}}=0.$

在域的底部，不存在流动边界条件，由于此原因，对于最低的单元的底部边界可能具有以下压力条件

${\mathrm{\&Phi;}}_{i,j,n_z}^{\mathrm{bot}}={\mathrm{\&Phi;}}_{i,j,n_z}.$

对于任意其他边界，总是必须存在顶部和底部相邻单元，而且由于超压连贯性，

${\mathrm{\&Phi;}}_{i,j,1}^{\mathrm{bot}}={\mathrm{\&Phi;}}_{i,j,k+1}^{\mathrm{top}},k=1,...,n_z-1.---\left(2.28\right)$

两个相邻单元之间的边界处超压的值利用上述段落中导出的通量连续性条件，也就是方程(2.21)限定，垂直于z方向的面的超压可以表示成

${\mathrm{\&Phi;}}_{i,j,1}^{\mathrm{bot}}=\frac{{\mathrm{\&Phi;}}_{i,j,k}{k}_{z,i,j,k}/{\mathrm{\&Delta;z}}_{i,j,k}+{\mathrm{\&Phi;}}_{i,j,k+1}{k}_{z,i,j,k+1}/{\mathrm{\&Delta;z}}_{i,j,k+1}}{k_{z,i,j,k}/{\mathrm{\&Delta;z}}_{i,j,k}+k_{z,i,j,k+1}/{\mathrm{\&Delta;z}}_{i,j,k+1}}.---\left(2.29\right)$

非线性方程系统

从方程组(2.15)、(2.23)和(2.24)，方程(2.13)的离散版本包含计算单元的未知厚度Δz和超压Φ的值，以及函数k_x、k_y、k_z和ρ_a，它们又取决于平均单元孔隙度φ、流体静压p_h和超压Φ。厚度Δz的值可以根据方程(2.10)确定，该方程包含未知量Δz和Φ以及静岩负荷L、流体密度ρ_a的值和函数k_x、k_y、k_z，考虑孔隙度和厚度之间的方程(2.4)，渗透率系数k_x、k_y、k_z可以重写成Δz的函数。为了闭合用于确定Δz和Φ的方程系统，需要L和p_h的两个额外方程，也就是方程(2.12)和(2.8)。因此，描述压实介质中流体流动的该组未知量包含四个变量，即超压Φ、单元厚度Δz、静岩负荷L和流体静压p_h。

引入包含如下四个变量的未知量的向量

X＝(Φ，Δz，L，P_h)，

其具有以下子向量

Φ＝{Φ_i，j，k}，Δz＝{Δz_i，j，k}，L＝{L_i，j，k}，P_h＝{p_{h；i，j，k}}，

其中Φ_i，j，k是超压，Δz_i，j，k是单元厚度，L_i，j，k是静岩负荷，并且p_{h；i，j，k}是流体静压。

之后，方程组(2.1)的离散化可以以非线性代数方程系统的形式书写为

F(X(n))＝0，

或以分量方式的形式(对于内部单元(i，j，k))书写为

$F_{\mathrm{\&Phi;};i,j,k}\&equiv;{\mathrm{\&Delta;x}}_i{\mathrm{\&Delta;y}}_j{\mathrm{\&rho;}}_{a;i,j,k}^n({\mathrm{\&Delta;z}}_{i,j,k}^n-s_{i,j,k}^n)-{\tilde{M}}_{i,j,k}^{n-1}+---\left(3.1\right)$

$\mathrm{\&Delta;t}\{{\mathrm{Tr}}_{i,j,k}^{i,j,k-1}({\mathrm{\&Phi;}}_{i,j,k}^n-{\mathrm{\&Phi;}}_{i,j,k-1}^n)+{\mathrm{Tr}}_{i,j,k}^{i,j,k+1}({\mathrm{\&Phi;}}_{i,j,k}^n-{\mathrm{\&Phi;}}_{i,j,k+1}^n)+$

${\mathrm{Tr}}_{i,j,k}^{i-1,j,k}({\mathrm{\&Phi;}}_{i,j,k}^n-{\mathrm{\&Phi;}}_{i-1,j,k}^n)+{\mathrm{Tr}}_{i,j,k}^{i+1;j,k}({\mathrm{\&Phi;}}_{i,j,k}^n-{\mathrm{\&Phi;}}_{i+1,j,k}^n)+$

${\mathrm{Tr}}_{i,j,k}^{i,j-1,k}({\mathrm{\&Phi;}}_{i,j,k}^n-{\mathrm{\&Phi;}}_{i,j-1,k}^n)+{\mathrm{Tr}}_{i,j,k}^{i,j+1,k}({\mathrm{\&Phi;}}_{i,j,k}^n-{\mathrm{\&Phi;}}_{i,j+1,k}^n)=0,$

$F_{L;i,j,k}\&equiv;L_{i,j,k}^n-L_{i,j,k-1}^n-\frac{g}{2}(({\mathrm{\&rho;}}_{s;i,j,k}-{\mathrm{\&rho;}}_{a;i,j,k}^n)s_{i,j,k}^n+---\left(3.3\right)$

$({\mathrm{\&rho;}}_{s;i,j,k-1}-{\mathrm{\&rho;}}_{a;i,j,k-1}^n)s_{i,j,k-1}^n)=0,$

$F_{p_h;i,j,k}\&equiv;p_{h;i,j,k}^n-p_{h;i,j,k-1}^n-\frac{g}{2}({\mathrm{\&rho;}}_{a;i,j,k}^n{\mathrm{\&Delta;z}}_{i,j,k}^n+{\mathrm{\&rho;}}_{a;i,j,k-1}^n{\mathrm{\&Delta;z}}_{i,j,k-1}^n)=0.---\left(3.4\right)$

第一组方程(3.1)中的项是两个表达式(2.14)与(2.16)之和

${\tilde{M}}_{i,j,k}^{n-1}={\mathrm{\&Delta;x}}_i{\mathrm{\&Delta;y}}_j({\mathrm{\&Delta;z}}_{i,j,k}^{n-1}{\mathrm{\&rho;}}_{a;i,j,k}^{n-1}{\mathrm{\&phi;}}_{i,j,k}^{n-1}+{\mathrm{\&rho;}}_{a;i,j,k}^{*}{\mathrm{\&phi;}}_{i,j,k}^{*}\frac{q_s\left(t_{n-1}\right)}{{\mathrm{\&rho;}}_{s;i,j,k}(1-{\mathrm{\&phi;}}_{i,j,k}^{*})}),$

其中标号^*表示该值分别在表面(输入数据)或从之前的用于沉积或侵蚀的时步t_n-1获得。流体密度由分别用于不可压缩或可压缩的流体流动的方程(2.6)或方程(2.7)限定。所述方程限定过压、单元厚度和单元的沉积负荷。这三个方程可以用于定义包括不可压缩流体的域。如果流体可压缩，则需要流体静压的方程描述域。

传输率

由具有针对边界单元的修改的(2.24)来限定，如边界条件部分描述的。

第二组方程(3.2)中的和

的值由表达式(2.11)定义，而和

的值利用表达式(2.28)和(2.29)计算。方程(3.3)和(3.4)通过以下方式在顶部边界被扩充

$F_{L;i,j,1}\&equiv;L_{i,j,1}^n-\frac{g}{2}\&CenterDot;({\mathrm{\&rho;}}_{s;i,j,1}-{\mathrm{\&rho;}}_{a;i,j,1}^n)s_{i,j,1}^n=0,$

$F_{p_h;i,j,1}\&equiv;p_{h;i,j,1}^n-p_{h;i,j,\mathrm{surf}}^n-\frac{g}{2}\&CenterDot;{\mathrm{\&rho;}}_{a;i,j,1}^n\mathrm{\&Delta;}{z}_{i,j,1}^n=0.$

本发明的实施例使用相容的方程组同时描述域的压实和域的流体流动。本发明的实施例平衡质量、动量和本构关系以确定压实和/或松散的域。本发明的实施例描述域中的流体流动。本发明的实施例引入未知量描述孔隙度。孔隙度可以依赖于有效应力，有效应力为物理行为，依赖于压力和负荷并且是压实的结果。

请注意本发明的其他实施例可能涉及其他未知变量。例如，本发明的另一实施例可能利用总压力、流体静压、厚度和有效应力描述域的流体流动和压实。可以使用任何未知量组，只要该组未知量在域上是有解的(相容的)。可以将其他变量添加到该组方程，例如，温度和描述其空间和时间分布的附加方程或多个方程。通常，方程系统(3.1)-(3.4)涉及的系数并不强依赖于其他变量，例如温度，因此为了使描述简单，未考虑这些附加的变量。

以上概述的各种过程和方法可以根据本发明的各实施例组合在一个或更多不同方法中、用在一个或更多不同系统中或用在一个或更多不同的计算机程序产品中。

例如，一个示例性方法500可以在计算机上建模物理区域，如图5所示。如上所述，物理区域可以通过发生在区域内的一种或更多过程或现象来建模，块501。例如，在地下地质盆地中，流体流动和沉积物压实可以用于建模盆地。因此，通过建模沉积物的累积和/或侵蚀，以及流体如何使沉积物流动，可以获得盆地的精确模型。对这种现象建模可能是困难的，因为流体流动和压实是相关的，其中流体流动依赖于压实，反之亦然。

根据本发明的实施例，模型使用一组方程来描述现象，块502。例如，如果流体是不可压缩的，例如水或油，表示区域的过压(或超压)、区域的厚度和沉积物负荷的一组方程可用于描述关联的流体流动和压实现象。如果流体是可压缩的，例如气体或天然气，则可以使用表示流体静压的附加方程。

可以通过在模型上施加一个或更多假设来简化方程，块503。尽管在对比模型和实际的物理盆地时假设可能引入误差或不准确性，但是假设允许以有效的计算方式求解方程组。这些假设可以施加到现象或方程组上。例如，一个假设可以是沉积物累积速率是已知的。物理盆地中的实际速率可能不是已知的，因此可以为模型假设速率。另一假设可以是压实只发生在竖直方向。换句话说，在横向方向上不发生压实。另一假设可能是压实是相对不可逆转的。这意味着在沉积物的侵蚀过程或抬升过程中沉积物将只是大部分压实，同时伴有一定数量的松散发生，但是并非完全返回到初始状态。本发明的实施例可以使用其他假设。

在简化了方程组后，可以对其求解以利用数据同时描述这两种现象，块504。通过求解方程组，模型将准确描述区域中现象的发生。之后，模型可用于帮助对物理区域进行更改。例如，模型可以用于从盆地更为高效地开采地下石油或气体。

请注意本文描述的任意功能可以在硬件、软件和/或固件、和/或其任意组合中实施。当在软件中实施时，本发明的要素实质是执行必要任务的代码段。程序或代码段可以存储在计算机可读介质中或通过计算机数据信号来传递。“计算机可读介质”可包括可存储或传递信息的任意介质。计算机可读介质的示例包括电路、半导体存储装置、ROM(只读存储器)、闪存、可擦写ROM(EROM)、软盘、压缩盘CD-ROM、光盘、硬盘、光纤介质等。计算机数据信号可包括可在传输介质诸如电子网络通道、光纤、空气、电磁、RF链路等上传播的任意信号。代码段可以通过计算机网络诸如互联网、企业内部网等下载。

图6示出了适于使用本发明的计算机系统600。中央处理单元(CPU)601连接到系统总线602。CPU 601可以是任意通用CPU，诸如英特尔奔腾处理器。但是，只要CPU 601支持本文描述的发明的操作，本发明并不限于CPU 601的架构。总线602连接到随机存取存储器(RAM)603，随机存取存储器可以是SRAM、DRAM或SDRAM。ROM 604也连接到总线602，ROM 604可以是PROM、EPROM或EEPROM。RAM 603和ROM 604保存用户和系统数据和程序，这是本领域公知的。

总线602也连接到输入/输出(I/O)控制器卡605、通信适配器卡611、用户接口卡608和显示卡609。I/O适配器卡605将存储装置606(诸如硬盘驱动器、CD光盘驱动器、软盘驱动器、磁带驱动器中的一个或更多)连接到计算机系统。I/O适配器605可以连接到印刷机，印刷机将允许系统打印信息的纸件，诸如文档、图片、文章等。注意，印刷机可以是打印机(例如喷墨、激光打印机等等)、传真机或复印机。通信卡611适于将计算机系统600连接到网络612，网络612可以是电话网络、局域网(LAN)和/或广域网(WAN)、以太网和/或因特网网络中的一个或更多。用户接口卡608将用户输入装置诸如键盘613和指点装置607连接到计算机系统600。用户接口卡608也可以通过扬声器(一个或多个)提供声音输出至用户。显示卡609由CPU 601驱动以控制在显示装置610上进行显示。

尽管已经详细描述了本发明及其优势，但是应当理解可以做出各种改变、替换和变更而不偏离本发明在所附权利要求书限定的精神和范围。而且，本申请的范围并不受此说明书描述的过程、机器、制造、物质构成、手段、方法和步骤的限制。本领域的普通技术人员将易于理解根据本发明的公开内容，可以利用根据本发明的现有或以后开发的过程、机器、制造、物质构成、手段、方法或步骤，其执行与本文描述的对应实施例基本相同的功能或达到基本相同的效果。因此，所附权利要求书意在将这些过程、机器、制造、物质构成、手段、方法或步骤包括在其范围内。

Claims (28)

1.一种建模物理区域的方法，其包括：

接收数据，所述数据定义所述物理区域的至少一种物理特性；

选择第一现象和第二现象，其中所述第一现象和第二现象在所述物理区域上关联以用于建模；

定义一组方程，该组方程描述所述第一现象和第二现象，其中所述方程在所述物理区域上是有解的；

通过在所述第一现象、所述第二现象和该组方程中的至少一个上施加至少一个假设来简化该组方程；和

求解该组方程以利用数据同时描述所述两种现象。

2.根据权利要求1所述的方法，其中所述物理区域是地下地质盆地并且所述两种现象是盆地中的流体流动和物质压实，其中所述流体位于所述盆地内。

3.根据权利要求2所述的方法，其中所述流体是石油、天然气、水、液体、气体和具有放射性同位素的流体中的至少一种。

4.根据权利要求2所述的方法，其中所述物质是沉积物。

5.根据权利要求4所述的方法，其中所述至少一个假设是以下中的至少一个：

沉积物累积速率是已知的；

所述压实只发生在垂直方向上；和

所述压实是相对不可逆转的。

6.根据权利要求2所述的方法，进一步包括：

在所述物理区域的模型上提供网格，其中所述网格包括多个单元。

7.根据权利要求6所述的方法，其中所述求解针对所述网格的每个单元执行。

8.根据权利要求6所述的方法，其中在建模过程中所述网格的每个单元在垂直方向发展以对随时间的物质累积建模。

9.根据权利要求8所述的方法，其中在建模过程中，至少一个单元掩埋在所述模型中，因为其他单元在所述一个单元上发展。

10.根据权利要求8所述的方法，其中每个单元是平行六面体单元。

11.根据权利要求8所述的方法，其中定义单元的水平面的x方向和y方向与地层时间线对齐。

12.根据权利要求6所述的方法，其中所述流体是可压缩流体，并且该组方程包括：

定义每个单元的过压的第一方程、定义每个单元的单元厚度的第二方程、定义每个单元的物质负荷的第三方程和定义每个单元的流体静压的第四方程。

13.根据权利要求6所述的方法，其中所述流体是不可压缩流体，并且该组方程包括：

定义每个单元的过压的第一方程、定义每个单元的单元厚度的第二方程和定义每个单元的物质负荷的第三方程。

14.根据权利要求6所述的方法，进一步包括：

向单元应用至少一种转换；

其中所述转换是沉积、下举、抬升和侵蚀之一。

15.根据权利要求6所述的方法，进一步包括：

在和所述区域的边缘相邻的单元上施加至少一个边界条件。

16.根据权利要求1所述的方法，其中所述物理区域是地下地质盆地，并且所述模型涉及地下石油，并且所述求解有助于从所述盆地开采石油。

17.根据权利要求1所述的方法，进一步包括：

从来自传感器的信息导出所述数据，所述传感器测量所述物理区域的所述至少一个物理特性。

18.一种具有计算机可读介质的计算机程序产品，该计算机可读介质具有记录在其上用于在计算机上建模地下地质盆地的计算机程序逻辑，其包括：

定义一组方程的代码，该组方程描述流体流动和沉积物压实，其中所述方程在所述盆地上是有解的，并且其中代码通过在所述流体流动、沉积物压实和该组方程中的至少一个上施加至少一种假设被简化；和

用于求解该组方程的代码以同时描述所述盆地中的所述流体流动和沉积物压实。

19.根据权利要求18所述的计算机程序产品，进一步包括：

在所述盆地的模型上提供网格的代码，其中所述网格包括多个单元。

20.根据权利要求19所述的计算机程序产品，其中该组方程包括：

描述第一方程的代码，所述第一方程定义每个单元的过压；

描述第二方程的代码，所述第二方程定义每个单元的单元厚度；和

描述第三方程的代码，所述第三方程定义每个单元的物质负荷。

21.根据权利要求19所述的计算机程序产品，进一步包括：

向单元应用至少一种转换的代码；

其中所述转换是沉积、下举、抬升和侵蚀之一。

22.根据权利要求19所述的计算机程序产品，其中所述至少一个假设包括：

沉积物累积速率是已知的第一假设；

所述压实只发生在垂直方向上的第二假设；和

所述压实是相对不可逆转的第三假设。

23.根据权利要求18所述的计算机程序产品，其中所述流体是石油。

24.一种在计算机上建模地下地质盆地的方法，包括：

接收数据，所述数据定义所述盆地的至少一种物理特性；

定义一组方程，该组方程描述所述盆地中的流体流动和沉积物的压实，其中所述方程在所述物理区域上是有解的；

通过施加所述压实只发生在垂直方向的假设来简化该组方程；和

求解该组方程以利用所述数据同时描述所述两种现象。

25.根据权利要求24所述的方法，其中所述模型涉及地下石油，所述方法进一步包括：

从来自传感器的信息导出所述数据，所述传感器测量所述物理区域的所述至少一个物理特性；和

使用已求解的方程帮助从所述盆地中开采石油。

26.根据权利要求25所述的方法，其中所述物理区域是地下地质盆地，并且所述两种现象是所述盆地中的流体流动和物质的压实，所述流体位于所述盆地中。

27.根据权利要求25所述的方法，进一步包括：

基于所述已求解的方程组产生所述地下地质盆地的盆地模型；

基于所述盆地模型预测所述物理区域内碳氢化合物的位置；和

布置生产设施以基于所预测的碳氢化合物的位置开采所述物理区域内的碳氢化合物。

28.根据权利要求27所述的方法，进一步包括基于所述盆地模型评估所述物理区域的生产潜能。

Images