Novel basin modelling concept for simulating deformation from mechanical compaction using level sets

As sedimentation progresses in the formation and evolution of a depositional geologic basin, the rock strata are subject to various stresses. With increasing lithostatic pressure, compressional forces act to compact the porous rock matrix, leading to overpressure buildup, changes in the fluid pore pressure and fluid flow. In the context of petroleum systems modelling, the present study concerns the geometry changes that a compacting basin experiences subject to deposition. The purpose is to track the positions of the rock layer interfaces as compaction occurs. To handle the challenge of potentially large geometry deformations, a new modelling concept is proposed that couples the pore pressure equation with a level set method to determine the movement of lithostratigraphic interfaces. The level set method propagates an interface according to a prescribed speed. The coupling term for the pore pressure and level-set equations consists of this speed function, which is dependent on the compaction law. The two primary features of this approach are the simplicity of the grid and the flexibility of the speed function. A first evaluation of the model concept is presented based on an implementation for one spatial dimension accounting for vertical effective stress. Isothermal conditions with a constant fluid density and viscosity were assumed. The accuracy of the implemented numerical solution for the case of a single stratigraphic unit with a linear compaction law was compared to the available analytical solution [38]. The multi-layer setup and the nonlinear case were tested for plausibility.


Introduction
Basin modelling has become an increasingly important component to the identification of hydrocarbons in various geological settings, see e.g., [15,26,34,35]. The thermal history and pore pressure of a basin as it evolves over millions of years is vital to estimating the generation, migration and trapping of gas and oil deposits.
The motivation of the present work is to focus on handling significant changes in the stratification of the sample system. Extensional and compressive forces can be large. For example, uplift and salt diapirs are potential causes of large deformations in the layering of the basin. These large deformations require specific computational techniques to render their inclusion in the model tractable. A well-known approach to computing the dynamical evolution of a system is the Lagrangian framework. In this setting, the coordinate system of the grid or mesh is updated to track the changes in the material as it moves. The grid nodes are identified with material positions. If the material moves significantly, then the grid nodes will as well. This can lead to warps and distorted grid geometry, which degrade the numerical properties of the discrete solution to the governing differential equations. To correct this, one needs to re-grid, which, in general, is an expensive procedure. An alternative to the re-gridding process is to hold the coordinate system of the grid fixed, while moving the material through. Known as the Eulerian framework, this fixed grid approach has the numerical advantages of a well-behaved, regular grid, but has the complication that a material element will, in general, not be identified with a single node. Locating the material elements on the grid leads to its own challenges. For a general review of the two computational approaches, see [5]. A popular method to address large deformations, the Arbitrary Lagrangian-Eulerian technique, employs a combination in order to seek the best of both worlds, see e.g., [11,13,23].
In contrast, the level set method [25,28] exploits the fixed grid, while introducing a potential function to implicitly locate the interfaces between material elements. Level set methods have been applied to a wide array of domains in the last decades, from computational fluid dynamics (e.g., [31,37]), groundwater flow [12], crystal growth [30] to image processing [20]. This study uses the level set method to investigate the issue of deformation of the geometry of a sedimentary basin due to mechanical compaction.
As basin modelling is a multi-component and multi-scale approach to complex physico-chemical processes, the focus of the current study is restricted to the problem of fluid flow and compaction that a basin undergoes as layers of sediment contribute to an evolving pressure profile of the basin. The compaction, i.e., reduction in porosity, contributes to changes in pore pressure, thus driving flow, as a result of increasing sediment load and displacement of strata. Furthermore, in this paper, we consider only a one-dimensional model of single phase fluid flow under isothermal conditions. The present work draws on the following studies in compaction and fluid flow.
Schrefler and Scotta [27] present a fully coupled model of multiphase flow in a deformable poroelastic medium. They use a finite element method to solve for the displacements of the solid matrix and the pressures of the two phases. The other variables are related through constitutive equations. A Galerkin method is used for the spatial discretization and a Newmark method is used for the temporal discretization. They validate the numerical results with three physical examples over small spatial scales and short time scales. In Minkoff et al. [22], coupling is studied through two 3-D finite element codes. With the mixed finite element flow simulator, a black oil model is run with initial values for porosity and permeability to compute a pore pressure. This is, then, communicated to the quasi-static geomechanics simulation, which computes porosity and permeability. These results are passed to the flow simulator and the subsequent time step flow is calculated, and so on. The authors argue for a loose coupling approach essentially since it is computationally less expensive than a fully coupled method. One advantage is that the domains can be quite different for each system. In Bethke [6], the author provides a numerical model of compaction-driven flow for a 2-dimensional space. The conservation of rock mass, fluid flow and heat transfer each contribute an equation to the coupled system, which is solved iteratively using a finite difference discretization. The problem is formulated using a Lagrangian reference frame that moves with the rock mass. Material added to the top of the basin will form new nodes when sufficient material is added. See also [7] and [3]. Chen et al. [9] presents a similar model for compacting soils, but differs in that only compaction and fluid flow are considered. The triangulation of the basin is updated at each time step, i.e., the elements are functions of time. The coupled system is iteratively solved using the mixed finite element method. Kikinzon et al. [17] follows the model set forth in [9], while modifying the algorithm for the timedependent grid structure. They clearly divide the problem into physically distinct steps which are treated in sequence, an operator splitting approach. The presentation of the Theory section below draws on this exposition. Christopher [10] studies the Ursa region of the Gulf of Mexico by using a 1-D compaction model with finite difference solution method for the fluid pressure. In the present study, as in the above series of papers, an effective stress principle with only vertical compaction is assumed. The same basic model of compaction-driven fluid flow is considered here. Also, the sedimentation step is split from the fluid pressure solution step, pursuing an iterative approach to the coupled partial differential equations.
In contrast to the above series of papers, however, the level set method is introduced in this paper to implicitly track the evolving locations of the rock layers in the sedimentary basin. Work of Longoni et al. [18,19] has already used implicit interface tracking in the context of sedimentary basin modelling. This approach does not explicitly model mechanical compaction, as each layer is thought of as an immiscible non-Newtonian fluid. Solved with an adaptive finite element method, the Stokes equations form the core of their model, with a level set advection framework to keep the fluid layers separate and an arbitrary Lagrangian-Eulerian framework to represent the boundary motion of the basin under deposition. To the best of the authors' knowledge, the current work is the first application of the level set method to the problem of modelling a sedimentary basin as a porous media fluid-structure interaction problem.
In our application, we use the level set method in the following way. Adopting an operator splitting approach, first the layer interfaces are propagated in response to sedimentation. The state and material properties at the nodes are changed to reflect the adjusted interfaces. In other words, the rock properties are transported through the domain according to the level set interfaces, rather than the grid nodes being moved around. Then, the fluid pressure equation is solved and a compaction rule is applied to determine the contribution of compaction to the movement of the interface in the next time step.
Thus, the rock is effectively advected through the domain via the level set interfaces, with the goal of avoiding the complex and expensive update rules for moving grids. This is the novelty of the approach considered here. It has two main upshots: (1) simplicity of the grid and (2) flexibility of the speed function. The grid remains regular, which is advantageous from an implementation point of view. Secondly, the setup of the coupling between physics and geometry, allows for free choice of the physical models that move the interfaces. What is required is to have a speed function which determines the movement of the front at each point on the interface.

Effective stress models
The pores, or void spaces, in a solid material allow for fluid, liquid or gas, to occupy the space. The fraction of the pore volume to the total volume is named the porosity. As stress exerted on the material affects its porosity, this in turn has consequences for the fluid flow through the porous medium. A central tool to modelling flow in porous material is the concept of effective stress introduced in [33]. The effective stress, σ [ML −1 T −2 ], is given by the total stress, σ [ML −1 T −2 ], minus the pore pressure, This approach assumes that the stress on the solid matrix is a result of solid particles contacting other particles. The pressure of the fluid in the solid matrix exerts a force on the constituent solids that is homogeneous and so does not contribute to deformation of the solid matrix. For this reason, the fluid pressure is subtracted from the total stress to isolate the effective stress that the solid matrix experiences. In the context of the one-dimensional model described in this paper, we consider effective stress in only one direction, a scalar, which we refer to as the vertical effective stress.
The importance of the effective stress hypothesis is that it directly relates the fluid pressure, essential for describing the flow dynamics, to the stress which is exerted on the solid skeleton of the medium, the mechanical aspect. For a historical perspective on the development of porous media mechanics until recently, see [8].

Consolidation
Consider a consolidation problem, where sediment is added to the top of a soil column, as in [14]. The sediment stack deforms continuously under its own growing weight. The fluid pressure, above the hydrostatic pressure, in the saturated pore space holds open the pores, resisting compaction. Extending this to different types of sediment, i.e., with different lithological properties, that are added at different stages, creates the stratification that we conceptualize as rock layer interfaces.
Following [15], the Darcy flow relation between the discharge velocity of the fluid (relative to the solid) and the overpressure gradient, including gravity, is assumed, where v [LT −1 ] is the discharge velocity of the pore fluid, is the density of the fluid, z [L] is the spatial position, and g [LT −2 ] is the acceleration due to gravity. The mass balance condition requires that any divergence be compensated by a reduction in the amount of mass contained, with φ [L 3 L −3 ] the porosity. The necessary reduction in mass comes from a change in the porosity, i.e., a reduction in volume for the fluid, or a change in fluid density. Neglecting changes in the fluid density, the last term in Eq. 3 is dropped, and inserting Eq. 2 into Eq. 3 gives, The change in porosity can be expressed as where C [M −1 LT 2 ] is the compressibility, a lithological property dependent on σ z , which is the effective stress in the z direction, as introduced previously. Recalling that effective stress is the total stress minus the pore pressure, σ z = σ z − p, Eq. 5 is combined with Eq. 4 to arrive at the pore pressure equation, Recall that for the purposes of this paper we consider only a one-dimensional model.
Introducing the overpressure, defined as pore pressure minus the hydrostatic pressure, p o = p −ρ f gz, for constant fluid density, and rearranging, we have The total overburden at a depth z is given by the weight of the material above, where the bulk density, ρ b , is the porosity-dependent average density of the material, i.e., with ρ s the density of the rock matrix and ρ f the fluid density. The right hand side of Eq. 8 becomes, This is the source term for the overpressure generation, the rate of change of the difference between the overburden and the hydrostatic pressure. The equation includes three lithological properties: the compressibility, the permeability and the rock density, which can vary vertically.
For mathematical closure, boundary and initial conditions must be specified. The top boundary condition is chosen to be a fixed atmospheric pressure, p 0 , for all t > 0. The bottom boundary condition is chosen to be a no flow condition, for all t > 0. For the initial condition, t = 0, In addition to the pore pressure equation, it is necessary to estimate the porosity as a function of depth or effective stress. A number of phenomenological models are presented in [15]. A good selection is the simple exponential relation known as Athy's Law [2], where φ 0 is an initial porosity and β is a compaction parameter. Bahr et al. [4] provide an analysis of the suitability of assuming this exponential relation. The main variables for the physical model considered here are summarized in Table 1.

Level set methods
When using level sets, the changing geometry is represented by evolving interfaces which are implicitly defined over a fixed grid [29]. Again, the motivating idea is that a fixed grid has the virtue of computational simplicity relative to frequently remapping the grid. Of course, there is a price to be paid, namely the introduction of additional partial differential equations to solve. As mentioned above, the key idea is to implicitly track the evolution of interfaces through the use of an additional equation, namely that for a higher dimensional function, the level set potential function, [24]. The curve that is to be tracked, Γ , is embedded in the level set potential function as the 0 level set, meaning all the points where the potential function is 0 represent Γ . The level set function, u, is defined such that at t = 0, its value at x is the signed distance, d, from x to Γ , i.e., The time evolution of u is given by an advection equation for the level set function, ∂u ∂t where w is the velocity field for the interface motion and u is the level set function . The interface is moved along in the direction of the velocity field. So, at any time going forward, Γ remains identified by all the points where the level set function has value 0. Defining the speed function, F , to be the the magnitude of the vector field normal to the front surface, F = w ·n, and with the definition of the unit normal vector,n = ∇u |∇u| , Eq. 17 can be rewritten as The velocity field for the interface motion, w, or alternatively the speed function, F , plays a determinative role in this scheme. For any particular physical application, F must come from the dynamics of the system. For applications that solve the Navier-Stokes equation, the speed function can be taken to be the fluid velocity, [39]. Vila et al. [36] use the level set method to track the evolution of resin being infused into a porous material. This problem is formulated on the assumption of Darcy flow, and the average fluid velocity is used as the speed function for the level set evolution.
The movement of the level set is determined by the relevant physics of the problem considered. This is a major feature of the method in that F can be constructed to take into account various physical phenomena at different levels of description. In the current section, we propose a speed function to capture the dynamics of the consolidation problem considered. The interfaces to be tracked are defined by the different rock types. The speed at which these interfaces move is mostly determined by the sedimentation rate, i.e., the rate at which material is added to the column. Further, there should be a contribution from the compaction of the layers. The following procedure constructs a speed function that satisfies this dependence on sedimentation rate and rate of compaction of the layers. It is not offered as a rigorous derivation, but rather as a concrete formula for the coupling.
Beginning with Equation B.2 in [15], we have from the conservation of rock mass, where v s is the velocity of the rock relative to the fixed coordinate system. Integrating in one dimension from the top of the column to a point z in depth, On the left hand side, we have the velocity at the depth z minus the velocity at the surface. With a constant time step and grid spacing, substituting the discrete differences for derivatives, the integration is then approximated as summation over the quantities evaluated on the discrete grid points, leading to The rock velocity is then used in the level set equation as the speed function. For numerical reasons, the level set function deviates from being a signed distance function as the iterations of the advection equation increase [21,24]. Therefore, [32] proposed a procedure to reinitialize the level set function back to being a signed distance function. In general, there will be some error introduced into the position of the interface from solving the reinitialization equation. For example, in a two-phase fluid problem with different fluid densities, this affects the mass conservation. Therefore, the recommended frequency with which this procedure should be applied depends on the application [16]. However, in our one-dimensional model, reinitialization does not affect the zero level set, but only the values of the level set function off the interface (since we always know the distance to the interface on the line). When the slope around the interface becomes too flat or steep, depending on some tolerance, the level set is reinitialized.
Another technical aspect of the level set method that is used in our model is extending the interface speed [1,24]. A speed must be assigned to regions neighboring the interface, as well as directly on it. In our case, the speed function computed within the layer above the bounding interface is extended into the layer below. This is required for smoothness reasons, since we have rock layers with heterogeneous properties, and is implemented by taking the slope in the speed function above the interface and extending the speed function linearly for a few nodes (a heuristic number has been the ratio of rock length added in a time step to the dicretization step). This amounts to something similar to upwinding of the speed of the rock layer above into the top of the layer below.

One dimensional example
As a demonstration, for the one dimensional case of an interface moving with speed F , one time step in the propagation of the front labelled Γ 1 , is illustated in Fig. 1. The interface between the blue section and the brown section, Γ 1 , starts at a position, z 0 , and moves to z 1 in the next time step. z 1 is computed through solving the advection equation for the level set function, u. The level set function, u, is plotted beneath a graph of the real coordinates. Given an intial level set function, u t 0 = z − z 0 , Eq. 18, can be used to find the level set function at the next time step, u t 1 . Using a simple forward discretization in time gives, recalling that u was initially defined as a signed distance function, |∇u| = 1. Rearranging gives, Inspection of u for u t 1 = 0 gives the location of the interface. So, z 1 = z 0 + FΔt, and the interface is moved forward appropriately. The position, z 1 , of the interface, Γ 1 , is found through solving the advection equation for u, hence an implicit tracking method. The linear level set function gives the zero dimensional point of the interface. In two spatial dimensions, the level set function is also two dimensional, with the interface taking values in a one dimensional space. This is the essence of the level set method, i.e.,

Model outline
The work flow that is implemented to carry out the consolidation simulation is shown in Fig. 2  This operator splitting approach separates the sedimentation and the compaction steps into sequential events. In order to enforce consistency between the porosity distribution and the pore pressure distribution, a Picard iteration loop was implemented over the pore pressure solution and the compaction law update. The solution to Eq. 8 is used in computing the effective stress, Eq. 1, which determines the new porosity from Eq. 15 (or similar). This is done because of the nonlinearity introduced into Eq. 8 by considering the porosity to be a function of effective stress, as in Eq. 15. By using a known value from a previous step, the Picard iteration is used to linearize the nonlinear appearance of the pore pressure in the compaction law.
In the case of one layer with a linear compaction law, the accuracy of this model can be compared to an analytical solution available in the literature, and in particular as presented by Wangen [38]. The details of which are laid out next.

Verification with analytical solution
Assuming a linear dependence of the void ratio on the effective stress and a linear dependence of the permeability on the void ratio, there is one available analytical solution for the consolidation problem considered here. Wangen [38] uses a coordinate transform, depending on the local porosity, to remove the pore space from the equation. In these "compaction-free" coordinates, the solution can be expressed analytically.
The transform between the real coordinate, z, and the Lagrangian compaction-free coordinate is This serves to scale the spatial coordinate by the porosity.
In the ζ coordinates, there is no pore space; therefore, the quantities of bulk density and bulk deposition rate need to be rendered as just the solid rock properties. Therefore, ρ b and ω b are changed to ρ s and ω s , the density of rock and the deposition rate (still length/time) of pure rock, respectively. The void ratio was chosen to scale linearly with the effective stress, e = e 0 −ασ z , where α [M −1 LT 2 ] is a compaction coefficient and e 0 [L 3 L −3 ] is the void ratio at the surface. Since e = φ 1−φ , one can see that the porosity dependence on effective stress is slightly nonlinear. For the comparison with the numerical model considered, which is formulated in terms of the porosity, it is important to note that and ∂e ∂σ z Therefore, the compressibility is The permeability is taken to be a function of the void ratio. The choice of k = k 0 1+e 1+e 0 , simplifies the pressure equation, and mimics the reduction in permeability due to a decrease in the void ratio.
With these assumptions in place and exploiting the coordinate transform, the pressure equation considered is Equation 29 is a constant coefficient diffusion equation with a constant source that has a solution of the form: where γ = (ρ s − ρ f )g and c = k 0 (1+e 0 )αμ . The coordinate transform can then be used to convert back to the real coordinates, z, using the porosity distribution computed from the pressure solution.

Dimensionless form
It is possible to perform a variable transform to nondimensional coordinates. Scaling the Lagrangian depth coordinate by ζ = ζ ω s t tot , where t tot is the total time period simulated and ω s is still the deposition rate of pure rock. Time is scaled as t = t t tot . The overpressure is scaled as Similarly, the VES is scaled as σ z = σ z (ρ s −ρ f )gω s t tot . Carrying out the variable substitions in Eq. 29 leads to, is the dimensionless diffusivity, which is time independent. One more dimensionless quantity is necessary to describe the compaction law, e = e 0 − ασ z . With the VES variable substitution, we get e = e 0 − ασ z (32) where α = α(ρ s − ρ f )gω s t tot is the dimensionless compaction coefficient.
In this formulation, the behavior of the dynamical evolution can be studied based on the two dimensionless quantities, diffusivity and compaction coefficient, thus reducing the parameter space of the of the system. In Fig. 3, across five orders of magnitude in dimensionless diffusivity, we can see the effective behavior of the overpressure, spanning from almost hydrostatic to almost lithostatic. Similarly, in Fig. 4, one can see large compaction when the overpressure is small and essentially no compaction when the overpressure is large. Thus, for a chosen dimensionless compaction coefficient, we can see the "full" range of behavior through considering the dimensionless form of the equations.
In Wangen [38], a similar transformation is carried out, though it differs in the treatment of time. Wangen considers the characteristic pressure that would be needed to see noticeable compaction. From this characteristic pressure, a characteristic length and time are derived. In particular, Wangen introduces the so-called gravity number, N = k 0 (ρ s −ρ f )g (1+e 0 )αμ and the dimensionless time, τ = α(ρ s −ρ f )gω s t. One can see that, up to the definition of time used, τ is the same as α and, furthermore, N τ = D. In [38] the formulation uses a characteristic time, thus making the scaling of the dimensionless variables time dependent. We use the total time of simulation to scale the variables, leading to a time-independent dimensionless diffusivity and compaction coefficient.

Plausibility checks
For the multi-layer case, there is no analytical solution available. Therefore, we introduce some criteria for evaluating these simulations. First, given the assumed dynamics of consolidation, except when a layer is actively being deposited, the layers should always decrease in size. This translates into a constraint on the speed function. The speed at a given layer interface should always be less than the speed at the layer interface above it, in order to be physically consistent. This monotonic decrease of the speed function with layer interface depth can be violated, for example, with a choice of a very low porosity layer above a very high porosity layer. In this case, the bounding interface below the very high porosity layer would move faster than the bounding interface below the very low porosity layer, thus leading to an increasing thickness of the very low porosity layer. Second, as in the case of the single layer, if the pore pressure becomes larger than the lithostatic pressure, leading to a negative vertical effective stress, this represents an unphysical situation for the proposed model. Physically, if the pore pressure approaches the lithostatic pressure, compaction stops and hydrofracturing, the opening of the pore space, occurs. No faulting or fracturing processes are considered in our model. Furthermore, by varying the compaction coefficients and the permeabilities of the layers, one can create conditions for overpressure buildup. A seal forms when the permeability of an overlying layer is significantly less than that of the layer beneath (e.g., an order of magnitude). The pore fluid in the lower layer has difficulty flowing upward, and thus, with continual deposition, the pore pressure increases.

Results and discussion
The reference example for the following numerical experiments is based upon a similar set-up as in [38]. The test parameters are one kilometer of uncompacted rock deposited over one million years. The single rock type is deposited with an initial porosity of 0.61, and a surface permeability of 10 −18 m 2 , with a rock density of 2720 kg m 3 . The compaction coefficient was chosen to be 5·10 −8 Pa −1 . This parameter plays a crucial role, and essentially determines the coupling strength between the vertical effective stress (VES) and the porosity. For the linear void ratio compaction law, the depth of the basin was compared to the analytical value for different values of the dimensionless quantities defined in Section 3.2.1.

Single layer with linear compaction
First we present a plot of the relevant physical variables after deposition, using a discretization of 100 time steps and a grid spacing of 2.5 m. In Fig. 5, we see that the pore pressure is bounded between the lithostatic pressure and the hydrostatic pressure, which is physically consistent. Even though the solution variable is the overpressure in the implementation, adding the hydrostatic pressure results in the pore pressure, which is plotted here.
Using the solution for the overpressure in Lagrangian coordinates as put forward in [38], one can compute the VES, and the porosity distribution. Using the coordinate transform, one then arrives back at the real coordinates to find the height of the column. Thus, for the parameters used here, the height of the column after one million years should be 907 m. The level set depth was 911 m for an error of 0.4%, which is acceptable in our opinion. The convergence of the numerical solution to the analytical solution for decreasing spatial discretization size is shown in Fig. 6.
To understand the parameter space of the simulation domain, the error of the level set height calculation for selected values of the dimensionless diffusivity and compaction coefficient are shown in Table 2.
It is necessary to consider the relationship between the diffusivity and the compaction coefficient in looking at the limits of the model using the linear compaction law. If the diffusivity is low, then the pore pressure will increase more quickly with compaction coefficient, leading to a situation where the pore pressure becomes higher than the lithostatic pressure, i.e., negative VES. On the other hand, when the overpressure is close to zero, then the VES is near its maximum and there is strong compaction. Depending on the choice of the compaction coefficient, this can lead to negative void ratios, hence invalid simulations. While that is clearly non-physical, even very small positive void ratios are outside the bounds of reasonable domain applicability of the model. The last entry in Table 2 Fig. 6 Error in the level set height calculation. Using 100 time steps per simulation, the spatial resolution, dz, was increased from 10 m to 1.25 m. The initial porosity was 0.61 and the initial permeability was 10 −18 m 2 and compaction coefficient of 5 · 10 −8 Pa −1 to a simulation where at the bottom of the sediment stack, according to the analytical solution, the void ratio is 0.0347 or 3% in porosity. (α = 1.56 corresponds to α = 2.4 · 10 −7 P a −1 , given the other parameter values.) This porosity is clearly much smaller than pure mechanical compaction alone can reach. This unphysical situation would not occur with a compaction law where the compaction coefficient changes suitably with the void ratio, e.g., an exponentional compaction law. Therefore, since such large compaction coefficients are outside the scope of physical applicability, Discretization of 2.5 m and 100 time steps used. The first five entries are the same α and D values as the analytical solution shown in Figs. 3 and 4. Note, a positive percent difference is when the level set height is bigger than the analytical solution (i.e., less compaction). The two negative percentages indicate that the numerical solution was shorter (i.e, more compaction) than the analytical solution. The last entry in the table is a simulation at the limits of the analytical solution's validity that simulation, and its large error, is largely irrelevant, and has been included for completeness. In summary, the compaction coefficient plays a fundamental role in determining the strength of the nonlinear coupling. The method is more accurate the less compaction that occurs. This coheres with the fact that there is no iteration over the coupling between the level set movement and the pressure solution. Thus, the stronger the compaction, the more the divergence. For compaction coefficients smaller than say, 1 − 5 · 10 −7 P a −1 , depending on choice of diffusivity, we conclude the model shows agreement with the analytical solution. This includes the regime of typical physically reasonable lithological compaction coefficients.

Multi-layer with linear compaction
In order to perform a simulation with four rock types, at least four level sets are required. Each interface is assigned its own level set potential function, with the 0 value defining the location of the interface. When considering four rock types, we have the three interfaces between them and one between the "basement" and the deepest rock. Each level set is solved independently, though the same speed function relation, Eq. 22, is shared among them. For a consistency check, we compare the multiple layers, initialized all with the same rock type, in order to compare to the single layer case. For the setup with an analytical solution of 907 m, the multiple layer, single rock numerical model simulates 912 m, compared to 911 m for the explicitly single layer case above. The physical variables are plotted in Fig. 7. We conclude this is successful for the multiple layer implementation.

Multi-layer with linear compaction
Introducing a heterogeneous permeability distribution in the column, (Fig. 8c) leads to differences in the pore pressure shown in Fig. 8. If a layer of significantly lower permeability is above a layer of higher permeability, then overpressure will build up, as the pore fluid is trapped. We can see that the pore pressure approaches the lithostatic pressure (overburden) under the seal, located at circa 500 m, which comes from a 2 order of magnitude permeability difference between the two middle layers. This contrasts to to the pore pressure profile in Fig. 7, where there is a smooth permeability profile.
The different layers can also have different compaction coefficients. The differences in compaction are clearly seen in the void ratio plot in Fig. 9. We can see in Fig. 10 how the four layers grow during their period of deposition and how the layers are reduced at different rates, depending on their Fig. 7 Multiple layers, single rock type: linear compaction law. A uniform compaction coefficient of 5 · 10 −8 P a −1 was used.The initial porosity was 0.61 and the initial permeability was 10 −18 m 2 with a discretization of 2.5 m compaction coefficient. In the simulation, the second layer from the bottom has a larger compaction coefficient than the other layers.
The thicknesses of the different layers are given in Fig. 10. The second layer ends up with the smallest thickness, as we might expect given the compaction coefficients. The oscillations that can be seen in the layer thicknesses come from the relative definition of the layers. As the 0 level  set of a layer approaches a node, it might reach it slightly before or after another layer's 0 level set reaches its next node. The size of the oscillations are a discretization effect dependent on the ratio of length of rock added per time step over the spatial step size. Over the simulation time, the trend From the first layer on the left, the compaction coefficient is 5 · 10 −8 P a −1 , 9 · 10 −8 P a −1 , 3 · 10 −8 P a −1 and 5 · 10 −8 P a −1 clearly shows a reduction in thickness. This decrease looks plausible when considering the compaction coefficients and permeabilities of the different layers.

Nonlinear compaction law
More realistic situations can be approached by modelling compaction using Athy's law, Eq. 15. As there is no known analytic solution for the consolidation problem using this compaction relation, we present the results of a simulation of multiple layers of one rock type to show the consistency of our modelling approach. In Fig. 11, we see the effect of the exponential compaction law. For comparison, the analytical solution, of course with linear compaction, using the same parameters shows the significant increase in compaction. The total length of the column after the simulation is over 100 m shorter than when the linear compaction law is used. For similar pore pressures and effective stresses, the geometrical difference as a result of the compaction law used is significant.
As a concluding test case, we demonstrate results from a simulation using Athy's law with layers defined by distinct depositional periods, initial porosities, permeabilities and compaction coefficients. The simulation parameters are included in Table 3. Looking at Figs. 12 and 13, one can see the fundamental interplay between compaction coefficient and permeability. Layer 3 compacts significantly, despite having the smallest compaction coefficient, since Fig. 11 Multiple layers, single rock type: Athy's compaction law. Using the multiple layer framework, a uniform compaction coefficient of 5 · 10 −8 Pa −1 was used. The initial porosity was 0.61 and the initial permeability was 10 −18 m 2 with a discretization of 2.5 m. As there is not an analytic expression in this case, the analytic solution for the linear compaction law is plotted for comparison only φ 0 is initial porosity, k 0 is initial permeability, β is compaction coefficient layer 4 above it has a large permeability. Conversely, layer 2 has a low permeability, which stops layer 1 from compacting significantly even with its high compaction coefficient. Figure 14 provides a burial plot of the simulation. In closing, the usefulness of the level set method in representing a sedimentary basin will be considered. By design, the level sets introduce an additional, though simple, advection equation for each interface to be tracked. Of course, this adds to the computational cost. However, there are two primary features: (1) simplicity of the grid and (2) flexibility of the speed function. Firstly, in adopting the Eulerian framework, the grid remains static. There is no need to regrid during the computation as the sedimentary layers evolve. The numerical advantages of a regular grid can be maintained even for larger and longer computations, wellsuited for high performance computing. Secondly, the speed function which controls the evolution of the fronts can be modified ad hoc. It is possible to consider the physical processes one wishes. The form of the coupling of the physics to the geometry remains the same through the level set Fig. 12 Multiple layers, multiple rock types: Athy's compaction law. Description of rock properties is found in Table 3 Fig. 13 Multiple layers, multiple rock types: Athy's compaction law. Thickness of layers in time. Description of rock properties is found in Table 3 evolution. This assumes that the physics solved on the regular grid can be sufficiently well-represented on the static grid for all subsequent material arrangements. Burial history. Layer 1 has a high compaction coefficient and layer 2 has a low permeability, while layer 3 has a low compaction coefficient and layer 4 has a high permeability. Description of rock properties is found in Table 3

Conclusion
In this paper, a model of single phase fluid flow under isothermal conditions is coupled to the movement of porous rock undergoing compaction due to sedimentation over long time scales. The novel ingredient is the use of the level set method to track the interfaces between rock types while keeping the numerical discretization grid unchanged. We proposed a speed function, dependent on sedimentation and compaction. Numerical experiments have been conducted to verify the model against the one available analytical solution, describing a restricted case of a single rock layer with a linear compaction law. In the more realistic cases of greater interest to basin modelling, namely Athy's nonlinear compaction law and multiple rock layers, we have presented results from simulations that are consistent with the plausibility checks discussed. We have shown the flexibility of the speed function. Once the numerical model is set up, with the level set equations and speed function, we could change the compaction behavior by modifying only the compaction law used (linear in the void ratio or Athy's law) while the rest of the code was unchanged. It is also possible to add other terms into the speed function, such as for erosion or cementation.
Looking forward, in order to directly address large deformations, such as those due to folding or buoyancy effects, a multidimensional model that includes lateral stresses is needed. This would mean relaxing the vertical effective stress assumption. How the effective stress tensor influences the interface motion is determined through the speed function. It remains open how to construct a speed function which would represent the strains due to these stresses. Further elaborations could also include introducing heat flow and the temperature dependence of viscosity and density.
One area where the level set coupling model could be advantageous is overthrusting. Overthrusting and lateral stresses lead to a modification of the vertical effective stress into a mean effective stress. Thus, allowing lateral movement, the multiple layer occurrences that result are treated in practice with a domain decomposition into blocks. In principle, the topological flexibility of the implicit surface method would not require this block concept for the thrust belts.