{'fraction_non_alphanumeric': 0.034906588003933134, 'fraction_numerical': 0.00983284169124877, 'mean_word_length': 4.639787430683919, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Space Time Domain Decomposition Approach using Enhanced Velocity Mixed Finite Element Method

A space-time domain decomposition approach is presented as a natural extension of the enhanced velocity mixed finite element (EVMFE) [Wheeler et. al] for spatial domain decomposition. The proposed approach allows for different space-time discretizations on non-overlapping, subdomains by enforcing a mass continuity argument at the non-matching interface to preserve the local mass conservation property inherent to the mixed finite element methods. To this effect, we consider three different model formulations: (1) a linear single phase flow problem, (2) a non-linear slightly compressible flow and tracer transport, and (3) a non-linear slightly compressible, multiphase flow and transport. We also present a numerical solution algorithm for the proposed domain decomposition approach where a monolithic (fully coupled in space and time) system is constructed that does not require subdomain iterations. This space-time EVMFE method accurately resolves advection-diffusion transport features, in a heterogeneous medium, while circumventing non-linear solver convergence issues associated with large time-step sizes for non-linear problems. Numerical results are presented for the aforementioned, three, model formulations to demonstrate the applicability of this approach to a general class of problems in flow and transport in porous media.


Introduction
Complex multiphase flow and reactive transport in porous medium is mathematically modeled using a system of non-linear partial differential equations. In fact, a multitude of physical processes in the porous medium are modeled by varying the nature of these non-linearities. For example, in multiphase flow the phase relative permeabilities are modeled as non-linear functions of saturation [17]. On the other hand, polymer flow is modeled with an apparent viscosity term dependent on phase velocities (or phase pressure gradient) [15]. Each of these aforementioned two examples introduces a non-linear term, function of saturation or pressure, as a coefficient in the otherwise linear relationship between Darcy velocity and gradient of pressure. A consistent spatial and temporal discretization of these PDEs result in a non-linear algebraic system of equations that can be linearized using Newton method. In order to accurately resolve these non-linear terms it is often necessary that a small time-step size be used during numerical solve. For example, the timescale for reactive transport with kinetic reactions controls the flow and transport timescales [11]. This requires that small timestep sizes be used for the entire computational domain, to be consistent, wherein a majority of the reactions are restricted to only a subdomain. Consequently, these conventional approaches are often computationally prohibitive.
Several space-time domain decomposition approaches have been proposed in the past that allow non-matching spatial and temporal subdomain discretizations to address these problems. [8,9] proposed space-time finite element methods for elastodynamics with discontinuous Galerkin (DG) discretization in time. [6] present two, space-time domain decomposition methods for a parabolic, diffusion problem using a mixed formulation differing in treatment of the non-matching subdomain interface condition. This was later extended to handle advection-diffusion type problems in [7] using operator splitting to handle the advection and diffusion terms separately. Another space-time domain decomposition approach for parabolic problems was proposed by [1]. Here, a comparison between a mortar mixed finite element method and finite volume method with ad hoc projections for coupling the non-overlapping, space-time, subdomain problems was also presented. A number of variational space-time methods for elastic wave propagation and diffusion equations were also proposed by [10,2] using DG in space with discontinuous, continuous, and continuous differentiable Galerkin methods in time.
In this work, we restrict ourselves to flow and transport problems in subsurface porous media with an intent to develop this approach for a more general class of problems in the near future. We extend the EVMFE method [20], for spatial domain decomposition, to encompass the temporal dimension and propose an enhanced velocity space-time domain decomposition approach. The EVMFE method for spatial domain decomposition has been applied to a wide range of flow and transport problems [17,18,19] in porous medium such as slightly and fully compressible multiphase flow as well as equation of state compositional flow. Each of these aforementioned problems differ only in the nature of the non-linearities introduced by the model formulation representing the physical process.
These non-linearities occur as relative permeability and capillary pressure defined as functions of saturations, local equilibrium problem associated with phase behavior or reactive equilibrium, as well as non-linear source/sink terms in the case of reactive kinetics.
We begin by describing a linear, single-phase flow model formulation, in section 2, with appropriate initial and boundary conditions. We then describe a mixed weak variational formulation with an appropriate choice of functional spaces for the pressure and velocity unknowns. We require that the velocity space be a subspace of H(div) for local mass conservation at the non-matching space-time interface by enhancing the trace of the velocity space both spatially and temporally. With this enhanced space-time velocity space we then describe a mixed weak variational problem problem for the proposed space-time domain decomposition approach. This is followed by a fully discrete form, in section 3, of the aforementioned weak variational form to derive the algebraic equations associated with the model formulation. A numerical solution algorithm is then presented in section 4.
Here, we present a space-time monolithic solver that eliminates the need to iterate between the two subdomain problems with the interface condition serving as an evolving internal boundary condition for each of the subdomains. We then present a preliminary numerical convergence analysis for the proposed scheme in section 5. This is followed by numerical results in section 6 where two practical problems of interest to flow and transport in subsurface porous media communities are considered.

Single phase flow formulation
We consider the following single-phase flow model; for simplicity of discussion, represented by a first order parabolic partial differential equation system, Here, J = (0, T ] is the time domain of interest, follows, d is used to indicate spatial dimensions with d + 1 as the temporal dimension. A Dirichlet boundary condition is considered for the sake of simplicity of description and more general boundary conditions can also be treated. The subdomains Ω i × J m are non-overlapping. Let us define Γ m,n i,j = (∂Ω i × J m ) ∩ (∂Ω j × J n ) as the interface between the space-time subdomains Ω i × J m and Ω j × J n .
are the interior interfaces. Note that this description of space-time subdomain interfaces does not take into account interfaces normal to the time-direction. Since there are no arguments for pressure regularity in the time-direction, no special considerations are required.
The functional spaces for the mixed weak formulation of Eqns. (1)-(4) are, The mixed weak formulation of Eqns. (1)-(4) is the pair u ∈ V , p ∈ W such that, Here, ν is the normal in the spatial direction. We use the lowest order Raviart-Thomas space The RT 0 spaces for any element E ∈ T t h are defined follows: Note that the definition of element E is in a space-time sense and hence the measure of E inherits this sense accordingly. The degrees of freedom for the vector v ∈ V t h (E) can be specified by the values of the normal component v · ν at the midpoint of all edges (faces) of E in the spatial direction.
Here, ν is the outward unit normal vector; along the spatial dimensions, on ∂E. The pressure finite element space on Ω × J is defined as, As in the original work on enhanced velocity mixed finite element method, since W h is a discontinuous space, no special consideration is needed on the interfaces. We want to construct a velocity finite element space V t, * h ⊂ V on the multiblock partition T t h of (Ω × J). Let, for m = 1, ..., q and i = 1, ..., r, however is not a subspace of H(div; Ω × J) since the normal vector components do not match on Γ.
We therefore need to modify the degrees of freedom on Γ. Let T t,m,n h,i,j be the rectangular partition of T m,n i,j obtained from the intersection of the traces of T t,m h,i and T t,n h,j . Again, as in the original work, we assume that the subdomains do not meet at an angle for a three dimensional problem. We then force the fluxes to match on each element e ∈ T t,m,n h,i,j . Consider any element E ∈ T t,m h,i such that E ∩ Γ m,n i,j = ∅. The interface grid divides the boundary edge (face) of E. This division can be extended inside the element in the space and time sense as shown in Figure 1 (right). On each subelement E m i we define a basis function v E l k in the RT 0 space V h (E l k ) which has a normal component in each spatial dimension v E l k ν equal to one on e l k and zero on the other edges. Let V t,Γ h be the span of all such basis functions. We then define the space-time multiblock mixed finite element space to be, With the space-time enhanced velocity and piecewise constant pressure spaces defined, described above, the mixed finite element discretization of Eqns. (1)-(4) is: find u t h ∈ V t, * h and p t h ∈ W t h such that,

Fully Discrete Form
In this section, we present the fully discrete form of the mixed variation formulation discussed in section 2. We begin by describing the RT 0 × DG 0 basis functions as piecewise constants for pressure and saturation (or concentration) unknowns and piecewise linear in space for velocity (or flux) unknowns as follows, Here, e n j+ 1 2 is an edge of a space-time element.
The quadrature rules remain the same in space, as for the original enhanced velocity mixed finite element method [20], augmented only by a midpoint rule in the time dimension as follows, As mentioned before, d = 1, 2, 3 only represents the spatial dimension of the problem with d + 1 as the temporal dimension. For the ease of describing the fully discrete form below below we choose one spatial dimension d = 1 with d + 1 as the temporal dimension. The solution of Eqns. (1)- (4) can then be written as, We will now construct an algebraic system of equations by testing the variational forms of the discrete constitutive and conservation equations with w n j and ϕ n j+ 1 2 , respectively. Here, we evaluate most of the integral terms in the variational problem on a matching grid, bifurcating to the non-matching grid only for the integrals where the non-matching, space-time interface enters the evaluation. This allows us to easily delineate the differences arising in the algebraic system of equations due to the choice of our discretization. For the constitutive equation we have, It is easy to see that, The second term in Eqn. (15) can be expanded as, Let us now consider, for a given j 0 ∈ I, a non-matching grid with fine a domain at j 0 + 1 2 − and a coarse domain at Here,Ẽ n j0+ 1 2 is a subelement of the coarse elementẼ n and ϕ n−1

and,
p, ∇ · ϕ n−1 respectively. Finally, the boundary or the third term in Eqn. (15) for the coarse and fine domains is given by, For the conservation equation testing with w n j we get, The first term in Eqn. (23) can be expanded as, For the coarse domain (j 0 + 1 2 ) + , the above integral remains unchanged. However, for the fine The second term in Eqn. (23) is, For a fine domain element, with an edge at the the non-matching space time interface, we can write this term as, Similarly, for a coarse domain element, The third term in Eqn. (23) is easy to evaluate for both the coarse and fine domains as, We now have a system of algebraic equations from the fully discrete form of the variational problem (7)- (8). For a non-linear, multiphase flow problems the above fully discrete form results in a nonlinear, algebraic, system of equations.

Solution Algorithm
We propose a space-time monolithic solver for the algebraic system obtained in the previous section. Here we first define define matching times as the location along the time dimension where the coarse and fine domain boundaries (∂Ω c ∩ ∂Ω f ) overlap such that ∆tc ∆t f = l ∈ I, where l is a non-zero integer. Figure 2 shows a schematic of the space-time domain decomposition used for the construction of the monolithic system with l = 3.  . . .
Note that here the flux unknowns δu are eliminated by taking a Schur-complement of the original linear algebraic system in pressure and flux unknowns δp and δu, respectively. We therefore avoid the saddle-point system associated with the original, linear algebraic system. Although not restrictive, for the ease of description we assumed non-matching grids in the time-dimension only with matching grids in space for the coarse and fine domains. The spatial sub-matrix then has a known sparsity pattern of three, five, and seven non-zero diagonals for one, two, and three spatial dimensions, respectively for RT0 mixed finite element discretization in space. As for the original EVMFEM scheme in space, the spatial sparsity pattern alters when non-matching spatial grids are considered.
The temporal sub-matrix is always diagonal for DG0 discretization in time.
Since DG0 in time is closely related to the backward Euler scheme [1], the solution scheme presented here is fully-implicit in the space-time unknowns. Figure 3 shows this fully-implicit solution algorithm for the space-time domain decomposition approach proposed in this work. Here, n and k are the time-step and Newton iteration counters, respectively. For the single phase flow problem (1)-(4), a single Newton iteration (k=1) is required for the linear system to converge, as expected.
For non-linear multiphase flow and transport problems, we require that the max norm of the nonlinear |R nl | residuals be less than a desired tolerance . We initially rely upon direct-solvers for the purpose of testing and benchmarking the solution algorithm. However, the monolithic spacetime solver allows us to utilize parallel in time linear solvers and preconditioners presented in [5] with relative ease. This also renders us a massively parallel, time-concurrent, framework for solving general, sub-surface, non-linear flow and reactive transport problems.

Numerical convergence analysis
In this section, we consider a linear parabolic problem to verify numerical convergence of the proposed space-time domain decomposition approach as described below, Here, I is the identity matrix with other symbols with their usual meanings. Further, the boundary and initial conditions are chosen to be, and respectively. The following choice of the forcing function, then allows us to obtain an analytical form of the solution as, p = e c1t sin(2πx)sin(2πy).
An L 2 norm in space and time for the error in pressure ||p − p h || L 2 is defined as follows, We select a unit domain with two spatial (d = 2) and one temporal (d+1 = 3) dimension of size 1 × 1 × 1. The fine subdomain is refined by a factor of 4 with respect to the coarse subdomain, both spatially and temporally. Equipped with the above norm, we compute the errors in the coarse and fine subdomains for increasing mesh refinement while maintaining the aforementioned refinement factor of 4 between the coarse and fine subdomains. Here, we only show a sense of convergence for the proposed scheme and reserve optimal convergence rate studies using appropriate norms for a future work along with derivation of rigorous a-priori error estimates.  Table 1 shows the coarse and fine subdomain solution errors for increasing mesh refinement.
Here, h f and h c are the fine and coarse subdomain h-refinements, respectively in both spatial and temporal dimensions. As expected, for this space-time monolithic, linear system at hand only one linear solve is required to reach the solution. Table 1 also shows the degrees of freedom associated with each of the three cases and the corresponding CPU time. Figure 4 show the pressure solutions for the linear parabolic problem at the three refinement levels.

Single Phase Flow and Transport Problem
We first present a numerical experiment for a slightly compressible, single phase flow and nonreactive tracer transport in porous medium. The transported component is assumed to be a tracer suchthat the fluid density remains invariant with changes in tracer concentration. The model formulation for a single phase, slightly compressible flow and tracer transport in porous medium is given by, with the Darcy velocity u and diffusive flux j defined as, , respectively. Further the boundary and initial conditions are given by, , respectively. Here, u and z are the Darcy and diffusive fluxes, respectively, ρ the density, c and p the pressure and concentration unknowns, respectively, D is the constant (spatial and temporal) diffusion coefficient, q the source/sink term,ĉ the injection concentration at the source, ν the unit outward normal, K is the spatially varying diagonal permeability tensor, µ the viscosity, φ the porosity, and p 0 and c 0 the initial conditions for pressure and concentration, respectively. The density is defined as a non-linear function of pressure given by, Here, ρ ref is the reference density at the reference pressure p ref , and c is the fluid compressibility.
The single phase flow and tracer transport model formulation presented here assumes a constant molecular diffusion for ease of description of a fully discrete mixed formulation. However, a spatially and temporally varying diagonal tensor can also be considered as the diffusion coefficient to account for hydrodynamic dispersion. A space-time, fully discrete form of the above model formulation is described in Appendix A. The fluid compressibility, viscosity, and density are taken to be 10 −6 psi −1 ,1cP, and 64lb/f t 3 , respectively.
The computational domain is (110ft × 30ft × 1ft × 100days) with the fine and coarse subdomains discretized using grid elements of size 0.5ft × 0.5ft × 1ft× 1day and 5ft × 5ft × 1ft × 5days, respectively. The fine subdomain is refined by a factor of 10 and 5 times with respect to the coarse subdomain in space and time, respectively. A source and sink term is considered at the bottom left and top right corners of the domain using an rate specified injection and pressure specified production wells, respectively. The injection well is water-rate specified at 4 STB/day with a non-reactive tracer injection at concentration 1. The production well is pressure specified at 1000 psi. A homogeneous porosity and permeability of 0.2 and 50 mD, respectively are assumed for the entire domain with a constant diffusion coefficient of 0.1. Further, the initial reservoir pressure and concentration are taken to be 1000 psi and 0.0, respectively.

Two Phase Flow Problem
We describe the immiscible, two-phase, slightly compressible flow in porous medium model formulation with oil and water phase mass conservation and constitutive equations along with the boundary and initial conditions. The fully discrete space-time formulation, similar to single phase slightly compressible flow in the previous subsection, is presented in Appendix B. Further details regarding the two-phase flow model formulation can be found in [16]. The mass conservation equation for phase α is given by, where φ and K have their usual meanings as described before, and ρ α , s α , u α and q α are density, saturation, velocity and source/sink term, respectively of phase α. The constitutive equation for the corresponding phase α is given by Darcy's law as, Further, k rα , µ α and p α are the relative permeability, viscosity and pressure of phase α. Although not restrictive, for the sake of simplicity we assume no flow boundary conditions.
Here, p 0 α , s 0 α are the initial conditions for pressure and saturation of phase α. Furthermore, the phase saturations s α obey the constraint, We assume capillary pressure and relative permeabilities to be continuous and monotonic functions of phase saturations, The oil and water phase are assumed to slightly compressible with phase densities evaluated using, Here, c f α is the compressibility and ρ α,ref is the density of phase α at the reference pressure p α,ref .
For the numerical experiment, the computational domain is kept the same as before (110ft × 30ft × 1ft × 40days) with the fine and coarse subdomains discretized using grid elements of size 0.5ft × 0.5ft × 1ft × 1day and 5ft × 5ft × 1ft × 5days, respectively. The fine subdomain is refined by a factor of 10 and 5 times with respect to the coarse subdomain in space and time, respectively. The fluid and reservoir properties are adapted from the SPE10 [4] dataset with an assumed homogeneous, spatial distribution for porosity of 0.2. The log (natural log) scale, permeability distribution in the horizontal (y) and vertical (x) directions are shown in Figure 6. The oil and water phase compressibilities are taken to be 1×10 −4 and 3×10 −6 psi −1 , respectively, and densities 53lb/f t 3 and 64lb/f t 3 , respectively. Further, the fluid viscosities are assumed to be 3 and 1 cP for the oil and water phases, respectively. Additionally, a Brook's Corey model [3], Eqn. (52), is considered for the two-phase relative permeabilities with endpoints s or = s wirr = 0.2 and k 0 ro = k 0 rw = 1.0, and model exponents n o = n w = 2. Figure 7 shows the relative permeability and capillary pressure curves as functions of saturation.
The capillary pressure function is defined using the van Genuchten model [12] given by Eqn.

Conclusions
A space-time, domain decomposition approach is presented for flow and transport problems in subsurface porous medium using an enhanced velocity approach. We also described an efficient spacetime monolithic solver for this approach that does not require subdomain iterations for computational efficiency. In fact, any linear problem requires only one linear solve to reach the solution, as expected.
A preliminary numerical convergence analysis for a linear parabolic problem allows us to confirm that the proposed approach converges as the mesh is refined in both space and time. A rigorous derivation of the a-priori error estimates equipped with appropriate norms in order to determine exact convergence rates is reserved for a future work. We presented two, non-linear flow and transport problems differing in the degree of non-linearity to demonstrate the general applicability of this spacetime, domain decomposition approach to a wide range of problems in subsurface porous medium.
The two-phase, non-linear flow problem with space-time, mesh refinement in the vicinity of injection well was able to circumvent non-linear solver convergence issues without imposing restrictively small time-step sizes over the entire computational domain. Further, the proposed domain decomposition approach along with the numerical solution algorithm renders us a massively parallel, concurrent in time, computational framework allowing us to use parallel in time, linear solvers and preconditioners without compromising computational efficiency.
that differ from the description in Section 3. The first term in the phase mass conservation Eqn.
(A.1) can be expanded as, As before, for the coarse domain (j 0 + 1 2 ) + , the above integral remains unchanged. However, for the fine domain (j 0 + 1 2 ) − we have, Similarly, the accumulation (first) term in the transport equation Eqn. (A.2) is expanded for the coarse and fine subdomains as, , and ∂ ∂t φρc, w  ∂ ∂t φρ w s t w,h , w + ∇ · u t w,h , w = (q w , w) , (B.2) with w ∈ W and v ∈ V . Please note that s o and p w are eliminated in the above formulation using the algebraic constraints Eqns. (48) and (49), respectively. Further, λ α is defined as the mobility of phase α as, An expanded mixed formulation [13,14], with additional auxiliary phase fluxesũ α , is used to avoid inverting zero phase relative permeabilities (k rα ). As before, the solution can be written as,