A constrained proper orthogonal decomposition model for upscaling permeability

Reservoir modeling and simulation are vital components of modern reservoir management processes. Despite the advances in computing power and the advent of smart technologies, the implementation of model-based operational/control strategies has been limited by the inherent complexity of reservoir models. Thus, reduce order models that not only reduce the cost of the implementation but also provide geological consistent prediction are essential. This article introduces reduced-order models based on the proper orthogonal decomposition (POD) coupled with linear interpolation for upscaling. First, using POD-based models, low rank approximate (LRA) are obtained by projecting the high dimensional permeability dataset to a low dimensional subspace spanned by its trajectories to decorrelate the dataset. Next, the LRA is integrated into the interpolation algorithm to generate upscaled values. This technique is highly scalable since low-rank approximations are achieved by the variation in the number of modes used for reconstruction. To test the validity and reliability of the model, we show its application to the practical dataset from SPE10 benchmark case2. From statistics of the error analysis, the classical POD algorithm seems to be more preferred for LRA; however, since non-negativity of the permeability data set is a priority, the constrained POD (non-negative POD) algorithm described in this article is more appropriate. Finally, we compared the POD-based models to a traditional industry-standard upscaling technique (e

simulations remain challenging due to the high-resolution heterogeneous properties of the reservoirs. 1 Reservoir models described on fine-scale discretization are required to capture the rock heterogeneity.These fine-scale reservoir models result in nonlinear high-order systems of equations known as the full-order model.Performing simulations on the full-order model (FOM) is expensive, as vast computational resources and time are often required. 2Therefore, it is of great interest to reparameterize reservoir properties like permeability and porosity through reduced-order models.To this end, having a practical approach to represent the FOM of high dimensionality with a reduce order model (ROM) of coarse-scale and low dimensionality is important.
Generally, the main objective of reduced-order modeling techniques is replacing the FOM with a reduced order model (ROM) while preserving essential features of the original system, such as input-output behavior, stability, and passivity.These reduce order models are simplifications of the FOM, having a much lower dimension on a coarser scale. 3Among the ROM techniques, the grid-based upscaling method is popular for reparameterizing.Upscaling refers to the computation of effective coarse properties based on fine-scale features.In upscaling, a heterogeneous property region consisting of fine-grid discretizations is replaced with an equivalent homogeneous region, made up of a single coarse-grid cell with an effective property. 1pscaling techniques for subsurface reservoirs can be classified in terms of the parameters to be upscaled: single or two-phase flow parameters. 4In single-phase parameter upscaling techniques, the fine-scale porosity  and permeability tensor K (or transmissibility T ) is upscaled to coarse-scale effective porosity  * , and effective permeability K * ( or T * ) respectively.In other cases where the two-phase flow parameters, such as relative permeability or/and capillary pressure, are upscaled, the techniques used are referred to as the two-phase parameter upscaling, or simply two-phase upscaling. 5uring practical reservoir simulations, single-phase static upscaling techniques are essential tools and can successfully be applied to develop reasonably accurate coarse-scale models for multiphase flows. 6The volumetric weighted arithmetic mean is commonly used to obtain the upscaled coarse-scale effective porosity  * .For permeability upscaling, several methods have been introduced in the literature.Some of which are analytical and others numerical. 2,7,8Analytical methods (also called averaging methods) such as the arithmetic, harmonic and geometric methods have the advantage that they are usually faster than numerical algorithms; however, they have the disadvantage that when applied to a sector outside their strict domain of validity, they quickly become inaccurate.On the other hand, numerical methods are based on the numerical solution of the pressure distribution and are more accurate but may become computationally too intensive and slow. 4n this article, an effective and efficient model for permeability upscaling will be proposed and investigated.We explore an adaptive approach where a reduced basis method (i.e., proper orthogonal decomposition (POD) and nonnegative proper orthogonal decomposition (NPOD)) is combined with an interpolation algorithm.Here, the central idea is to find an orthogonal transformation between the original permeability space (high fidelity permeability data) and a low dimensional uncorrelated permeability space (low fidelity data set).Then, the low fidelity data is remapped to the coarse simulation grid using an interpolation algorithm.The interpolation algorithm is embedded within the flow simulator.
The main aim of this paper is to extract the primary information representative of the typical characteristics of the permeability distribution from the high-fidelity data and represent it as a new set of independent variables of the POD/NPOD modes from which we obtain low-fidelity data (LRA).As an advantage, the upscaling method presented is data-driven.Like the analytic methods, our models are not based on the numerical solution for the pressure distribution; hence they are faster.Also, there is minimal redundancy within the low fidelity data used during the dynamic coarse-scale simulation since the method ensures that significant attributes of the high fidelity data are recognized, and the dimension is reduced such that information loss is minimal.The process of LRA generation is essential because it is not very likely that all the row/columns in high-fidelity data K are linearly-independent.Such collinearity may lead to instabilities in the solution space, producing inconsistent results.Also, in zero permeability domains, simple averaging methods are known to produce errors in the upscaled values.The LRA seeks to overcome this problem since the dimension reduction is made in the transform domain, not the spatial domain.Finally, it can be applied to upscaling tensors.
This proposed methodology presents two significant improvements over analytical methods: (i) It is scalable and allows for qualitative and quantitative analysis of the discrepancies between reference and scaled-up data of the geological reservoir even before dynamic simulations are performed.(ii) Dimension reduction using the POD model is in the transform domain rather than the spatial domain, making it possible to overcome the problem of geological discontinuities at the boundaries associated with analytical methods.The transform domain method can overcome these issues by introducing geological spatial correlations into the scale-up process.This paper is organized as follows.In Section 2, we present the governing flow equations.Section 3 introduces mathematical tools that will be used to build the upscaling model.Next, we describe the architecture and experimental setup of the model in Section 4. Results and discussion is presented in Section 5. Finally, in Section 6, we present the conclusion.

DESCRIPTION OF THE FLOW EQUATIONS
Fluid flow in porous media are described by a set of partial differential equations representing conservation of mass and momentum, as a function of pressure, saturation, and velocity.
For an immiscible multiphase flow, the momentum conservation equation is described by the extended Darcy's law in the form of, where  designates a phase.u is the saturated Darcy velocity.K is the absolute permeability tensor of the porous medium.K r , , and p are the relative permeability, viscosity and pressure respectively.Finally  mom is the source/sink term.
Here, fluids are assumed to be incompressible, then from the continuity equation for multiple phases, the saturation equation can be obtained in the form where  is the porosity.S  and S cty, are the saturation and source/sink term of the th phase, respectively.The saturation equation is bounded by the constraint ∑ N  =1 S  = 1, 0 ≤ S  ≤ 1, where N  denotes the number of phases.A control volume finite element method (CVFEM) formulation based on dual consistent CVFEM representation, embedded in the families of triangular and tetrahedral finite element-pairs P n DG − P m and P n DG − P m DG, is used to discretize and solve the set of equations.The numerical formulation is fully described in Gomes et al. 9 When discretizing the equations above, fine-scale mesh resolutions are required to fully capture the reservoir heterogeneities with each sub-grid containing reservoir parameters like permeability and porosity.Running simulations with these fine-scale models is computationally expensive.Therefore, some form of upscaling is required to transfer relevant features from the fine-scale model to a coarse-scale model to reduce the computational burden.
Upscaling volumetric (additive) properties like porosity is relatively straightforward.The effective porosity is the bulk volume weighted arithmetic average of the sub-grids.On the other hand, upscaling permeability is a complicated matter, as it is not an additive variable (i.e., the equivalent permeability in the reservoir scale cannot be represented by arithmetic means).
The permeability K may be defined as the ability of the fluid to flow through the pores connectivities in the porous media and is also the most important parameter when simulating in a reservoir.Therefore reliable permeability upscaling methods are required to accurately capture the heterogeneous structure of the fine-scale reservoir model when represented by a coarse-scale model.
Thus this paper will focus on developing a fast and reliable upscaling technique for the permeability tensor K.When upscaling the permeability of a grid block Ω, our goal will be to obtain an upscaled (effective) permeability K * which gives a very similar total flux in the grid block as the original permeability K, given the same pressure gradient.

REDUCE ORDER MODELING VIA POD
This section briefly introduces the proper orthogonal decomposition method and other essential tools used to buildup our proposed reduce order model.

Proper orthogonal decomposition
Proper orthogonal decomposition (POD) is a powerful and elegant modal decomposition method of data analysis aimed at obtaining low-dimensional approximate descriptions of high-dimensional processes.POD has been developed by several authors and known by different names depending on the field of studies. 10The POD technique provides an algorithm to decompose a given data into a minimal number of basis functions or modes to capture as much energy as possible.
In fluid dynamics, it was first introduced by Lumley 11 as a mathematical technique to extract coherent structures from turbulent flow fields, it has also been used in vibration analysis, 12 data analysis or compression, 13,14 damage detection 15 and as tool in deriving reduced order models for the shallow water equations. 16n Reference 17, an enhanced POD model known as weighted proper orthogonal decomposition was introduced and used to numerically investigate the evolution of the swirling flow exiting the hydraulic turbine runner as the turbine discharge is modified.Reddy et al. 18 presented a constrained POD-based ROM to investigate propagation and compounding of error through the time domain of the first-order wave equation. [21][22]

POD setting
Given a non-negative data matrix X > 0 containing n rows and m columns, that is, The rank of X is defined as the maximum number of linearly independent column/row vectors in the matrix.We let N be the rank of X, such that N ≤ min(n, m).POD seeks to factorize X ∈ R m×n in the form where U = (u 1 , u 2 , … , u d ) ∈ R m×d contains the orthogonal basis vectors called the modes (or factors) and R n×d is the associated coefficient matrix.d is a specified rank for matrices U and V which is normally much lower than the rank of X (i.e., d ≪ N).The matrices U and V are obtained by solving the following minimization problem; Here is the square Frobenius norm and X is the low rank approximate (LRA) of X.The solution to the minimization problem relies on the covariance matrix R ∈ R n×n defined by Note that, R is a symmetric positive semidefinite matrix with real, nonnegative ordered eigenvalues  1 ≥ … ≥  n ≥ 0 and corresponding eigenvectors denoted by u j (j = 1, … , n).Due to the structure of R we can choose the eigenvectors as the orthogonal basis of X, that is, u j (j = 1, … , d) are the modes in Equation ( 3) above and V = U T X.
When the minimum is taken over all U, V of dimension d, 23 then it holds that: min When applying POD to practical problems, the choice of the d most relevant basis vectors is of central importance.A strategy commonly used for choosing d is the energy criterion (E C ). 24 Here we assess the variance, which is the ratio between the modeled variance and the total variance contained in An adhoc rule frequently applied, the closer E C is to 1, the better.Here, the realization of POD is identical to PCA (principal component analysis), and its application in this work is to reduce the dimension of the static model.We used POD in this work to reduce the dimension of the geophysical property.

A constrained proper orthogonal decomposition model
Although the POD model is optimal for low rank approximation, a major drawback is that the POD mode do not necessarily conform to specific constraints, such as non-negativity.So when non-negativity of the data is desirable because the underlying factor has a physical interpretation as in the case of upscaled permeability, an additional constraint will have to be imposed on Equation ( 6) to ensure X ≥ 0.
Adding the non-negativity constraint (U, V ≥ 0 ) to Equation ( 6) we arrive at the following optimization problem min A direct byproduct of the combination of U T U = I and U ≥ 0 entails that U is disjoint, meaning that each column of U contains at most one non-zero element.While having disjoint modes may be considered as a kind of sparseness, it is too restrictive for most problems (i.e., as in the case of permeability upscaling a sparse solution is not desirable).We therefore wish to allow some overlap among the modes (i.e., column of U) to relax the disjoint property.The degree of coordinate overlap can be represented by an orthogonality distance measure ||U T U − I|| 2 which vanishes only when U is orthogonal. 25he relaxed version of Equation ( 10) becomes: min where  > 0 is a balancing parameter between reconstruction and orthogonality.Now, the factor matrix U is a near-orthogonal nonnegative matrix, that is, U T U ∼ I. 26 Equation ( 11) describes the constrained POD method, which we refer to henceforth as the non-negative POD (NPOD) model.The NPOD algorithm straightforwardly yields an algorithm identical to the known orthogonal non-negative matrix factorization (ONMF). 27To solve Equation (11), an iterative two-step strategy was adopted by alternatively optimizing U and V.At each iteration, one of the matrices is optimized while the other is fixed.The resulting multiplicative update rules are obtained as where ⊙ is the Khatri-rao product.The correctness of these solutions for U and V is assured by the fact that at convergence they will satisfy the Karush-Kuhn-Tucker (KKT) conditions for Equation (11).

Karush-Kuhn-Tucker (KKT) conditions
The Karush-Kuhn-Tucker (KKT) conditions are used to obtain the solution for optimization problems constrained to one or more equalities.Consider an optimization problem: where f is the objective function, h i (x) is the equality constraint and l j (x) is the inequality constraint.Assuming the equation is differentiable, there exist constants u i (i = 1, … m) and v j (j = 1, … , l) called KKT multipliers that satisfy the Equations ( 15)-( 18) are the stationary, complementary slackness, primal feasibility and dual feasibility conditions, known as the KKT conditions. 28In the absence of the inequality constraints in Equation( 14), the KKT conditions decompose to the Lagrange conditions, and the KKT multipliers are called the Lagrange multipliers.

Derivation of multiplicative updates
The core idea for deriving the multiplicative updates in Equation ( 12) and ( 13) is based on the KKT conditions described above.We start by defining the objective function 29 Then where tr(A) represents the trace of A. Let Ω ij and Λ jl be the Lagrange multiplier for the constraints U ij ≥ 0 and V jl ≥ 0 respectively.The Lagrangian function  for the two multiplier metrics Λ = Λ ij and Ω = Ω ij is expressed as For optimality, the KKT conditions requires Also from the KKT complementarity slackness condition we have Incorporating Equation ( 22) and ( 23) into ( 24) and (25), respectively leads to Finally separating negative and positive parts of the gradient leads to the multiplicative update rules in Equation ( 12) and ( 13).

Re-mapping via interpolation
During upscaling, the LRA X is remapped from the fine grid to the coarse grid using nearest neighbor interpolation (NNI).Any interpolation methods aim to use a source dataset as the reference to construct a new scaled/interpolated dataset.The new or constructed dataset depends on the interpolation ratio selected.Nearest neighbor is one of the simplest interpolation techniques and is very useful when computational speed is an area of concern.This method sets the value of an interpolated point to the value of the nearest existing data point.Therefore, for two successive data points (x k , x k+1 ) ∈ X, this method first find the mid value using these data points. 30Then values of x, smaller than the mid value, corresponds to the value of y k and values greater than the mind value corresponds to the value of y k+1 .The kernel function f (x) which define the nearest neighbor interpolation technique, is a zeroth-degree polynomial and is given by,

SIMULATION MODEL SETUP
This work introduces a novel ROM for permeability upscaling.Here, we explored an adaptive approach where the NPOD model combined with the nearest neighbor interpolation results in a coupled model (i.e., NPOD-NNI reduce order model).
The model involves two stages.First, the permeability tensor K, a diagonal tensor {K x (i, j), K y (i, j)} with i = 1, 60 and j = 1, 220, is obtained from fine-scale reservoir model (SPE10 model2 31 ).The NPOD model is applied to K to generate low rank approximation K.The objective of the NPOD model is to extract the primary information from a large amount of data and represent it as a new set of independent variables (modes).Thus using the NPOD model, the key attributes of the permeability data are recognized, and the dimension is reduced, ensuring minimal information loss.This process is essential because it is not very likely that all the row/columns in K are linearly-independent.Such collinearity may lead to instabilities in the solution space, producing inconsistent results.
In the second stage, the LRA K is re-mapped from the fine-grid to the coarse-grid using the NNI algorithm.The NNI algorithm is embedded in the simulation codes within the ICFERST simulator.When running the simulation, high-order numerical methods are used to solve the flow equations described in Section 2. Figure 1 illustrates the process of the experimental model setup.

SPE 10th comparative solution project: Model 2
The SPE10 reservoir model 31 was proposed as a benchmark for upscaling techniques and is, therefore, a good test case for our methodology.Although the model is structurally simple, it is highly heterogeneous (see Figure 2).The geological model consists of part of a Brent sequence which includes two units.The top part of the model (made up of the first 35 layers) belongs to the Tarbert formation.Its lower part (made up of the remaining 50 layers) is an Upper Ness sequence (fluvial).Both formations are characterized by large permeability variations, but are qualitatively different as shown in Figure 3.
The Tarbert formation has relatively smooth permeability variations and has good communication in the vertical and horizontal directions.However, the Upper Ness formation is fluvial and contains channels making the formation particularly hard to upscale.From each formation, a layer will be extracted and studied.

Goodness of fit
"Goodness of fit" (gof) statistics are presented to validate the model's performance.For the static model, qualitative analysis on the LRA K reconstructions are performed.In the dynamic model, we analyse the responses (such as flow behaviors) of the fine-scale and corresponding coarse-scale model.The relative root-mean-square error is used to measure the quality of LRA K reconstructions and evaluate the performance of the different models.
where X and X are the observed and expected data, respectively, and N t is the total number of elements in X.In order to measure the amount of information loss during the LRA reconstructions, we used the Kullback Leibler (KL) Divergence, defined as The significance of the D KL (P||Q) can be tested through its relationship to the minimum discrimination information (MDI) statistics.
where N = ∑ x.The MDI is Asymptotically chi-square distributed. 32inally, to measure the performance of the model as the number of the modes increases, we define the explained variance score by where Var is the variance.To select the number of retained modes to use in the dynamical simulation, we compare the cumulative variance ratio, and the relative root mean square error of LRA reconstructions.We determine that the minimum number of the POD modes acceptable is such that the conditions EVS ≥ 97%, and RRMSE ≤ 10E − 1 (34) hold.

SIMULATION RESULT
To validate the NPOD-NNI reduce order model for upscaling, two 2-dimensional heterogeneous permeability data from the top and bottom layers of the SPE 10th comparative solution project 31 are considered.For both test cases, the computational domain entails 6 m × 11 m with 13,258 triangular elements at fine-scale.Two wells (injector and producer) are located at opposite ends (i.e., left and right end) of the reservoir.In all test cases at initial conditions, the computational domain contains two fluids: fluid 1 with a saturation of 20% and fluid 2 with a saturation of 80%.Fluid 1 is injected with a constant flow velocity (0.3 ms −1 ) through the left boundary displacing fluid 2 towards the right boundary.Pressure is kept constant at the production end.The relative permeability is evaluated using the modified Brooks-Corey model. 33Results from low-rank approximation using the POD-based (i.e., POD and NPOD) models are first analyzed.Next, numerical results from the POD-based upscaling models and an analytic upscaling method will be compared to those obtained from fine-scale simulations.

Test case1: SPE10-Layer20
In this test, we used permeability data obtained from the 20th layer of the SPE10 model.This layer, as shown in Figure 4, is an extract from the Tarbert formation having good communication in the vertical and horizontal directions.The data described by a 60 × 220 matrix is decomposed using the POD and the proposed NPOD algorithm.The modes are extracted and used in reconstructing the upscaled permeability.Analysis of the quality of reconstructions using different numbers of modes is observed.In Figure 5A, as the number of modes used in reconstruction increases, the error asymptotically tends to zero.Theoretically, all information and characteristics of the original data would be recovered when all modes are used in reconstruction using the POD algorithm; however, this is not the case for the NPOD model.This is because the optimization problem in the NPOD model described in Section 3.2 has no exact solution, and only approximate solutions are obtained using the KKT conditions.Hence, the error from the NPOD model will only approach zero but is never equal to zero even if all modes are used in the reconstruction.Figure 5B shows the information loss during reconstructions.As the number of modes increases, the information loss is minimized, tending asymptotically to zero.
Although from statistics of the error analysis for the statics model, the POD algorithm might seem more preferred for LRA, as stated before, when non-negativity of the data set is a priority, the NPOD should be more appropriate, as will be shown in Figure 6.The NPOD model can use a minimal number of modes to reconstruct the permeability image while retaining the geological details from the original data as opposed to the POD model.Reconstructed images using different modes are generated using the POD-NNI (coupled POD with nearest neighbor interpolation) and NPOD-NNI models, respectively.The LRA's obtained from the POD algorithm are refined by taking the absolute value of the LRA data to ensure non-negativity.The permeability distribution in Figure 4 has 13,258 triangular elements, which is the reference image.Figure 6 shows images of reconstruction with 8200 triangular elements obtained from the POD and NPOD models using five, twenty-five and forty-two modes.The reference image has a rank of 60.In low-rank approximation generation, the number of modes used in reconstruction is equivalent to their ranks.Thus the LRA reconstructed with 5, 25, and 42 modes correspond to having ranks 5, 25, and 42, respectively.
Clearly, by increasing the number of modes (rank), the image gradually changes as they become sharper, smooth out at the contours, and tones are marked, highlighting regions of high and low permeability more precisely.
The choice of the optimal rank for the LRA reconstruction is essential since if we go below a certain threshold, the features of the original map will not be recovered.In this work, the optimal rank is selected by performing POD-based models for a range of ranks and then choosing a rank based on the metric of goodness described in Section 4.2.
From Figure 7, the explained variance score is at EVS = 99% when the mode number reaches 42, and from Figure 5A, the RRMSE for the LRA reconstruction using 42 modes is 4.627E − 2. Therefore, from Equation (33), it can be confirmed that a high reconstruction accuracy is achieved if the top 42 modes of layer 20 are adopted in computations for the upscaling.
Table 1 shows the basic statistics of the LRAs reconstructed with 42 modes.As statistics of the permeability dataset produced by the POD-based models indicate similarity with the original data set (static model), the next step is to conduct numerical simulations to investigate the flow dynamics using the full (i.e., original ) and reduced(upscaled) data set.
For running simulation, the ICFERST 1 flow simulator is used.From Layer 20 data set, LRA K using 42 modes are obtained and mapped to the coarse-grid using the nearest neighbor algorithm.We simulate the fine-scale test case and compare this to the coarse-scale test case.Figure 8 shows the permeability distribution and saturation flow profile of the reference(fine-scale) test case and coarse test case using the NPOD-NNI model.
Here, we scaled up the reference test case from 13,258 elements to 1004 elements.Clearly, at the coarse level, the structure of the fine-scale permeability distribution is still preserved.The high and low permeable zones remain detectable.Also, we note a similar flow behavior between the responses at the coarse-scale and the corresponding fine-scale test case.
Figure 9 shows the comparison plot of the cumulative production rate of fluids 1 and 2 using the POD based models and the arithmetic mean (AM) upscaling model.Noticeably, from the curves, we observed a good match for the results from the forward simulations between the proposed upscaling models and the reference test case.Statistically comparing the performance of reduced-models in reproducing the original response of the references case, the RRMSE in table 2 shows that the NPOD model results in considerably less error than the POD model in the upscaling process, albeit both POD and NPOD models performed better than the classic arithmetic mean upscaling model.Also, as expected, the computational time is significantly reduced when using the upscaled models compared to the reference case.
Comparing the production rate curves in Figure 9 and the results in Table 2 clearly illustrate the advantage of using POD-based models over the classic arithmetic mean upscaling model as the performance of the POD-base model is more accurate with a better curve fit.Even more remarkably, Figures 8 and 9 indicate clearly that the NPOD-NNI algorithm could be used to represent the fine-scale model as fluid flow behavior are very similar.

Test case2: SPE10 layer 81
In the second test case, permeability data is obtained from the 81st layer of the SPE10 model.The layer, as shown in Figure 10, is an extract from the Upper Ness formation.This formation is hard to upscale because it contains channels resulting from rivers or running water deposits in a delta-plain environment.Like in the first test case, similar analyses are carried out to investigate the performance of the upscaling model on complex reservoirs.The permeability, described by a 60 × 220 matrix, will be decomposed using the POD and proposed NPOD algorithm.The modes are extracted and used in reconstructing the upscaled permeability.
Figure 11 shows the associated error and information loss during reconstruction with different numbers of modes.As the number of modes used in the reconstruction increases, both the error and information lost asymptotically tend to zero.The data can be said to have some sparseness since it is extracted from a channelized region with permeabilities having significant variations in the orders of magnitudes.Thus, the correlation between the row/column of the data matrix is minimal, with almost all modes required in the LRA generation to minimize the error and associated information loss.Using the POD-NNI and NPOD-NNI, the reconstruction images from the different modes are generated.Figure 10 shows the permeability distribution of SPE10 layer81 on a 13,258 triangular element, which describes the reference data set used in this test.Figure 12 shows the images of reconstruction with 8200 triangular elements, obtained from the POD and NPOD models using five, twenty-five, and fifty-five modes.
The reference image described in Figure 10 has a rank of 60.The LRA reconstructed with 10, 25, and 42 modes correspond to ranks 10, 25, and 42, respectively.By increasing the number of modes (rank), the image gradually changes, highlighting regions of high and low permeability more precisely.Estimating the number of modes to use LRA reconstruction in the dynamic simulation will be determined based on Equation (33).
From Figure 13, the explained variance score is at EVS = 97% when the mode number reaches 51, and from Figure 11A, the RRMSE for the LRA reconstruction using 51 modes is 1.07E − 1.Therefore, from Equation ( 33), it can be confirmed that a high reconstruction accuracy is achieved if the first 51 modes of layer 81 are adopted in computations for upscaling.We also note that when using the POD model, the cumulative energy contribution is at E c = 97% when the mode number reaches 44.
Table 3 shows the basic statistics of the LRAs reconstructed with 51 and 44 modes for the NPOD and POD models, respectively.As statistics of the permeability dataset produced by the upscaling models indicate similarity with the original data set (static model), the next step is to conduct numerical simulations to investigate the flow dynamics using the original (i.e., high-fidelity data) and reduced (low-fidelity data) data set.
From Layer 81 data set, LRA K using 51 modes are obtained and mapped to the coarse-grid using the nearest neighbor algorithm.We simulate the fine-scale test case and compare it with the coarse-scale test case results.The permeability distribution and saturation flow profile of the fine-scale and coarse-scale test cases using the NPOD-NNI are shown in Figure 14 10970363    Here, we coarsened the original mesh from 13,258 elements to 1004 elements.The original fine-scale permeability distribution structure is slightly preserved at the coarse scale.Also, there is a similar flow behavior between the saturation profile of the coarse scale and the corresponding fine-scale model.
Figure 15 shows the comparison plot of the cumulative production rate of fluids 1 and 2 using the different upscaling models.Noticeably, from the curves, we observed a good match for the results of dynamic simulations between the proposed models and the reference test case compared to the arithmetic mean upscaling model results.Statistically comparing performances of reduced-models in reproducing the original response of the references case, the RRMSE in Table 4 clearly shows that the POD-based models result in considerably less error than the arithmetic mean model in the upscaling process.Furthermore, the CPU times for running the dynamic simulations using the upscaling models are similar and much less than the reference model.
The cumulative production curves in Figure 15 and the results in Table 4 clearly illustrate the advantage of using POD-based models compared to the existing analytic model for upscaling.More specificly, the NPOD

CONCLUSION
Proper orthogonal decomposition based reduce order models have evolved during the past decades as a promising model order reduction technique for multiphase flow in large-scale heterogeneous systems.This work presents a novel reduce-order modeling technique for upscaling heterogeneous domains.The main objective was to extract the primary information representative of the typical characteristics of the permeability distribution from the high-fidelity data and represent it as a new set of independent variables of the POD/NPOD modes from which low-fidelity data (LRA) are determined.Next, we conduct simulations involving high-order accurate numerical methods using this low-fidelity data on coarse grids to investigate the multifluid flow behaviors.The new ROM technique was evaluated using extracts from the SPE10 benchmark, and the results were compared with an existing analytical upscaling method.As the numerical tests show, the upscaling model preserves the geological structures and flow behavior from the fine-scale model.Moreover, the technique introduced significantly reduces the redundancy of the high-fidelity data and provides qualitative and quantitative analysis of the static data.In future work, we will extend the POD/NPOD to multilinear algebras to investigate high-order POD/NPOD reduce order models for upscaling 3D domains.

F I G U R E 1
Pseudocode describing the model setup F I G U R E 2 SPE 10 model showing the spatial distribution of permeability [Colour figure can be viewed at wileyonlinelibrary.com]

F I G U R E 3
Histogram of rock permeability (log K x ) for the SPE 10 model.The Tarbert formation is shown in blue and the Ness formation in red [Colour figure can be viewed at wileyonlinelibrary.com]

11 F I G U R E 4
Permeability distribution ( in m 2 ) of layer 20 (6×11 m) of the SPE10 Benchmark [Colour figure can be viewed at wileyonlinelibrary.com]F I G U R E 5 The RRMSE and Information loss using different modes for LRA reconstruction of the permeability distribution of layer20 [Colour figure can be viewed at wileyonlinelibrary.com]Porosity and viscosity are set at  = 0.2 and  = 1cP respectively.Boundary conditions of no-flux across upper and lower borders were set.

6
Images of reconstruction with 8200 triangular elements obtained from the POD and NPOD models using some selected number of modes for layer20 [Colour figure can be viewed at wileyonlinelibrary.com]

F I G U R E 7 F
Variiance explained by each of the modes of K x [Colour figure can be viewed at wileyonlinelibrary.com]TA B L E 1 The mean, variance and standard deviation of the different datasets to be used in the dynamic simulation of layer20 I G U R E 8 Top: fine-scale permeability distribution on a 13258 mesh grid and corresponding flow profile.Bottom: coarse-scale permeability distribution on a 1004 mesh grid calculated by the NPOD model using 42 modes for LRA reconstruction and corresponding flow profile [Colour figure can be viewed at wileyonlinelibrary.com]

F I G U R E 9
Cumulative production curves for fluid 1 and 2 for the fine-scale reference test-case and coarse-scale test-cases using NPOD-NNI and POD-NNI [Colour figure can be viewed at wileyonlinelibrary.com]TA B L E 2 The CPU times for the dynamical simulation carried out on an x86 64 quad-core machine running Linux and RRMSE of the cumulative production rate of fluid 1 and fluid 2 from performances of reduced models in reproducing the original response of the references case for layer 20

F I G U R E 10
Permeability distribution ( in m 2 ) of layer 81 (6×11 m) of the SPE10 Benchmark [Colour figure can be viewed at wileyonlinelibrary.com]F I G U R E 11 The RRMSE and Information loss using different modes for LRA reconstruction [Colour figure can be viewed at wileyonlinelibrary.com] U R E 12 Images of reconstruction with 8200 triangular elements obtained from the POD and NPODbased models using some selected number of modes for layer81 [Colour figure can be viewed at wileyonlinelibrary.com]F I G U R E 13 The variance explained by varying numbers of modes for K x in layer81 [Colour figure can be viewed at wileyonlinelibrary.com]TA B L E 3 The mean, variance and standard deviation of the different datasets to be used in the dynamic simulation of layer81 U R E 14 Top: fine-scale permeability distribution on a 13,258 mesh grid and corresponding flow profile.Bottom: coarse-scale permeability distribution on a 1004 mesh grid calculated by the NPOD model and corresponding flow profile [Colour figure can be viewed at wileyonlinelibrary.com]F I G U R E 15 Cumulative production curves for fluid 1 and 2 of the fine-scale reference test-case and coarse-scale test-cases using upscaling models [Colour figure can be viewed at wileyonlinelibrary.com] 10970363, 2023, 6, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/fld.5171by University Of Aberdeen, Wiley Online Library on [15/05/2023].See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions)on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 10970363, 2023, 6, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/fld.5171by University Of Aberdeen, Wiley Online Library on [15/05/2023].See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions)on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 10970363, 2023, 6, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/fld.5171by University Of Aberdeen, Wiley Online Library on [15/05/2023].See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions)on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License , 2023, 6, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/fld.5171by University Of Aberdeen, Wiley Online Library on [15/05/2023].See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions)on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 10970363, 2023, 6, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/fld.5171byUniversityOf Aberdeen, Wiley Online Library on [15/05/2023].See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions)onWileyOnline Library for rules of use; OA articles are governed by the applicable Creative Commons LicenseTA B L E 4The CPU times for the dynamical simulation carried out on an x86 64 quad-core machine running Linux and RRMSE of the cumulative production rate of fluid 1 and fluid 2 from performances of reduced-models in reproducing the original response of the references case for layer 81 better choice for upscaling as the performance of the NPOD model is more accurate with a better curve fit.