Dynamic adaptive mesh optimisation for immiscible viscous fingering

Immiscible fingering is challenging to model since it requires a very fine mesh for the numerical method to capture the interaction of the shock front with the capillary pressure. This can result in computationally intensive simulations if a fixed mesh is used. We apply a higher order conservative dynamic adaptive mesh optimisation (DAMO) technique, to model immiscible viscous fingering in porous media. We show that the approach accurately captures the development and growth of the interfacial instability. Convergence is demonstrated under grid refinement with capillary pressure for both a fixed unstructured mesh and with DAMO. Using DAMO leads to significantly reduced computational cost compared to the equivalent fixed mesh simulations. We also present the late-time response of viscous fingers through numerical examples in a 2D rectangular domain and in a 3D cylindrical geometry. Both problems are computationally challenging in the absence of DAMO. The dynamic adaptive problem requires up to 36 times fewer elements than the prohibitively expensive fixed mesh solution, with the computational cost reduced accordingly.


Introduction
The term viscous fingering refers to the unstable displacement of one fluid by another in a porous medium. This instability was first described by Saffman and Taylor [1] and has since attracted considerable attention in a variety of different applications. It can occur in both miscible [2][3][4] and immiscible [5][6][7][8][9] displacements, although the behaviour and modelling of miscible viscous fingering has received considerably more attention in the literature. This is because miscible viscous fingering depends on relatively few parameters (viscosity ratio and diffusion/dispersion) and thus can be more easily analysed mathematically and modelled numerically. Immiscible viscous fingering depends upon the A. E. Kampitsis a.kampitsis@imperial.ac.uk 1 Novel Reservoir Modelling and Simulation Group, Department of Earth Science and Engineering, Imperial College London, Prince Consort Road, London, SW7 2BP, UK relative permeabilities of the fluids as well as the viscosity ratio and instead of diffusion/dispersion depends upon capillary pressure, which in turn depends upon saturation.
Although immiscible viscous fingering is less well studied, it is important in a variety of subsurface applications including enhanced oil recovery (water displacing viscous oil, immiscible gas injection) [10] and carbon dioxide sequestration [11,12]. In enhanced oil recovery, it can result in reduced sweep efficiency meaning that oil is bypassed by the injected fluid whilst in CO 2 sequestration it may result in the CO 2 plume travelling further than expected through the formation, potentially reducing the security of the storage [13]. Conversely, viscous fingering in conjunction with gravity may enhance the dissolution of the CO 2 in the aquifer, increasing the security of the storage [14]. Due to the strongly non-linear behaviour of viscous fingers, the only way to investigate the potential impact of fingering on oil recovery or security of CO 2 storage is via numerical simulation, although it is possible to derive analytical solutions to describe the very early-time behaviour (e.g. [5] and [15]).
The numerical modelling of immiscible viscous fingering is challenging because the dynamics depend on the balance between the finger growth (which is driven by the viscosity ratio of the fluids and the relative permeability functions) and capillary pressure. To capture the correct fingering pattern, the numerical simulation must be dominated by capillary pressure rather than numerical diffusion and dispersion. Consequently, very fine meshes are required which make the simulations computationally intensive, especially when exploring the late-time behaviour in realistic 3D displacements.
Over the last two decades, the continuum equations that govern viscous fingering have been solved in the literature using many different numerical methods including finite volume [4,16], spectral [17,18], and continuous and discontinuous Galerkin finite element (FEM) [19][20][21][22] as well as mixed control volume finite element (CVFEM) [23,24]. Early work focussed on the use of higher order numerical schemes in association with finite volume methods to ensure that physical diffusion dominated over numerical diffusion [4,16]. However, these schemes required a fine mesh resolution across the whole physical domain making the simulation of fingering patterns computationally expensive, particularly in 3D. In addition, they increase grid orientation error [25] which does not reduce with mesh refinement. Grid orientation errors tend to be smaller when a CVFEM approach is used in conjunction with an unstructured mesh and appropriate initial conditions; however, the fine meshes required to obtain converged solutions can still make these approaches prohibitively expensive to apply.
Adaptive mesh refinement (AMR) has the potential to reduce computational effort when simulating viscous fingering by increasing mesh resolution around the fingers to ensure that capillary pressures dominate over numerical diffusion and dispersion whilst coarsening the mesh in other parts of the domain. Initial applications to viscous fingering involved the use of an underlying, fixed, structured mesh which was locally refined or coarsened depending on specified error metrics. These enabled the modelling of viscous fingering in both miscible [26,27] and immiscible [28] 2D systems. More recently, Lee and Wheeler [29] used adaptive enriched Galerkin methods on structured Cartesian grids to model miscible fingering in linear and radial displacements. They presented results from large 3D simulations and argued that AMR reduced computational cost by allowing resolution to be focussed exclusively along finger tips but did not present any data to support this assertion.
Dynamic adaptive mesh optimisation (DAMO) is an improvement over AMR in that it allows free movement and repositioning of nodes in an arbitrary unstructured mesh as well as bisection/joining of existing elements. It is thus potentially less prone to the grid orientation errors associated with structured meshes and discussed above. It has been used for many years in computational fluid dynamics (CFD) [30] but has only recently been applied to porous media flows [24,31,32]. Both [24] and [32] applied DAMO to the modelling of immiscible viscous fingering in 2D rectangular systems but neither included capillary pressure and so were unable to demonstrate convergence under mesh refinement. Adam et al. [24] showed that DAMO could reproduce the pattern seen in a very fine, fixed mesh simulation of one immiscible viscous fingering case whilst [32] demonstrated that DAMO using an unstructured mesh gave similar results to finite volume simulations using a fixed Cartesian mesh but the simulations were less affected by grid orientation errors. More recently, Abdul Hamid et al. [33] showed that the DAMO model used by [32] predicted the same early-time growth rates of viscous fingers as predicted by two different, finite volume numerical models. All these agreed well with the values predicted by perturbation analysis for very small wavenumbers, thus providing some validation of the approach.
In this paper, we significantly extend what has been done previously in the literature by applying DAMO to the realistic problem of immiscible viscous fingering with capillary pressure in 2-and 3D domains. We build upon some preliminary results presented by Adam et al. [34], using similar data sets but using the double control volume finite element method (DCVFEM) [35], which allows us to use a higher order representation for velocity and pressure in conjunction with the adaptive mesh. These enable the model to achieve converged solutions with fewer elements compared to the results presented in [34]. Higher order interpolation is used to minimise the growth of errors during each mesh adapt [24]. To illustrate the practical utility of unstructured adaptivity, we focus on fingering simulations that are very challenging to carry out using fixed mesh approaches due to the prohibitively high computational cost. We demonstrate convergence of viscous fingers in a 2D simulation where the non-wetting phase is displaced by the wetting phase using both fixed and adaptive meshes. We quantify the speed-up due to mesh adaptivity. DAMO allows us to consider two otherwise highly intensive simulations, a late-time simulation of a 2D displacement in a high aspect ratio domain as well as a 3D displacement in a cylindrical domain that is, to the best of our knowledge, the first of its kind.

Governing equations
The governing equation for multiphase flow in porous media is the generalised form of Darcy's law, written for a phase α as where the subscript α = w, nw indexes the wetting and non-wetting phases respectively, q α is the volumetric fluid flux, K r α is the relative permeability, K is the permeability tensor, ρ α is the density and μ α the viscosity of phase α.
In the formulation with capillary pressure, we write Darcy's law in a slightly modified form as where v is the force density, u is the phase saturationweighted Darcy velocity, p and p c are the pressure and the capillary pressure, respectively, and σ α is defined for phase α as with S α being the saturation (phase volume fraction) of phase α. Gravitational effects are assumed to be negligible in the cases studied in this paper, and hence the gravitational terms have been omitted from Eqs. 2 and 3. The saturation equation for incompressible flow is written as where φ is the porosity of the medium. Finally, in order to ensure a closed system of equations, the saturation constraint is imposed as To discretise the governing equations, the DCVFEM [35] is used in this paper, which is a variation of the commonly used CVFEM approach. The main improvement in DCVFEM is that although the velocity is discretised using finite elements as in CVFEM, the pressure and saturation are discretised using the control volume mesh. Both saturation and pressure are expanded using the same shape function ensuring a consistent representation. Therefore, fields calculated from the saturation but resolved in the pressure space, such as the capillary pressure, are consistently represented.
The DCVFEM improves the quality of the pressure matrix for the highly distorted meshes often required in high aspect ratio domains found in subsurface reservoir problems [36]. Moreover, using the DCVFEM enables us to use large-angle elements at the displacement fronts when using DAMO. Consequently, fewer elements are used leading to the reduction of the computational demands and ultimately allowing solutions to be obtained for systems where the classical approach fails. Further details of the discretisation method can be found in [31,35,37].
In this paper, the element pair P n−1 DGP n CV is used for the numerical simulations. P n−1 refers to the polynomial order of velocity discretisation, DG denotes the use of discontinuous Galerkin, P n refers to the order of discretisation of the pressure and CV stands for the use of the control volume shape functions. The θ-method is used for time discretisation, where θ smoothly varies between 0.5 (Crank-Nicolson) and 1 (implicit Euler) based on a total variation diminishing criterion [38]. The numerical method presented here is implemented in the open-source code IC-FERST (Imperial College Finite Element Reservoir Simulator).

Dynamic adaptive mesh optimisation
Dynamic adaptive mesh optimisation is a method of automatically refining the computational mesh where properties are changing rapidly in space and coarsening elsewhere. An overview of the method and its applications in computational fluid dynamics over the last decade can be found in [30].
The DAMO approach applied in this paper utilises the anisotropic mesh optimisation techniques presented in [39], in which elements edges may collapse, split or swap and element vertices may be moved. It can be shown that the interpolation error between a smooth field ψ(x i ) and its linear interpolation over a given fixed mesh ψ = i q i N i , where N i are the finite element basis functions and q i the nodal values of the field, is bounded by a function of the Hessian matrix [39] Mesh adaptivity proceeds by constructing a functional I dependent on this interpolation error bound that measures mesh quality for a given domain [39]: where and v i are vectors describing the element edge lengths on the finite element mesh, is a normalisation constant, γ is the polynomial degree of the finite element interpolation and δ is the number of dimensions in the problem. The mesh adaptivity process therefore amounts to minimising the functional I and hence generating a mesh with a minimum interpolation error estimate. Mesh adaptivity also accounts for other possible constraints such as the geometric configuration of the problem at hand, the maximum required number of elements, and the mesh anisotropy or gradation. Mesh anisotropy refers to the maximum aspect ratio of the element's edge-lengths, whilst gradation refers Allowed mesh adaptivity operations [39] to the variation of the size on consecutive elements, i.e. it controls how fast the mesh size may change. The new mesh is obtained by applying the following techniques ( Fig. 1) to the original mesh, provided that the defined error criterion is improved: 1. Edge splitting: refinement via splitting existing elements and adding an additional node along an existing edge and regenerating the elements which share it 2. Edge collapsing: coarsening via removing an existing node by collapsing an existing edge to zero length and thus replacing two nodes by a single one lying at the edge midpoint 3. Edge and face-edge swapping: reordering the connectivity of existing elements and introducing an edge between the two nodes of two elements that are not shared 4. Node movement: repositioning nodes within the convex hull spanned by the elements which share it, to improve mesh quality For the CV fields (pressure and saturation) a CV-Galerkin interpolation is used to map the data from one mesh to another. CV-Galerkin is a three-step Galerkin technique [40], fully conservative and bounded, and has second-order re-mapping [24]. The steps of the CV-Galerkin interpolation method are summarised as -Mapping of the CV field on the old mesh onto a FE representation via Galerkin projection -Mapping of the resulting FE mesh to a supermesh (i.e. the intersection of the old and target meshes) via FEM Galerkin projection -Project the new FE mesh back into a CV representation on the new mesh For the FE field (velocity) interpolation, the Galerkin technique is employed. In this paper, DAMO is used to capture variations in saturation and pressure fields, leading to increased mesh resolution where required and coarsening elsewhere. Material properties such as permeability and porosity are constant and uniform so there is no need for interpolation. Figure 2 shows the flow diagram of the methodology described in Sections 2 and 3. Three major loops are considered. The dotted line denotes the fixed-point iteration (FPI) method [41], used to solve the non-linear system of equations. Next, the accumulated time is calculated. If the final time has not been reached, the time is increased and the algorithm may enter the loop in which the mesh is adapted. The DAMO loop is denoted with dashed line. Ultimately, the process is repeated in the time loop, denoted with solid line, until the final time-level is reached.

Numerical simulations
We now discuss the set-up of the immiscible displacement simulations that we consider in this paper. We start with a 2D problem that is computationally tractable on a fixed mesh. The wetting phase (such as water) is injected into a square geometry initially saturated with the nonwetting phase (such as oil) at the irreducible wetting phase saturation.
Displacement instabilities can be triggered numerically by introducing a perturbation to the saturation, pressure or permeability fields. Here, viscous fingering is triggered by a wetting phase saturation perturbation of Fig. 3 along the inlet boundary at the first time step. The perturbation can where L is the domain length transverse to flow, ξ = z − v shock t defines a coordinate system moving with the advancing linear shock front, p(ξ ) and s(ξ ) are pressure and saturation eigenfunctions respectively that have no explicit time dependence in the moving frame, σ = 2nπ/L is the wavenumber of the perturbation and ω is the growth rate of the perturbation. In order to be consistent, in the 3D cylindrical case, the same explicit saturation perturbation is used as a linear superposition along the diameter transverse to the flow. We begin by performing a mesh convergence analysis to determine the fixed mesh resolution needed for a physically converged solution in 2D. As the metric to measure convergence, we use the time taken for the wetting phase to reach the outlet boundary (breakthrough time).
Having found a converged fixed mesh fingering solution, we then repeat the process for simulations using DAMO. In the adaptive simulations, the initial mesh resolution is the same as in the respective fixed ones. The mesh is initially kept fixed so as to allow the fingers time to grow. Once the dimensionless time of pore volume injected (PVI)∼0.03, the mesh adapts to the water saturation field in every timestep. This delay ensures maximum resolution of the finger growth and that the early-time growth of the fingers is not influenced by mesh adaptivity. Adapting the mesh too early (before the finger pattern is established) or too late (after the fingers have progressed into the coarser region of the mesh) will result in finger growth being controlled by numerical diffusion rather than capillary pressure. In these The converged fixed mesh solution is then compared to the converged DAMO solution to determine how accurately the adaptive solution was able to reproduce the fixed mesh results. In order to fully utilise the capabilities of DAMO, we further investigate two alternative adaptive meshes with higher mesh anisotropy. In these tests, resolution is controlled by changing the maximum element edge length and the aspect ratio of the element size, whilst keeping the converged minimum element edge length. In this way, we manage to reduce mesh resolution where it is not required whilst maintaining the same level of accuracy.
We further assess the computational efficiency of the adaptive simulations compared to the respective fixed ones. We show that even when adapting and interpolating at every time step, the computational overhead is easily outweighed by the reduced computational cost resulting from the significantly smaller number of elements employed when implementing DAMO. More information about the computational cost of the DAMO can be found in [31,[42][43][44].
We move on to two examples that are much more challenging to simulate using fixed meshes and which rapidly become intractable at high resolution. First, we considered the above 2D displacement problem but in a much higher aspect ratio domain (10:1) which is closer to those encountered in real reservoirs. Finally, we considered a 3D simulation in which the wetting phase is injected into a cylindrical geometry saturated with the non-wetting phase and 3D fingers are allowed to form. This set-up is typically used in core flood experiments such as those of Riaz et al. [45]. These geometries are much harder to model in fixed mesh simulations using cuboid grid blocks.
The set of rock and fluids properties used for all simulations are shown in Table 1, based on the experimental data of Riaz et al. [45] and simulations of these experiments by Jaurè et al. [46]. The model dimensions for the three cases are shown in Table 2. A Corey correlation [47] is used for the relative permeability curves, which are the same as those used by Jaurè et al. [46] where S wr = 0.3 and S or = 0.4 are the immobile fractions of the displacing/displaced fluids (the fractions of the two fluids that cannot be displaced). At the specified viscosity ratio M = 303, the above displacements have a shock front mobility ratio M * = 1.8 (the theoretical stability limit is 1) and are expected to be unstable to fingering. For the capillary pressure, the functional form in [46] is adopted: In all 2D simulations, a P 1 DGP 2 CV (linear discontinuous velocity, quadratic pressure) element pair is used, whilst the 3D simulations use P 0 DGP 1 CV . In all cases, a small CFL ∼ 0.05 is used to minimise numerical diffusion. Table 3 summarises the different cases simulated and their associated mesh parameters.

Fixed mesh 2D model
We now discuss the results of the viscous fingering simulations outlined in Section 4. We start with a mesh refinement analysis to demonstrate numerical convergence and to determine the fixed mesh resolution required for a converged solution. The immiscible viscous fingering results are presented in Fig. 4 for different mesh resolutions, in the presence (Fig. 4a-d) or in the absence (Fig. 4e-h) of capillary pressure.    We observe that convergence under grid refinement is achieved only in the case of solutions with capillary pressure. In the absence of capillary pressure, it is not possible to achieve convergence, because refining the mesh without capillary pressure reduces numerical diffusion so fingers with higher wavenumber can grow with higher rates (Fig. 4e-h). Inclusion of physical diffusion (capillary pressure) regulates the problem, suppressing the growth of high frequency modes and leads to a converged solution (see, e.g., the linear stability analyses [15,48]). This is also verified quantitatively in Fig. 5, where the saturation profile along a vertical slice through the fingers at x = 0.0325 and a horizontal slice along the fingers at y = 0.025 are presented at time P V I ∼ 0.07, for different mesh resolutions. From Fig. 5 a and c, we can infer that in the case of solutions with capillary pressure the error is reduced as the number of elements is increased, and thus the viscous fingering pattern and the discontinuity in the saturation are accurately captured given the appropriate grid resolution.
In the no capillary pressure simulations, we observe a different flow pattern that is controlled by numerical diffusion. As the mesh resolution increases, the number of fingers increases from 6 for mesh #3 to 20 for mesh pressure. Non-monotonic curves of (d) are due to horizontal sampling across curved fingers where the tip is becoming disconnected from the fluid behind . Even for small changes in the grid refinement, the finger pattern changes significantly, as shown in Fig. 5b. This is also evident in Fig. 5d, where we observe sudden jumps in the saturation profile along the fingers due to the tips becoming disconnected. Finally, Fig. 6 shown (a, b, c, d), together with the corresponding dynamic mesh at these times (e, f, g, h). Notice that the mesh is initially held fixed and that it is the same as that used in the fixed mesh simulations mesh simulations stabilises at ∼ 40,000 elements, indicating convergence at this resolution (Table 3, mesh #7).
Having demonstrated convergence, we now look at the converged fingering solutions themselves. The converged fixed mesh simulations were performed on mesh #7 with ∼ 40, 000 elements. Figure 7 shows the time evolution of the viscous fingers in the presence (a-d) or in the absence (e-h) of capillary pressure. This is the first time DCVFEM has been used to simulate immiscible viscous fingering. The solutions with no capillary pressure (Fig. 6e-h) resemble those found in [24,49] using a control volume Galerkin formulation, indicating the correctness of the method. We note in particular that the flow pattern is dominated by many small fingers. This is in contrast to the case with capillary  (Fig. 6a-d) where we see fewer, thicker fingers. This behaviour is expected as capillary pressure diffuses nearby fingers, effectively joining them together.

DAMO 2D model
We now repeat the simulations with capillary pressure using DAMO to demonstrate that the same solution is obtained with significantly lower computational cost. Similar to the fixed mesh case, it is important to demonstrate convergence. Figure 6 (red line, diamonds) also shows a plot of breakthrough time versus mesh resolution for the DAMO simulations. The adaptive mesh simulations converge to the same value of the breakthrough time as predicted from the fixed mesh simulations. Here, it is possible to obtain converged solutions using DAMO with fewer than 5000 elements (more than 8 times reduction in the number of elements compared to the respective fix mesh simulation). This is because DAMO concentrates mesh resolution at the finger front, coarsening elsewhere, making the initial fixed mesh resolution irrelevant as far as computational cost is concerned. We now consider the converged adaptive solution with capillary pressure. This adaptive solution has the same minimum element edge length as that found for the converged fixed mesh (e = 0.0004). In Fig. 8a-d, the time evolution of viscous fingering using DAMO is presented together with the corresponding dynamic meshes (Fig. 8eh). The initial mesh (Fig. 8e) is the same as that used for the fixed mesh simulations and is held in place until PVI∼0.03 before being allowed to evolve. At that time, fingers have reached one-third of the way across the domain. After that time, the mesh adapts in every time step. Notice how mesh resolution automatically tracks the evolving finger front, putting resolution where it is most needed. The mesh is also refined near the inlet. This latter refinement is due to the rapid change in the water saturation at the trailing edge of the rarefaction behind the fingering (Fig. 6). The DAMO parameters of meshes #11 and #12 in Table 3 are now considered. In Fig. 6 (green line, cross), a plot of breakthrough time versus mesh resolution for adaptive simulations of various element aspect ratios is presented. We see that for all adaptive cases the solution stabilises at the converged value of the breakthrough time with less than 0.1% variability. Remarkably, a converged solution is obtained without loss of accuracy, with fewer than ∼2600 elements. This leads to a 16 times reduction in the number of elements compared to the fixed mesh simulation. Figure 9 shows a close-up of the adaptive mesh for various element aspect ratios.
There is no visual difference between the adaptive results and the fixed mesh results (Fig. 7) suggesting that the adaptive mesh is sufficient to capture the key flow behaviour. This is verified quantitatively in Fig. 10a, where we plot the saturation profile across a vertical slice through the fingers at x = 0.0325 at time P V I ∼ 0.07, from both adaptive and fixed mesh simulations. In Fig. 10b, we compare the wetting phase flux across the outlet boundary for the respective simulations. In all cases, there is very close agreement between the fixed and adaptive results indicating that adaptive mesh simulations are capable of accurately capturing the immiscible viscous fingering behaviour whilst significantly reducing the computational cost.

Computational speed-up using DAMO
The significant reduction in the number of elements needed for convergence leads to a major computational speed-up on the adaptive mesh compared to its fixed counterpart. Figure 11 shows the normalised CPU time against the minimum element edge length for the fixed meshes and the adaptive ones with fixed maximum edge-length (meshes #1-7). The CPU time is normalised by the time taken to run the converged fixed mesh simulation (mesh #7). We see that the converged adaptive mesh simulation runs 3.5 times faster than the fixed one. The normalised CPU time for the meshes #11 and #12 are also plotted in the same figure.
The higher mesh anisotropy leads to further reduction in the number of elements required. Consistently, lower computational effort is demanded and speed-up factors of 6.8 and 10.9 are achieved, respectively. We see that the additional computational cost for adapting the mesh is insignificant compared to the reduced cost due to the  In all cases a close-up is presented to illustrate how mesh adaptivity automatically refine at the finger front

High aspect ratio 2D model
Whilst DAMO is very useful in the above simulations, it is still feasible to use fixed mesh simulations due to the relatively small domain size. We now consider an example that is much more challenging on a fixed grid. Since DAMO focusses resolution mainly at the front, it is ideal in cases where one needs to simulate viscous fingering in high aspect ratio domains. This is needed to probe the late-time behaviour of fingering where the transverse dimension is much smaller than the direction along the flow (typical of laboratory core floods, oil reservoirs and aquifers used for geological CO 2 sequestration).
As a proof of principle, we show in Fig. 12a-f viscous fingers in a 2D domain with L = 0.5m, giving an aspect ratio of 10. We see that the number of fingers in the system has reduced from 8 at early time (P V I ∼ 0.01) to 4 just before breakthrough (P V I ∼ 0.1). This is due to fading of the smaller fingers. This fading and merging of fingers will continue and it is likely that the fingers will eventually join to form a single finger (see [50]). In this case, it seems likely that the overall behaviour of the displacement will be independent of the initial conditions.
These results represent a step towards a systematic study of non-linear late-time immiscible viscous fingering. Fixed mesh simulations in these geometries are extremely costly but this is much less so in the adaptive case. The equivalent fixed mesh resolution would require approximately 415,000 elements whilst the adaptive one uses on average around 22,000 elements leading to a 19 times reduction in the required number of elements. In Fig. 11, the normalised CPU time against the minimum element edge length is Fig. 15 Cylindrical viscous fingering pattern shortly before wetting phase breakthrough with the adaptive mesh overlaid presented. We show that a speed-up factors of approximately 18 is achieved.
In Fig. 13, the adaptive mesh resolution is plotted over the normalised simulation time. We observe a constant number of elements for the initial time steps where (similar to above) the mesh is held in place to ensure proper finger growth. We notice a sudden drop of the required number of elements since DAMO reduces them by half, only using smaller elements where they are needed to track the evolving finger Solid lines represents curve fitting of the raw data presented with a faint gray lines, together with the respective error bars in black lines. In both cases adaptive mesh converge to ∼ 1.1 million elements front and the rapid change in saturation at the inlet. At breakthrough, the mesh resolution is further reduced since viscous fingers cross the outlet boundary and there is less variation in water saturation in the system itself.

Cylindrical 3D model
As a final demonstration of the advantages of DAMO for immiscible viscous fingering simulation, Fig. 14 shows the results of simulating viscous fingers in a 3D cylinder. This demonstrates the practical utility of our approach.
In Fig. 15, the cylindrical viscous fingering pattern is presented, shortly before wetting phase breakthrough overlaid with the associated dynamic adaptive mesh. The adaptive mesh has an average number of ∼ 1.1 million tetrahedra and was run in parallel on 20 cores using a Linux workstation with 2.6 GHz Intel Xeon processors. The equivalent fixed mesh resolution would have required approximately 40 million elements rendering this problem prohibitively computationally expensive. Utilising DAMO reduces the required number of elements by a factor of 36, making this problem feasible to solve. A computational efficiency of the same factor is estimated. In Fig. 16, the adaptive mesh resolution is plotted over the normalised simulation time for two different adaptivity settings. Similar to the 2D case, we observe a sudden drop in the mesh (a) (b) Fig. 17 a Viscous fingers and b the computational mesh before breakthrough as a longitudinal slice at the centre of the cylinder, through the fingering pattern. The finger growth is dominated by non-linear effects and mesh adaptivity continues to automatically refine in the regions where the viscous fingers are growing resolution once the initial mesh is allowed to adapt and a quick convergence of the maximum required resolution.
The cylindrical fingering pattern resembles that of the 2D case. Figure 17 shows the saturation pattern and the respective mesh seen in a longitudinal slice through the viscous fingers (plane normal to the x axis and passing through the centre of the cylinder). We observe that, as in the 2D case, the number of fingers in the system has reduced to 4 and that mesh adaptivity continues to automatically refine in the regions where the viscous fingers are growing. This similarity is due to the fact that we used the same saturation perturbation at the inlet boundary to control finger growth. Figure 18 shows the saturation profile across cross-sectional slices for various times at different distances along the cylinder. From this figure, the nonlinear late-time behaviour of the fingers is also verified in the cylindrical set-up, where small fingers tend to fade and/or merge together forming thicker ones that continue to grow. A similar response is observed experimentally in [45]. Further investigation is required to determine the physical

Conclusions
In this paper, we have demonstrated the benefits of using DAMO for the modelling of immiscible viscous fingering. The mesh adaptivity was implemented in a DCVFEM model using an unstructured 2D and 3D tetrahedral mesh. Whilst adaptivity techniques are now ubiquitous in other applications of CFD, their use in porous media flow remains limited. Conventional simulation of viscous fingering requires a combination of very fine mesh resolution and higher order discretisation to ensure physical diffusion dominates numerical diffusion.
This paper shows that DAMO is ideally suited to viscous fingering applications, that would be prohibitively expensive on fixed grids. Using DAMO fine resolution is only needed along the fluid interface where the fingers are growing and not throughout the domain. We showed that using DAMO in a small 2D case, leads to a speedup of factor of ten, with no loss in accuracy. The compositional cost was decreased accordingly, resulting from the significantly smaller number of required elements.
Two challenging examples were also discussed. A 2D simulation of the late-time non-linear regime of fingering in a high aspect ratio domain and a 3D simulation of viscous fingering in a cylindrical geometry. These examples showed that the DAMO can be used to begin a systematic study of non-linear fingering which has not been attempted thus far in the literature and could yield important insights into the scaling behaviour of viscous fingers at late times. The 3D example is of particular interest as it is a step towards simulating viscous fingering in more realistic, non-trivial geometries that more closely match experiments. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.