Incomplete mixing in porous media: Todd-Longstaff upscaling approach versus a dynamic local grid refinement method

Field-scale simulation of flow in porous media in presence of incomplete mixing demands for high-resolution computational grids, much beyond the scope of state-of-the-art simulators. Hence, the upscaling-based Todd and Longstaff (TL) approach is typically used, where coarse grid cells are employed with effective mixing fluid properties and parameters found by matching results obtained with fully resolved reference simulations. Dynamic local grid refinement (DLGR) techniques, on the other hand, only employ fine-scale grid resolution where the fully mixed assumption is not valid. The rest of the domain is then solved at coarser resolutions, where the fully mixed assumption is valid. Here, we assess the accuracy and the robustness of DLGR- and TL-based simulations of miscible displacements in homogeneous and heterogeneous porous media. Due to the intrinsic uncertainty within the unstable displacement nature of the studied incomplete mixing processes, the performance of the methods is also investigated based on a range of acceptable solutions rather than relying only on a single reference one. Systematic numerical results illustrate that the DLGR method is much more robust and accurate than the upscaling-based TL approach, and employs only a small fraction of fine-scale reference grids. Especially, the TL upscaling results (though history matched with computationally expensive fine-scale results) are very sensitive to the change of the simulation parameters. Based on this study, we propose a dynamic multilevel simulation strategy for efficient and reliable large-scale simulation of the complex incomplete mixing processes.


Introduction
The miscible and immiscible displacements in subsurface porous media can develop instabilities at the interfaces of the fluids with different properties [1]. Moreover, reservoir rock heterogeneity plays a dominant effect on how fluids are displaced through the porous reservoir rock. Depending on its values, variance, and spatial correlation, the permeability distribution can either accentuate or suppress the viscous fingering phenomenon [2][3][4][5][6][7][8]. To allow for fast decisions Matteo Cusini M.Cusini@tudelft.nl Extended author information available on the last page of the article. with the available computational capacities, the state-ofthe-art reservoir simulations are typically performed on far lower grid resolutions than the scale at which the incomplete mixing, i.e., fingering and unstable interfaces, occur. Apart from introducing additional numerical diffusion, these lowresolution grid cells may suppress small-scale physical phenomena that contribute to unstable displacements, leading to inaccurate simulations [9]. There has been a significant effort to account for these sub-grid features through effective upscaled fluid properties [10][11][12]. Among them, Todd and Longstaff's (TL) model [11] is widely used in the petroleum industry.
Similar to alternative upscaling-based models, a major disadvantage of TL approach is that the effective fluid properties are described by an undefined scalar mixing parameter, which must be obtained from either highresolution simulation or experimental data. The use of full-field high-resolution simulations, in order to find these tuning (effective) parameters, is impractical. Moreover, the tuned parameters obtained for a given simulation settings may not be valid if the simulation inputs are changed. As such, the applicability of these approaches to properly address incomplete mixing processes is questionable. Of particular interest is to investigate the applicability of the TL approach, and, if it is found impractical, to develop an alternative simulation approach applicable for field-scale studies.
In this work, the accuracy and robustness of DLGR and the widely used TL model for the simulation of incomplete mixing are investigated. The sensitivity of DLGR simulations to different error estimate strategies in order to refine or coarsen the grid is also assessed. Both DLGR and TL are implemented in a commercial-grade simulator for quality benchmarking purposes. As the studied fluid flow involves unstable displacements, the variations of the fine-scale reference solution based on different perturbations are used to define "acceptable solution" range rather than a specific fine-scale solution map.
Numerical experiments show that DLGR is capable of simulating incomplete mixing phenomena with only a fraction of the fine-scale grid blocks. Important is that although TL upscaling approach was employed after tuning parameters based on the fine-scale reference solution, its results are sensitive to changes of the simulation inputs (e.g., grid resolution and mobility ratio). The DLGR strategy, on the other hand, is consistent when simulation parameters are changed. Based on the results presented here, DLGR should be preferred over the TL model, since it allows to accurately and efficiently simulate the incomplete mixing displacement, and it shares the same sensitivities with respect to the input parameters compared with the fine-scale fully resolved simulation.
The manuscript is organized as follows. The fine-scale governing equations are presented in Section 2. Then, the DLGR simulation strategy is briefly described in Section 3, with details of several error criteria approaches that have been considered in this work. The formulation employed by the upscaled TL model is described in Section 4. To investigate the accuracy and consistency of DLGR and TL models, numerical test cases are presented in Section 5 for both homogeneous (5.1) and heterogeneous (5.2) media. Concluding remarks are presented in Section 6.

Fine-scale model equations and solution strategy
The equations describing a single-phase incompressible system of two miscible components (oil and a solvent) flowing in a porous medium [27] are and φ ∂c ∂t where φ is the porosity of the porous medium, t is the time, c is the volume fraction (or concentration) of the solvent, and q t and q c denote, respectively, the total and the tracer source term with [1/s] dimension (only through injection and production wells). Note that since the source term is only induced by the wells, the maximum concentration value cannot exceed 100% void volume of a grid block [28]. Also note that at c = 0, the convective term and concentration source in Eq. 2 become zero, and at c = 1, q c = q t . Thus Eq. 1 enforces the net outgoing fluxes in Eq. 2 to balance the tracer source term. Furthermore, the Darcy velocity u in the absence of gravitational effects reads where μ mix is the viscosity of the mixture, K is the permeability tensor, and p is the pressure. Equation 3 is substituted in Eq. 1 to obtain a pressure equation, i.e., Here, μ mix is calculated based on a quarter-power mixing rule [10], i.e., Finally, D in Eq. 2 is the dispersion tensor representing molecular diffusion and mechanical dispersion [29]. In a 2D domain, the dispersion tensor reads where Here, u x and u y are the x and y component of the Thus, the longitudinal (Pe l ) and transverse (Pe t ) Péclet numbers are defined as and where L x and L y are the dimensions of the domain along the x and y axes respectively and A is the cross section. Finally, q inj is the volumetric flux of the injected solvent. An unstable front (i.e., viscous fingering) is generated throughout the miscible displacement process whenever the injected solvent reduces the viscosity of the fluid phase (thus μ s < μ o ) as shown by the linear stability analysis [30][31][32].
The system of governing equations (i.e., Eqs. 2-4) is discretized with a finite-volume method and solved using a sequential approach. First, Eq. 4 is solved to obtain the pressure field that is used to compute velocity by using Eq. 3. Then, Eq. 2 is solved using an operator splitting technique. The linear advection term is solved using an explicit time discretization whereas the dispersion term is solved implicitly due to the stability considerations. A second-order (two-point upstream) scheme along with high-resolution time step sizes is employed to limit the effect of discretization error that may interfere with the physical dispersion.

DLGR method
A nested DLGR technique with implicit evaluation of the grid resolution is employed [17] to solve (2)-(4). At the beginning of each time step n, the solution grid of timestep n − 1 is locally coarsened or refined based on a pair of user-defined refinement and coarsening criteria that aim at employing high grid resolution only at the advancing front location. Once the solution grid has been chosen, the non-linear system of equations is solved. If the converged solution violates the refinement criteria, the grid is updated and a new solution is found. This process is referred to as repeat time step and it is shown in Fig. 1. Previous studies have shown that the overhead inherent to the inclusion of the additional solution step is minimal [17]. Here, the maximum resolution jump between two neighboring cells is forced to be at most equal to one level of grid refinement to limit the discretization error [13,33].
A variety of front-tracking criteria have been presented in the literature both involving fluid properties [13-15, 34, 34-36] and mass or volume fluxes [15,16,37]. Here, given a variable x (e.g., saturation, concentration, and viscosity), the following quantities can be considered to determine the appropriate grid resolutions: For any given measure (m), two threshold values, δ r and δ c , are defined by the user to determine, respectively, when to refine (if m > δ r ) and coarsen (if m < δ c ) the grid. Note that the threshold values for refinement and coarsening can be set independently as long as δ c ≤ δ r . In this work concentration, mixture mobility and mixture viscosity have been considered as properties to locate the advancing solvent front. Whenever a grid cell is refined, piece-wise constant functions are employed to interpolate all variables (i.e., pressure and concentration) to the finer resolution. On the other hand, when coarsening is performed, the average values are assigned to the coarse blocks.
In heterogeneous domains, the finest level (or maximum) corresponds to the resolution of the given permeability field. Coarser permeabilities are obtained with a local flow-based upscaling procedure (Fig. 2) [38]. Here, the focus of the studies is on comparing TL and the DLGR approach.

Upscaled TL model equations and solution strategy
Our main motivation to revisit the TL model in this work is that it has been widely used in the literature to investigate the incomplete mixing process. As such, a complete investigation of any advanced method, including DLGR, would require a comparison with the outcome of this model.
The TL model [11] was originally introduced to represent mixing phenomena when a solvent is added to a water-oil system with the use of three-phase black oil equations [28]. The main advantage of this approach was that it allowed for the use of a black oil simulator (instead of a compositional one) only by adjusting the fluid properties. Furthermore, the TL model aims to represent the incompletely mixed fluid displacement at a much coarser grid resolution than the mixing (fingering) resolution. Consequently, the effect of mixing is captured by introducing the effective mixture coefficients. These coefficients are obtained, typically, by tuning the TL solutions against a few expensive training fine-scale simulations. Therefore, the TL model does not represent the details of the displacement physics (e.g., fingering front shape and velocity), despite being widely used to estimate the production history [39].
The equations of the TL-based approach for miscible oil (o) and solvent (s) flow with no gravitational effects read Here, S α , k r,α , and q α are saturation (volume concentration of oil and solvent), source term (wells fluxes), and phase relative permeability, respectively. Because the two components (oil and solvent) are assumed miscible, there exist no interfacial forces between them, and thus linear relative permeabilities are employed in the TL model [11]. Finally, μ α, eff is the effective viscosity computed as The mixture viscosity μ mix is calculated from Eq. (5) where c = S s . Here, ω is the TL-mixing parameter and its value determines the mixture effective properties. In particular, ω = 1 corresponds to the case of complete mixing and ω = 0 corresponds to the immiscible one. Todd and Longstaff reported that the use of ω = 2 3 gave satisfactory results for the simulation of miscible displacements with an oil/solvent viscosity ratio of 86 [11]. However, it should be noted that the value of the empirical mixing parameter, ω, highly depends on changes in properties such as grid block size, permeability, geometry, and fluid properties, and therefore, cannot be easily generalized. For example, Fig. 3 shows two possible concentration distribution scenarios within a simulation grid block. The left figure shows a stable dispersed front and the right one shows an unstable displacement with both frontal dispersion and viscous fingering. The ratio between the size of the grid block and the thickness of the mixing zone determines the value of ω that should be used. If the grid block size is small compared to the dispersed zone, the solvent and oil fluid parameters can be seen as completely mixed and ω = 1 can be considered. On the other hand, if frontal dispersion is negligible compared to the grid block size, pure component properties can be employed; thus, ω = 0 [11].
The constraint S o + S s = 1 can be used to eliminate one unknown so that Eq. 12 leads to a well-posed system of equation for two unknowns (i.e., p and S s ). Here, Eq. 12 is solved employing the finite-volume method on structured Cartesian grids and an implicit time discretization. A twopoint flux approximation is used to evaluate interface fluxes along with a first-order upwind for all fluid properties.

Case 1: homogeneous porous media
A 100 m × 100 m reservoir, with a permeability of 400 mD and uniform porosity equal to 0.3, is considered. No flow  [38] boundary conditions are considered at all boundaries and an injection well perforates all cells at the left boundaries whereas a producer is located in all cells at the right boundary. The solvent is injected at the fixed rate of 1 ft day (in the left boundary). At the right boundary, i.e., production cells, the fixed bottom hole pressure of 350 bar is imposed. The viscosities of the single components are μ s = 1 cP, for the solvent, and μ o = 100 cP, for the oil. The reservoir is initially saturated with oil and simulations are run until breakthrough occurs with time step sizes corresponding to a CFL number equal to 1. Instabilities are generated by randomly perturbing the saturation distribution in the cells perforated by the injection well (see the Appendix). Three different physical scenarios are analysed: a zero, a medium, and a high dispersion case, resulting in different longitudinal and transverse Péclet numbers as presented in Table 1.

Convergence of fine-scale simulations
The solvent injection process, for all three physical scenarios, is solved with different fine-scale grid resolutions. Starting from a 25 × 25 grid, 2 × 2 refinement is applied up to obtaining a 400 × 400 grid. Figure 4 shows the grid sensitivity of the solvent concentration profile at 0.2 pore volume injected (PVI) upon grid refinement. The left column of Fig. 4 corresponds to the solutions obtained without any physical dispersion, so that discretization errors are the only cause of smearing out the solution. Note that, since discretization error reduces with grid refinement, the solution does not converge. The second and third columns show reducing instabilities at the front, caused by the increased physical dispersion. As a result of the physical dispersion, small perturbations are smeared out, allowing only large fingers to grow. It can be observed in Fig. 4 that the two highest resolution simulations show a converging trend, which indicates that the physical dispersion dominates numerical diffusion.

Accuracy of the DLGR method
In this subsection, the accuracy of the DLGR method and its sensitivity to the choice of the refinement and coarsening criteria are evaluated for cases 1a, 1b, and 1c. The highresolution (400 × 400 grid) simulations discussed in the previous subsection are employed as the reference. In order to mimic the physical uncertainty, a range of possible results is obtained by generating ten statistically equivalent initial tracer distribution at the injector. The standard deviations of the oil production and solvent rates are presented in Table 2. This allows for quantification of the underlying uncertainty,  so to allow for an acceptable solution range for DLGR and TL simulations. The error of DLGR simulations is quantified using the following error indicators: -Relative error of the cumulative oil production (cp): -Relative error of solvent rate (SR) at the producer well -Second norm of the concentration difference: Here, the subscript f s stands for fine-scale (reference) and NOE is the number of grid blocks with a solvent concentration (c) higher than 0.01.
For the DLGR simulation, the grid resolution has been selected based on the values of the following quantities:  1a, 1b, 1c) are performed, each using the tolerance shown in Table 3. Employing the value of the time derivative of the concentration gradient ( c 3 ) as an indicator of the front location has shown promising results in the literature [17]. However, in the case of viscous fingering, this will lead to very early coarsening of the rear part of the fingers which is not subject to large changes in time. This affects the development of the instabilities. Therefore, the criterion c 3 > δ r is used for refinement whereas c 1 < δ c is used to determine where to coarsen. Figure 5a-c shows the DLGR errors for different refinement criteria as functions of the average number of the employed grid blocks (indicated as a percentage of the number of cells of the fine-scale grid). The red line indicates the allowable mean deviation from the reference solution, due to the variations within the fine-scale results for different perturbations. The average number of repeated time steps (for the implicit dynamic grid strategy) is shown in Fig. 6. Note that the DLGR performance is similar Fig. 7 Normalized gradients of phase viscosity (blue), concentration (magenta), and phase mobility (green) against dimensionless length. The corresponding concentration gradient is shown in red for difference-and gradient-based operators, but different between the quantities of μ, c, and λ.
In fact, viscosity-based criteria (shown in blue) perform slightly better than the mobility (in green) and the concentration (in magenta) ones. A possible explanation can be deduced when analyzing the normalized gradients illustrated in Fig. 7.
As described earlier by Eq. 5, viscosity is a strong function of the concentration. Figure 7 illustrates that the gradient of viscosity criterion tends to refine the grid relatively towards the tip of the front, compared with the gradients of mobility and concentration strategy. As such, the front is captured more accurately, when a viscosity gradient (or difference) is being used. Figure 8 illustrates this concept by showing the concentration solution and the grid at 0.2 PVI for the normalized viscosity and mobility gradient criteria using the tolerance of 0.05. It is observed that, as mentioned above, the viscosity criteria refine the grid more towards the front tips where the resolution is actually needed. It is also clear from Fig. 6 that the viscosity-based approach requires the minimum re-adjustment of the implicit dynamic grid strategy. It should be also stressed that these results obtained for the phase viscosity and mobility criteria are based on the quarter-power mixing formulation, as mentioned in Section 2. Therefore, if other functions are used, one has to first study which criteria tend to employ the fine grids towards the tip of the front.
Due to the unstable nature of the front, different grid resolutions impose different perturbations, which in the absence of physical dispersion can cause different fingering patterns. As such, the concentration error c , shown in Fig. 5, does not show any specific trend when more refined grids are being employed.
Both the space-time derivative and the third-order spatial derivative of concentration were able to perform within the allowable range of cumulative oil production and solvent rate of 1% and 11% respectively. Moreover, the mean repeat time steps needed per simulation time step strongly increase for both error estimators when the tolerance value is relaxed (as high as 5.3 repeat time steps per simulation time step shown in Fig. 6). Figures 9 and 11 show the accuracy evaluation for the medium and high dispersion cases, respectively. The acceptable error ranges are shown by the red lines, which indicate the uncertainty within the fine-scale solution. Figures 10 and 12 show, instead, the number of repeat time-steps for these two scenarios.
For the medium dispersion case, based on the cumulative oil production, acceptable results are obtained by employing 15% of the fine-scale grid cells. For the high dispersion case, all the studied values fall below the allowable error line; thus, only a few percent active grid cells are enough. The concentration error (mean l 2 − norm) and the number of repeat time steps, illustrated in Figs. 9c, 10, 11c and 12, show that DLGR is more accurate when using difference, i.e., Δ, based criteria. This can be explained by the different asymptotic behaviors, i.e., Here, f is an arbitrary smooth function of location x and Δx is the discretization step. Additionally, Δf = f (x + Δx) − f (x). When Δx tends towards zero, the ratio Δf Δx converges to the first derivative of f with respect to x (i.e., ∂f ∂x ) whereas Δf tends to zero. This means that the gradient criteria require an enforced maximum refinement level in order to limit the amount of refinements, because the true derivative, ∂f ∂x , may be larger than the refinement tolerance (so the refinement never meets the tolerance value). On the other hand, the gradients can be lowered by simply coarsening the grid, which then leads to smearing out the fronts. This is illustrated in Fig. 13 showing the solution for concentration difference (top) and gradient (bottom) criterion at 0.6 PVI using the normalized tolerance of 0.2. Figure 14 shows the amount of grid blocks used in each refinement level during the simulation time for both the concentration difference and gradient criteria, confirming that more intermediate refinement levels are used with the difference criterion. The distinct difference between gradient-based and difference-based criteria in the mean repeat time steps used during during the simulation, as shown in Figs. 10 and 12, emphasizes the added benefit of the difference-based strategy to refine the grid.
It is observed that the error criteria using the space-time derivative of the spatial gradient of the solvent concentration performed well for simulations including dispersion. This is reflected in the error in production data and the concentration solution, as it can be seen in Figs. 9 and 11. However, as the tolerances are loosened, the error in cumulative production strongly increases, suggesting that the space-time derivative criteria are smaller than the tolerances, and thus no local grid refinements are introduced.

Upscaled TL model simulations
In this subsection, the sensitivity of the upscaled TL model to the grid block size and the parameter ω are investigated for case 1b (medium dispersion). Figure 15 shows a surface plot of the error made by the TL model in the cumulative oil production as a function of the grid block size and the mixing parameter. At three different locations on the surface plot, the actual difference in cumulative production is visualized in order to relate the error.  Fig. 15, it can be observed that the range of accurate TL results (i.e., results within the "allowable deviation" of 1.2%) obtained for a given grid resolution are narrow.

TL model vs. DLGR method: sensitivity to domain size
Here, four homogeneous reservoirs, of equal heights (25 m) and lengths of 200, 100, 50, and 25 m, are considered. The same injection rate (1 ft/day) and fluid properties are used in all reservoirs and no dispersion is present. The TLmixing parameter is fit to match the fine-scale simulation of the largest reservoir. The cumulative oil production curves obtained with fine-scale, DLGR, and TL simulations for all four reservoirs are shown in Fig. 16. In each figure, the high-resolution simulation curves are depicted in red, DLGR in green, and TL effective model in blue. The top left (a) corresponds to a reservoir length of 200 m. DLGR uses the concentration difference across grid block interfaces as error criterion, using a threshold value of 0.05. On average, this resulted in a mean active grid block count of 31% of the fine-scale simulation grid block count. It is clearly observed that DLGR is capable of capturing the physics independent of the scale of the system. It is also observed that the number of grid blocks used by DLGR reduces when increasing the size of the domain. On the other hand, the TL model is not able to scale with the physical solution. The mixing parameter used shows an underprediction of local mixing, resulting in incrementally earlier solvent breakthrough for the downscaled reservoir sizes.

TL model vs. DLGR method: sensitivity to mobility contrast
The same 100 m × 100 m reservoir previously described is considered along with the three different physical scenarios (cases 1a, 1b, 1c). The injection process is simulated with fine-scale (400 × 400), DLGR, and TL model (40 × 40), for three different mobility ratios M = μ o μ s : M=10, M=100, and M=500. DLGR simulations use the concentration difference error criterion and the mixing parameter of TL is fit to match the fine-scale simulation cumulative oil production curve for M=10.  18 and 19 show the cumulative oil production curves obtained with the three simulation strategies for all scenarios. For all three cases, DLGR provides accurate cumulative oil production curves irrespective of the mobility ratio. Note that the mean active grid blocks employed by DLGR tend to increase with the mobility ratio. On the other hand, the TL model only provides accurate solutions for the mobility ratio for which the parameter is tuned (here, M = 10). Figure 20 illustrates the matched (optimum) mixing parameter values against the mobility ratio for the three different scenarios. Clearly, a linear relation exists on a semi-log scale, allowing for potential rescaling of the mixing parameter in order to fit the data consistently. The literature reported a relation of the mixing parameter as a function of the mobility ratio [40,41]. This was obtained by equating the fractional flow formulation derived from the TL effective (upscaling) model with the fractional flow obtained from Koval's effective (upscaling) model [10], assuming an average solvent concentration at the front. This relationship, depicted with yellow markers, shows a good fit with the zero dispersion case, however, is incapable of predicting an optimal mixing parameter for cases with added dispersion.

Case 2: heterogeneous permeability with small correlation length
A 100 m × 20 m 2D heterogeneous reservoir is considered with no flow boundary conditions at all boundaries. A rateconstrained injection and a production well perforate all the cells at the left and right boundaries, respectively. The solvent is injected with a constant rate of 1 ft/day. The solvent and oil viscosities are 1 and 10 cP respectively. Three permeability fields with an average value of 90 mD, a variation corresponding to a Dyskstra-Parson (V dp ) coefficient of 0.63, and a correlation length, λ D =0.01, are used, one of which is shown in Fig. 21. Table 4 illustrates all simulation parameters. Note that, in DLGR simulations, the highest refinement level (level 4) corresponds to the fine-scale permeability map. Permeability fields of coarser resolutions are obtained with a flow-based upscaling procedure.
The concentration solutions for all three simulation techniques as well as the grid generated by DLGR, for one permeability realization, at 0.4 PVI are shown in Fig. 22. DLGR employs a concentration difference tolerance of 0.05 as refinement criterion. Note that DLGR is capable of capturing all small-scale flow features using, on average, 35% of the grid blocks compared to the fine-scale simulation (see Fig. 22b, c).
The calibration of the TL model-mixing parameter for heterogeneous systems obtained by equating the fractional flow formulation of the Koval and TL models is the only   [40]. However, this leads to the estimation of an erroneous negative value (-0.043) for the mixing parameter. By tuning the mixing parameter to fit the data of the high-resolution reference simulation, an optimal mixing parameter of 0.6 is found. Figure 22d shows the concentration profile obtained with the empirical TL model. Figures 23 and 24 show the cumulative oil production and the oil rate. In Fig. 23, DLGR employs a refinement Fig. 20 Case 1: mixing parameter match against the logarithm of mobility for zero, medium, and high dispersed cases. The results from the mixing parameter relation [41,42] are shown in yellow, showing good results for zero dispersion cases criterion tolerance of 0.05, resulting in an average grid block count of 37% of the fine-scale grid blocks. Figure 24, instead, shows results with DLGR using a tolerance value of 0.15, resulting in using only in 19% of the fine-scale grid blocks on average. It can be observed that the TL model shows low sensitivity to the variation in production and local solvent concentration distributions compared to the high-resolution and DLGR simulations.

Case 3: heterogeneous permeability with long correlation length
The same reservoir as the one of the previous case is considered along with the same fluid properties and the  same boundary conditions. Here, a permeability field with a higher correlation length is employed as shown in Fig. 25.
As for the previous case, the fine-scale permeability map is upscaled (locally) for TL and the coarser grid resolutions used by DLGR. Figure 26 shows the solvent concentration map at 0.4 PVI for fine-scale, DLGR, and TL simulations. Here, DLGR employs, on average, 29% of the fine-scale grid cells. The TL model is run using the optimal mixing parameter ω = 0.375. Note that the long correlation length reduces the local mixing as shown by the optimal value of the TL-mixing parameter. As for the previous case, DLGR is capable of capturing the physics of fluid flow altered by the underlying permeability field. However, due to the increased correlation length of permeability, the frontal surface area is higher than that in the previously described case resulting in a higher amount of refinements during simulation. Figures 27 and 28 show the cumulative oil production and the oil rate. In Fig. 27, DLGR uses a refinement criterion tolerance of 0.05 and in Fig. 28, a refinement criterion tolerance equal to 0.15. Similar to the previous case, DLGR results are in good agreement with the finescale reference one for all permeability realizations whereas the TL ones show a considerable error for some of them.

Conclusion
In this work, the accuracy and sensitivity of a DLGR method were compared against those of the TL upscaling approach for incomplete mixing process.
Numerical results show that DLGR has a great potential to overcome the limitations of effective models like those of the TL one, both for homogeneous and heterogeneous media. In fact, it was shown that, for a homogeneous domain, the TL model is very sensitive to modifications of the simulation parameters (both size of the domain and mobility contrast). Consequently, TL simulations require adequate tuning of the mixing parameter ω every time the simulation input changes. In many cases, the correct value of the TL-mixing parameter can only be found by running high-resolution simulations and an a priori estimate cannot be found. On the other hand, DLGR provides accurate solutions for simulations of incomplete mixing processes by employing only a fraction of the grid cells used for fine-scale fully resolved simulations.
Accuracy and robustness of the described DLGR strategy and the TL model for heterogeneous media with different correlation lengths were also assessed. DLGR  allows for local preservation of the underlying geological features, clearly indicating the potential of DLGR to solve for incomplete mixing in heterogeneous media. This is illustrated by the fact that the production results of DLGR show the variations similarly observed in the high-resolution simulations that were used as a reference. Furthermore, it was shown that the error of DLGR simulations can be controlled by regulating the tolerances used for the refinement criteria.
As such, DLGR strategies represent a promising approach to obtain more accurate results and overcome the limitations of upscaled models for simulations of incomplete mixing processes. Note that efficient implementation of the DLGR method is an important factor for its successful implementation in field-scale commercial simulators. Additionally, the use of upscaled properties (i.e., permeabilities) in coarse grid blocks may negatively affect the accuracy of the dynamic multilevel simulations, especially in the presence of more complex physics. These issues have been resolved by the algebraic multiscalebased ADM method [43]. As such, the developed DLGR approach of this article can be implemented in a commercial-grade simulator with an ADM-based implementation strategy, where constant interpolators are used for both pressure and concentration across different scales. Such implementation would allow for computational speedup measurements, which is the subject of our future work.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http:// creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Appendix: Concentration perturbations at initial time step
The solvent (randomly generated) perturbations (for the initial concentration values) imposed in the first column of the domain are shown in Table 5. To perturb the initial concentration values at the finer grid resolution N y , these values are applied per patches of N y /25 cells. Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.