Regular Article
Hierarchical simultaneous upscaling of porosity and permeability features using the bandwidth of kernel function and wavelet transformation in two dimensions: Application to the SPE10 model
^{1}
Mining Exploration Engineering, Shahrood University of Technology, Shahrood 9761866466, Iran
^{2}
School of Mining, Petroleum and Geophysics, Shahrood University of Technology, Shahrood 3619995161, Iran
^{3}
The Center of Excellence in Mining Engineering, Shahrood University of Technology, Shahrood 3619995161, Iran
^{4}
Department of Statistics, Faculty of Mathematical Sciences, Ferdowsi University of Mashhad, 9177948974, Iran
^{*} Corresponding author: azad66.mohammad@gmail.com
Received:
7
November
2020
Accepted:
17
February
2021
In this paper, two methods of kernel bandwidth and wavelet transform are used for simultaneous upscaling of two features of hydrocarbon reservoir. In the bandwidth method, the criterion for upscaling is the cell variability, and by calculating the optimal bandwidth and determining the distance matrix, the upscaling process is performed in a completely nonuniform and unregularly manner. In areas with extreme variability, the bandwidth is considered small enough to maintain the fine scale characteristics of model. Conversely in homogenous areas, with the choice of large bandwidth, the maximum rate of upscaling will occur. The bandwidth upscaling algorithm is an iterative and hierarchical algorithm. The bandwidth method, unlike conventional scaleup methods, focuses on how to upgrid cells and, by determining the optimal averaging window, we will have the least loss information for the fine scale model. Upscaling is a preprocessing to building a simulator model with lower cell number, and thus, reducing volume and computational cost, while maintaining and retaining the basic information of the fine model. Due to the various variability of the reservoir features, the attribute upscaling pattern differs, and in order to show the variability of two features in the upscaling model simultaneously, it is suggested in this paper to upscale two features simultaneously. For simultaneous upscaling, we applied two different approaches; minimum and maximum bandwidth. Moreover, wavelet transformation is applied to upscaling the model. Then, as a result, the variance of the scaleup models based on wavelet is about onethird of the variance of the bandwidth method. Simulation results show that the bandwidth method is a good approach for upscaling the heterogeneous reservoirs.
© M.R. Azad et al., published by IFP Energies nouvelles, 2021
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 Introduction
Reservoir simulation plays an important role, and is used to predict and analyze fluid for decisions and management of oil production. Models for geological characterizations typically contain many cells. The number of cells is typically an order of 10^{7}−10^{10} cells, which are referred to as fine scale models. These models represent geological variation on fine scales. Upscaling of the geological fine model is required because computational costs can make it unpractical to perform compositional simulations in finescale models. Therefore, upscaling of finescale geological models is important and the use of coarsescale models in reservoir simulation is necessary to reduce computational time [1]. To develop a reducedorder model from a fine scale model, first, we must create a grid model with reduced spatial resolution, and then, capture reservoir features to a new coarse scale model using a suitable upscaling method. Finally, we have to allocate to each coarse cell an equivalent value [2–7].
Based on different factors, there are many upscaling methods in the literature. For example, according to the parameters that should be upscaled, upscaling techniques are classified in singlephase and twophase flow parameter upscaling. Singlephase methods, which are used for upscaling of static parameters such as porosity and absolute permeability, consist of analytical methods, flow based methods, wavelet, renormalization, etc. For example, Renormalization is a suitable upscaling technique when connectivity is considered as a feature [8]. However, in twophase flow case, relative permeability can be upscaled as a property. These methods include pseudo methods, dynamic pseudo, and verticalequilibrium [2, 4, 5]. Preux et al. considered some upscaling techniques and by comparing their results chose the appropriate upscaling method [8].
Geological models have some complexities in order to fluid flow simulation. In this situation, it is necessary to look at the issue from a multiscale perspective. In this way, many goals will be pursued: data integration, reduce uncertainty as well as reduce computational costs. Noetinger and Zargar (2004) have studied the importance of having a multiscale qualification of fluid flow simulation from the perspective of a reservoir engineer in a comprehensive review. They have considered all relevant problems from the geological model to the simulation model [9]. Considering the importance of identifying discontinuities and heterogeneity in reservoirs, traditional upscaling methods do not have the required efficiency in complex and heterogeneous reservoirs, and it is necessary for upscaling based on the variability. Among upscaling methods that work based on variability, we can refer to the wavelet transform and the bandwidth of the kernel function. Wavelet transformation works in a nonuniform and intelligent way with a multiresolution upscaling approach, in which areas with severe variability will remain fine and areas with moderate smoothness will be coarse [10–12]. The application of wavelet transforms to the coarse graining of flow in porous media is not a new concept and in these studies; a wavelet transform is applied to the reservoirs to map the important information at diverse scales with a set of approximations and detail coefficients [13–18]. Bouard et al. used with the separable 3D wavelet transform to decomposition properties such as porosity and they evaluated the simulation results in different decomposition levels [19]. The wavelet drawback is that the upscaling process is done in binary way. For example, in two dimensions, four cells are either merged or not, and there is no state between them. This problem can be solved by the bandwidth method.
Moreover, the main challenge in preparing upscale models is to optimize the upscaling method. The global upscaling method is composed of upgridding and upscaling. A procedure of selecting the appropriate option in merging the cells is called upgridding. Once the option is chosen, the standard procedure of upscaling is performed. Therefore, the main challenge in the upscaling process is how to upgrid the cells, otherwise the essence of any upscaling method is averaging [20]. If the upgridding done is based on cell variability, the coarse scale model will have the least difference with the fine scale model. In this paper, we propose a new method based on bandwidth of kernel function in which upgridding will be a function of cell variability. Bandwidth will determine the development of scaleup cells in a coarse scale grid.
In order to determine the number of cells in the coarse scale model, a proper division based on the cell variability must be made. This process can be likened to the determination of bandwidth in the kernel function method. In the proposed method, upscaling is a function of bandwidth, and the bandwidth is a function of cell variability [21]. The basis for upscaling will be the variability. The method is based on the distance of adjacent cells. This method searches for cells with the most similarity, and then, based on the defined bandwidth; the upscaling process will be performed. The main difference between this method and other upscaling methods is to scale up in the manner of the upgridding. When upgriding is right, scaling up will have the best result. Moreover, in this method, the scale does not follow a particular trend, and nonuniform coarsening is applied in every porous medium.
In other hand, in fluid flow simulation, the upscaled model for all reservoir properties, including permeability, porosity, saturation, type of rock, etc., is considered the same, regardless of the spatial distribution of different properties. The behavior and distribution of the variability of two distinct properties, such as permeability and porosity, are different in the reservoir, thus, their upscaling pattern must also be different. Since all parameters are used in fluid flow simulation, it is therefore necessary to calculate the upscaling pattern based on the variability of each property. In this paper, simultaneous upscaling of the two features or properties is used to show the effect of the variability of the two features of the reservoir in the upscaled model. Hierarchical multiresolution simultaneous upscaling of two features based on the kernel function bandwidth allows the ability to achieve the same pattern for both features, so that the number of cells and the placement of fine and coarse cells is the same for them, and only the difference between the models is the quantitative amounts of the features. In this regard, Azad et al. [22] first used the bandwidth method for simultaneous upscaling of two petrophysical properties of the reservoir in one dimension. Using different approaches, they calculated simultaneous upscaled model for two properties [22]. Moreover, Azad et al. [23] compared two upscaling methods: the bandwidth and wavelet transformation for simultaneous upscaling of two petrophysical properties of the reservoir in one dimension. They showed that the bandwidth method was more accurate than the wavelet method in order to maintain the basic information of the reservoir [23].
The main purpose of this paper is simultaneous upscaling of two features in reservoirs, particularly permeability and porosity, using the kernel bandwidth method. The proposed upscaling method in this paper is applied for SPE10 model [2]. The model has a simple geometry to provide maximum flexibility in comparing upscaled grids. At the fine scale model, it is described on a regular Cartesian grid. The model dimensions are 1200 × 2200 × 170 ft, which constitutes two different formations: Tarbert Formation, and Upper Ness Formation. The fine scale cell size is 20 × 10 × 2 ft, and the model consists of 60 × 220 × 85 cells (1.122 × 10^{6} cells). Figure 1 shows the porosity for the fine model, which varies in range from 0% to 50%. The reservoir contains channeling paths, with perpendicular permeability to the horizontal one equals 0.3 and in other areas; this value is up to 0.0001 declines [24, 25].
In summary, we have followed two main new strategies in this research: (1) we applied a new method in upscaling that is based precisely on the cell variability (Bandwidth) and then we have compared these results with a wellknown method in upscaling that is referred to wavelet transformation. Bandwidth solves all the disadvantages of wavelet as a common method of upscaling. This method allows the reservoir engineer to execute the scale model and then the simulator model, depending on the number of upscaling cells. As will mention in the following text, by selecting the appropriate bandwidth, the number of enlarged final cells can be controlled. Moreover, the upscaling error is always visible and based on the amount of error, and then we can calculate the upscaled model. (2) Another innovation of this paper, which is more important than the first, is the simultaneous upscaling of two reservoir features (porosity and permeability). Basically, in simulation topics, upscaling is not done simultaneously. Typically, for all attributes with varying variables, a unique upscaling pattern is considered, which in turn causes errors and information loss. In this paper, an intelligent algorithm is used to determine the model of Simultaneous Upscaling of Two Features (SUTF). This means that the local variability of both features is considered and portrayed in this final model. Our focus in this article has been on upgridding, otherwise upscaling is nothing more than averaging.
2 Research methodology
2.1 Bandwidth method
Kernel density estimation is an important smoothing technique. It has been applied to nonparametric problems, e.g. discriminant analysis, goodnessoffit testing, hazard rate estimation, intensity function estimation and regression, in which the main idea is to estimate the density function at a given point using neighboring observations [26–28]. A typical form for the multivariable kernel density function is shown in equation (1).
The data matrix X contains the number of n independent random ddimensional vectors drawn from the density f(x), that must be estimated. Here, the function K is called the kernel and h is the variable bandwidth. Determining bandwidth is a crucial step in estimation [29]. In estimation, in areas with extreme variability, the best estimate can be made by choosing a small bandwidth, and vice versa.
In this paper, we employ the bandwidth for upscaling. The upscaling in this method will be a function of variability. Similar to what is defined in the kernel function estimate for bandwidth, it can also be used to scale up. To determine the number of cells in the upscaled model, a proper distribution based on cell variability should be made. This approach can be used to determine the bandwidth of the kernel function. In areas of the reservoir where the variability is intense, considering the small bandwidth, minimum upscaling will occur, and these areas will remain fine scale. Conversely, in areas with low and smooth variability, considering the choice of large bandwidth, most cells will be merged together.
The main challenge in bandwidth approach is to determine optimal bandwidth. The domain bandwidth changes are determined by the feature changes. For each bandwidth value, the upscaled model will be obtained exclusively. Because the basis of the upscaling is bandwidth, and the cells are merged together when the bandwidth allows, therefore, determining optimal bandwidth is crucial. In order to determine the optimal bandwidth, we have introduced various methods in previous studies [18–20]. One of these methods uses a defined number of cells for the final upscaled model before starting the upscaling process; this is an important distinguishing feature of this method compared to conventional upscaling methods in which the number of cells in the final upscaled model before starting the upscaling process is unknown. Suppose the fine model has 1 million cells. If the goal of upscaling is to provide a simulator model with 100 000 cells, then, this can be done by determining the bandwidth appropriate to this number of cells. In fact, using this method, we can achieve 100 000 cells that have the most variability of the fine scale model.
In twodimensional (2D) case, the bandwidth is defined as a vector (λ_{x}, λ_{y}). For example, two cells in x direction will be upscaled if their difference is less than λ_{x}. After determining the bandwidth, the distance matrix must be calculated. The distance matrix represents the similarity of cells. Given that the purpose of upscaling is to identify areas with the least feature variations, as well as to identify the cells with the most similarity, the distance matrix will be a similarity pattern. The starting point of the upscaling process is the minimum point in this matrix. Each entry of the matrix represents the position of two neighboring cells. If this amount is less than the optimal bandwidth in the specified direction, then the two cells will merge together and a coarse cell will be built. After the first upscaling step, the distance matrix should be updated because the distances will change. The process is hierarchical and repeatable, and when the difference of two cells is more than the bandwidth, the operation will stop. Equation (2), given in the following, is used to calculate the distance matrix:
(2)here D_{x} and D_{y} are the distance matrices in x and y directions and P is desired feature quantity. The smallest value in the matrices is known as the starting point for the upscaling process. Each distance represents the position of two neighboring cells.
Obviously the scaleup process based on the distance matrix and bandwidth will be done in such a way that the cells with the highest degree of similarity will be merged. If the algorithm is stopped, this ensures that the cells are merged together when the difference of the cells is less than the bandwidth. The cells that remain in the fine scale are those that have been highly variable and their distance is far away in the matrix. The twodimensional upscaling algorithm is presented as a schematic flow diagram in Figure 2.
Fig. 2 Twodimensional upscaling algorithm based on the kernel bandwidth. 
2.2 Wavelet transformation
Wavelet Transform (WT) is a method for upscaling that have been used to upscale heterogeneous porous media by coarsening the property field such as permeability and porosity. The method upscales the highresolution geological model nonuniformly so that we can preserve the important information on the spatial distribution of the feature field. Thus, the number of grid blocks in the computational grid and, hence, the number of flow and transport equations to be solved are reduced drastically without loss information about the feature field [11–14, 17, 30]. In the wavelet transform, each signal is divided into two series of approximation and detail coefficients. These coefficients are determined based on scaling and wavelet functions.
Suppose we have a 1D signal with a resolution of 8 data as shown in Figure 3. The Haar wavelet transform is shown in Figure 3. First row is the signal. After decomposition, signal separates in two series coefficients: four approximations and four detail coefficients. Approximation coefficients are calculated based on the mean of samples pairwise, and then the difference between the pair is referred to as detail coefficients.
Fig. 3 (a) Wavelet decomposition and (b) 1D decomposition. 
The 2D Haar wavelet decomposition can be computed using 1D Haar wavelet decompositions because 2D Haar wavelet basis is separable. In two dimensions, we have three wavelet functions and only one scaling function. Wavelet transform attempts to upscale the feature by a factor 2. In 2D square grid, (i_{1}, i_{2}) is the center of block, which has four wavelet coefficients, and consists of three wavelet detail coefficients and one approximation coefficient. The detail and approximation coefficients are calculated using equation (3):
(3)where φ is called the wavelet scaling function and is orthogonal to ψ, which is called wavelet detail function or mother wavelet. In this article, we have used the Haar wavelet as the mother wavelet. In fact, the wavelet analysis is a measure of similarity between the basic functions (wavelets) and the signal itself. In this equation, measures the difference between f(x, y) in the coarser scale and their neighbors in the previous level, with d = 1, 2 and 3 in the x, y and diagonal directions. Moreover, S_{j}(i_{1}, i_{2}) measures the basic information about f(x, y) in the coarser level. The subscript j is the upscaling level.
The main idea of the application of wavelet transform is to minimize the averaging in upscaling. In the waveletbased upscaling, after analyzing the signal to the desired level, the last approximation coefficients and the remaining detail coefficients are introduced as the upscaled signal. For upscaling, we must specify two thresholds ε_{s}, ε_{d} for the approximation and detail coefficients respectively. These thresholds are the key to upscaling wavelet transform. These values determine whether scale coefficients or details are significant. ε_{s} is calculated based on the maximum approximation coefficients that is a measure of the cells’ feature (such as porosity), with which is related a corresponding wavelet scale coefficient, and also ε_{d} is the fraction of the largest detail coefficients which measures the contrast between the porosities or permeabilities of the neighboring cells. To implement the upscaling process, we compare the approximation and detail coefficients with the thresholds. If S_{j} (i_{1}, i_{2}) < ε_{s}, and also , then each of the four cells with the center (i_{1}, i_{2}) will be merged together. If S_{j} (i_{1}, i_{2}) < ε_{s}, there is no need to check the detail coefficients and the cells remain fine [11]. This confirms that the cell’s feature is large enough and we must move to the next cell. Therefore, in this method, only cells whose detail and approximation coefficients are less than the thresholds will be merged. In the later levels of upscaling, the cells that remained fine in the previous level will not enlarge in later levels.
The main challenge in waveletbased upscaling is to determine thresholds. It is crucial to recognize the importance of coefficients. If the threshold is set to zero, no upscaling will be made and conversely, if the thresholds are set to the maximum, the whole grid will be merged into a cell. Many methods have been developed in the threshold setting, one of which is the limit value defined by Donoho and Johnstone, which is known as the general threshold, which is calculated with . In this case, σ is the standard deviation of coefficients and n is the number of samples at each level [30].
3 Upscaling
In this paper, upscaling has been carried out using two methods: bandwidth and wavelet. The main purpose of the paper is simultaneous upscaling of the two features but first, the singlefeature upscaling is required. Therefore, Section 3.1 is presented for upscaling of a feature in the two methods, and then, in Section 3.2 simultaneous upscaling of the two properties of porosity and permeability is investigated. Finally, the results of the two methods are compared in Section 4 of the paper.
3.1 Singlefeature upscaling
3.1.1 Bandwidthbased upscaling
The bandwidth upscaling algorithm based on the distance matrix measures the distance of two adjacent cells in the x and y directions and stores them in matrices of (n × m − 1) in the x direction and (n − 1 × m) in the y direction. In this case, n and m are the amount of data in rows and columns, respectively. The purpose of upscaling is to combine cells with similar and interchangeable variability. The smallest element of the distance matrix is the starting point of the algorithm. Then, the distance matrix elements are compared with the optimal bandwidth vector. If the spacing of the two specified cells is less than the bandwidth, the two cells merge to one coarse cell. For each step of upscaling, the distance matrix will be computed again because the distance of adjacent cells to the coarsened cells will change. This process is an iterative and hierarchical process. The stopping condition of the algorithm is the bandwidth. If the difference between the two cells is greater than the bandwidth, the upscaling will not be performed. Using this algorithm, the upscaling initially starts from areas far from the reservoir heterogeneity and approaches in the next steps to heterogeneous areas.
The proposed scaleup algorithm is applied on SPE10 model data [2]. Porosity data has been used to singlefeature upscaling. The vector (0.1365, 0.1454) is the optimum bandwidth of the porosity data. In fact, two cells in the x direction can be magnified until their difference is not greater than 0.1365. The output parameters of the upscaling algorithm are the number of coarse cells and the upscaling error. We can control to prevent coarse cells from adapting to the small cells by adding a constraint to the size of the coarse scale cells. In fluid flow simulation, the dimensions of adjacent cells should not exceed 2–3 times of each other. Table 1 shows the upscaling parameters of the porosity based on the bandwidth method. Naturally, by selecting the smaller bandwidth, the scaleup model will be more similar to the original or fine model, and the larger bandwidth yields to an upscaled model in which more information from the fine scale model is lost. This can be seen in the error rate as well as the variance of the models in different bandwidths. As mentioned earlier, when the bandwidth increases the number of averaging between the property values of the cells increases, and hence, the number of cells in the upscaled model decreases (Tab. 1). Furthermore, as can be seen from Table 1, the variability and variance of the models will decrease by increasing bandwidth.
Upscaling parameters of porosity model.
Figure 4 shows the fine and coarse scale models for the SPE10 porosity model. Figure 4a is the fine scale porosity and Figure 4b shows the scaledup model based on the optimal bandwidth. Figures 4c and 4d illustrate the enlarged form of a specified part of the grid. These sections describe how to upscale the fine scale model.
Fig. 4 (a) Finescale porosity model, (b) upscaled model, (c) enlarged section of fine model and (d) enlarged section of upscaled model. 
3.1.2 Waveletbased upscaling
The purpose of upscaling is to create cells with specific properties for fluid flow simulation. Every upscaling method involves an averaging process; therefore, it is inevitable to remove some details in the coarser cells. On the other hand, if less averaging is done, minimum information from the fine scale model is lost. The main idea based on the wavelet transform is to increase the scale, reduce the number of averages, and thereby, preserve the basic and important information of the fine model. As mentioned earlier, the scaleup process using the wavelet transform approach is not uniform, and the grid will become coarse fully intelligent.
The basis of upscaling in wavelet transform is the calculation of approximation and detail coefficients. At each decomposition level, the grid with a factor of 2 will be upscaled in two directions. To implement the upscaling based on this method, the thresholds are determined for the coefficients. At each step, the cells will be coarse if the detail and approximation coefficients will be lower than the thresholds. Coefficients greater than the thresholds indicate the high importance of the cell and there is a reason for not to break and scale up. Areas that remain fine scale at current step will be fine at later stages. The key point in the waveletbased upscaling is to determine the thresholds for approximation and detail coefficients. By selecting the smaller thresholds, more information will be retained and the coarsening will start far from heterogeneities. The reason for choosing the approximation threshold is to maintain the fine structure of the grid in areas with high values of variables that must be finetuned to avoid loss of information.
The threshold can be applied in two ways of: (a) constant at each step and (b) variable and decreasing form. In order to better understand the effects of the desired feature distribution at different stages of coarsegraining, threshold values can be selected as a reduction at higher levels of decomposition. By decreasing the threshold at higher stages, the sudden increase in the cell size is controlled, and this can result in that the finescale cells not being placed next to very coarse cells. Thresholds are the fraction of the largest coefficients at special level. Figure 5 shows the result of upscaling based on the wavelet transform with decreasing threshold at four levels as {0.9, 0.8, 0.7, and 0.6}. The threshold concept means that, for example, in the first step, areas with approximation coefficients of less than 90% of the largest approximation coefficient at that level, and also, detail coefficients of less than 90% of the largest detail coefficient at that level, will be merged.
Fig. 5 Waveletbased upscaling: (a) Fine scale of porosity model, (b) upscaled model, (c) enlarged section of the upscaled model. 
As shown in Table 2, the fine scale data of the SPE10 porosity model after scaling up based on the wavelet transform to fourth level will convert to a model with 2881 coarse cells. The concept of decomposition up to level 4 is that in the process of upscaling, the maximum size of a coarse cell will be 16 times of the size of a fine cell. The upscaling parameters for the porosity model after scaling up based on the wavelet transform are presented in Table 2. The data variance can be a measure of the loss information in the fine scale model. In the scaleup model based on the wavelet constant threshold, the variance of the porosity model is reduced to 0.0015. In addition, the variance of model has been reduced to 0.0025 at the stepwise threshold but it gives better result than fixed threshold. The upscaling error in the two cases also confirms this result.
Upscaling parameters based on the wavelet method.
3.2 Simultaneous Upscaling of Two Features (SUTF)
3.2.1 Bandwidthbased upscaling
In this method, to determine the SUTF model, the porosity and permeability scaling models are directly compared. This method will have two different approaches: the minimum bandwidth approach and the maximum bandwidth approach. In the minimum bandwidth approach, the two models are compared and depending on how the cells are coarsened, the smaller bandwidth will be considered as the step criterion of scaling for the SUTF model. For example, suppose that along the x direction, 3 cells are merged in the permeability model and 8 cells are merged in the porosity model. Therefore, based on the minimum bandwidth, 3 cells will be considered as the bandwidth or dimension of a coarse cell in the x direction for the SUTF model. This process in the y direction will be done similarly. After determining the structure of the SUTF model, the porosity and permeability data will be calculated for the coarsescaled cells in the final model. Finally, an identical scaleup model will be calculated based on the minimum bandwidth approach for both features. The upscaling error of each feature for the SUTF model will be different from that of the single upscaled model due to decreasing the number of averages, and in particular, it is less than that of the single upscaled model. It is also obvious that the number of cells in the SUTF model is greater than the number of cells in the singlefeature upscaling models. The SUTF model is an identical model for the porosity and permeability features, and the difference between them is only the cell values, otherwise the dimensions and location of the coarser cells will be the same for both features. The results obtained for the SUTF model are shown in Table 3.
Upscaling parameters in minimum bandwidth approach.
The SUTF model has 5342 coarsescale cells for the porosity and permeability features. The SUTF model with these coarsened cells is the best model to show simultaneous variability of the two properties that most closely matches the finescale model. Figure 6 shows the results of the scale up and the SUTF model based on the minimum bandwidth approach.
Fig. 6 (a) SUTF model of porosity, (b) SUTF model of permeability, (c) enlarged section of porosity model, (d) enlarged section of permeability model, (e) data variability of section c, and (f) data variability of section d. 
As shown in the enlarged sections of the models (Figs. 6c and 6d), the upscaling has been done in a completely nonuniform way without any order. It is observed that the two models are quite similar and the only difference is in the type of property and the amount of cells. The cell distribution does not follow a certain order, and the upscaling has been carried out just based on variability. In the maximum bandwidth approach, the larger bandwidth will be the criterion for constructing the SUTF model. In this case, statistically in each area, the assembly of the two models will be selected as the coarse cell. It is therefore obvious that the number of coarsegrained cells in the SUTF model is lower than in the singlefeature models. Moreover, the upscaling error of the SUTF model will be greater than the errors of singlefeature models for the porosity and permeability properties. Table 4 shows the scaleup parameters of the SUTF model in the maximum bandwidth approach.
Upscaling parameters in maximum bandwidth approach
The SUTF model has 2448 coarse cells. It is clear that the number of cells in the SUTF model is lower than the singlefeature model, because, depending on the maximum bandwidth setting, some cells in the singlefeature model may be resized and merged again, which leads to an increase in the computational error. The upscaling error for porosity increases from 10.8 to 44.6 as well as it increases from 3.1 × 10^{7} to 9.09 × 10^{9} for permeability data. The results of upscaling on the maximum bandwidth approach are shown in Figure 7.
Fig. 7 (a) SUTF model of porosity, (b) SUTF model of permeability, (c) enlarged section of porosity model and (d) enlarged section of permeability model. 
In summary, the SUTF models that are obtained with maximum and minimum bandwidths have 5342 and 2448 coarsescale cells, respectively. The upscaling error of the SUTF model in the maximum bandwidth approach is far greater than that of the minimum bandwidth approach, as the number of the averaging operations is substantially increased. This is clearly seen in the enlarged sections of Figures 5 and 6. In a same area of 8 × 8 finescale cells, using the minimum bandwidth approach, 40 coarse cells are obtained as a result of upscaling, in which 62% of the fine scale structure is actually preserved. However, in the maximum bandwidth approach, in this area, the SUTF model has only 18 coarsescale cells, that indicate only 28% of the fine model is maintained in the upscaled model.
3.2.2 Waveletbased upscaling
Simultaneous Upscaling of Two Features is a new topic in the wavelet field, and so far, nothing has been done in this regard. Similar to the kernel function bandwidth method described above, it can also be used in wavelet transform, except that the scaleup model for the features is obtained on the basis of wavelet transform. The waveletbased upscaled models for porosity and permeability features have 2881 and 1218 coarsescale cells, respectively. The SUTF model with the same strategy as shown in the bandwidthbased upscaling method, with the specified inputs, is shown in Figure 8. Table 5 also indicates the scaleup parameters of different models. The SUTF model has 3005 coarsescale cells that comprise only 22% of the finescale model data. The upscaling error of the SUTF model for both features is less than that of the singlefeature upscaling because of increasing in the number of coarser cells of the SUTF model. However, the behavioral variance shows the opposite.
Fig. 8 Waveletbased SUTF, (a) upscaled model of porosity, (b) upscaled model of permeability, (c) SUTF model of porosity and (d) SUTF model of permeability. 
Waveletbased upscaling parameters in minimum bandwidth approach.
As can be seen from Figure 7, the coarse scale permeability feature model and its SUTF model are very different. This difference is also seen in the difference between Sum of Square Error (SSE) and data variance of models. In the case of permeability, the number of cells in the SUTF model is more than twice of that of the singlefeature model, as shown in Figure 7. In contrast, coarsescale models of the porosity in singlefeature and SUFT models are very similar. Generally, the number of cells in the SUTF models is approximately the same. However, the number of cells in the singlefeature upscaling models is different.
4 Comparison of bandwidthbased and waveletbased upscaling methods
In this section, the results of upscaling obtained by the two methods, i.e. kernel bandwidth and wavelet conversion, are compared. The purpose of comparing the results is to evaluate the performance of the bandwidth method against the wellknown wavelet method for upscaling. Based on what has been said, there is no idea about the number of cells in the scaleup model in the wavelet transform method, but in the kernelbandwidth method, it is possible to perform the scaleup process by determining the number of cells before starting the process. In the bandwidth method, since the upscaling is done based on the spacing matrix that refers to variability, there is no doubt that before the algorithm stops, all the cells that have the upscaling condition are enlarged.
Table 6 shows the upscaling parameters of the SPE10 porosity and permeability properties using two kernel bandwidth and wavelet transform methods. The scaleup model is based on the wavelet transform and has 2881 coarsescale cells. To compare this model with the result of upscaling based on the kernel function bandwidth, it is enough to select the bandwidth corresponding to this number of cells and calculates the scaleup model for the desired feature. The bandwidth upscaling algorithm is capable of terminating the scaleup process by applying a coarse cell number constraint and presenting the simulator model. In singlefeature upscaling, the kernel bandwidth error is less than the wavelet transform in both properties. The variance of data in the scaleup model based on the kernel bandwidth is about 3 times of the variance of the data in the scaleup model based on the wavelet transform. Therefore, it can be concluded that the kernel function bandwidth upscaling model is more similar to the finescale model and is able to retain the highest variability and variance of data with the same cell number as that of the wavelet transform method.
Comparison of the results of bandwidthbased and waveletbased upscaling methods.
The SUTF results obtained using different approaches can also be compared. The SUTF model obtained using minimum bandwidth approach and based on wavelet transform consists of 3005 coarse cells. As shown in Table 6, the calculated upscaling errors of SUTF model for porosity and permeability properties are 84 and 1.59 × 10^{10}, respectively. Furthermore, the simulator model based on the kernel function bandwidth method has 5342 coarsescale cells with upscaling errors of 41.96 and 6.1 × 10^{7} units for porosity and permeability features, respectively. It is clear that the upscaling error of kernel function bandwidth is lower than the wavelet transform method in both singlefeature upscaling and simultaneous two features upscaling modes. The reason for this difference that the upscaling algorithm is designed on the basis of data variability.
Figure 9 shows a comparison of part of the grid for the porosity and for the two upscaling methods. The waveletbased upscaling pattern is on a regular system but nonuniform. While the bandwidthbased upscaling pattern is very different from the waveletbased upscaling pattern; the results of both upscaling methods are seen in irregular and nonuniform form. Figures 9a and 9b illustrate the singlefeature upscaling and Figures 9c and 9d show the SUTF using both bandwidth and wavelet methods. Figures 9e and 9f also show the variability of the data for the SUTF models.
Fig. 9 (a) Bandwidthbased upscaling model, (b) waveletbased upscaling model, (c) SUTF model based on bandwidth method, (d) SUTF model based on wavelet method, (e) data variability of SUTF model based on bandwidth method and (f) data variability of SUTF model based on wavelet method. 
5 2D simulations
In order to considering the bandwidth upscaling method performance, we applied a water injection plan to simulating a heterogeneous reservoir. Then, simulation results compared in fine and coarse scale models. In this case, water flooding is done from one face. In this state, water enters from the west and oil leaves from the east face. The northern and southern borders are closed uncurrent. Permeability distribution is constructed as heterogeneous and layered (anisotropic). Reservoir parameters are shown in Table 7.
Reservoirs parameters in case study.
Simulation results are gathered in the following. For this case, upscaling results of wavelet transformation and bandwidth method compared with the fine geological model. Wavelet results are based on the decreasing thresholds and the bandwidth results calculated based on Section 3.1.1. Figures 10 and 11 show the comparison of the whole input water from the west border and the mean reservoir pressure. The performance of the two models is close to each other, but the bandwidth model, which is based on the distance matrix and identifies similarities, shows better and closer results than the wavelet transform method.
Fig. 10 Mean reservoir pressure variation based on the time in geological model and upscaled models. 
Fig. 11 Injected water flow variation based on the time in geological model and upscaled models. 
Figure 12 compared the total oil production of reservoir in different models. Despite the use of enlarge cells, both upscaling models show an acceptable trend for the total oil production of reservoir. As can be seen, over time and with decreasing production rate, the models get closer together and the upscaling error decreases. However, the bandwidth results have more correlation with the fine scale model and this can be justified according to the upscaling error calculated in Section 4 of the article. Figure 13 has investigated the average reservoir pressure in different models. Upscaling models have acceptable pressure behavior; however, the behavior of the bandwidth model is closer to that of the finescale model.
Fig. 12 Comparison of the reservoir oil production between upscaled models and fine model. 
Fig. 13 Comparison of the reservoir pressure between upscaled models and fine model. 
As mentioned earlier, increasing the speed of calculations as well as reducing the computational cost is one of the goals of upscaling of geological models. Although no direct simulation of the geological model has been performed in this study, but based on the number of cells, an estimate of the simulation time of the fine model can be obtained and compared with the computational time of the upscaled models. According to the estimation, the execution time of the model on a fine scale will be equal to 74 559 208 s [12]. Therefore, it is now possible to provide an estimate of the speed of the upscaling models shown in Table 8. As can be seen, the increase in speed in both methods of upscaling is very remarkable.
Comparison of simulation computation speed.
6 Conclusion
Upscaling of reservoir properties based on cell variability is the main idea of scaleup in the kernel function bandwidth method. Due to various degrees of variability in the cells, the upscaling using this method is made in nonuniform manner. Bandwidthbased upscaling is controlled by the bandwidth in different areas. The number of upscaled cells is a function of the heterogeneous level of the desired property. In areas with high heterogeneity, the selection of small bandwidths will preserve the basic information of the finescale model, and in areas with smooth or small variability, large bandwidths will reduce the number of cells. In this paper, in order to show the simultaneous variability of two features on a larger scale, we have employed SUTF method. The SUTF model can be calculated using two approaches: minimum bandwidth and maximum bandwidth. Each scaleup method can be evaluated by two important parameters, (a) accuracy of model and (b) computation time. If the purpose of upscaling is to maximize information retention, the minimum bandwidth approach for constructing the SUTF model has a better response. If computing time is the priority, the method of maximum bandwidth approach to build the SUTF model will have better results. The results of bandwidthbased upscaling method have been compared with the results of waveletbased upscaling method as a wellknown method in the upscaling with a multiresolution approach. This comparison confirms that under the same conditions, upscaling error based on the bandwidth method is far less than that of the wavelet method. In addition, the bandwidth method can perform the scale conversion process, depending on the number of cells as well as the expected error in the scaleup model. These criteria can be very effective in managing time and cost in fluid flow simulation calculations. By comparing the simulation results, it can be seen that the bandwidth method has much closer results than the wavelet method to the results of the fine scale geological model.
References
 Rios V.S., Santos L.O.S., Quadros F.B., Schiozer D.J. (2019) New upscaling technique for compositional reservoir simulations of miscible gas injection, J. Pet. Sci. Eng. 175, 389–406. [Google Scholar]
 Christie M.A. (1996) Upscaling for reservoir simulation, Soc. Pet. Eng. 48, 1–4. [Google Scholar]
 Durlofsky L.J. (2003) Upscaling of geocellular models for reservoir flow simulation: a review of recent progress, in: 7th International Forum on Reservoir Simulation, Bühl/BadenBaden, Germany, pp. 23–27. [Google Scholar]
 Durlofsky L.J. (2005) Upscaling and gridding of fine scale geological models for flow simulation, in: 8th International Forum on Reservoir Simulation, Stresa, Italy, June 20–24, 2005, pp. 1–10. [Google Scholar]
 McKee F., Preux C., Berthon C. (2012) Homogenization of relative permeabilities curves for twophase flow in porous media using an optimization method, in: ECMOR XIII, 13th European Conference on the Mathematics of Oil Recovery, Biarritz, France, Sep 10–13, 2012. [Google Scholar]
 Ma L., Dowey P.J., Rutter E., Taylor K.G., Lee P. D (2019) A novel upscaling procedure for characterizing heterogeneous shale porosity from nanometerto millimeterscale in 3D, Energy 181, 1285–1297. [Google Scholar]
 Liao Q., Lei G., Wei Z., Zhang D., Patil S. (2020) Efficient analytical upscaling method for elliptic equations in threedimensional heterogeneous anisotropic media, J. Hydrol. 583, 124560. [Google Scholar]
 Preux C., Le Ravalec M., Enchéry G. (2016) Selecting an appropriate upscaled reservoir model based on connectivity analysis, Oil Gas Sci. Technol. – Rev. IFP Energies nouvelles 71, 5, 60. [Google Scholar]
 Noetinger B., Zargar G. (2004) Multiscale description and upscaling of fluid flow in subsurface reservoirs, Oil Gas Sci. Technol.  Rev. IFP Energies nouvelles 59, 2, 119–139. [CrossRef] [Google Scholar]
 Panda M.N., Mosher C.C., Chopra A. (2000) Application of wavelet transforms to reservoirdata analysis and scaling, SPE J. 5, 92–101. [Google Scholar]
 Rasaei M.R., Sahimi M.R. (2008) Upscaling and simulation of waterflooding in heterogeneous reservoirs using wavelet transformations: Application to the SPE10 model, Transp. Porous Med. 72, 311–338. [Google Scholar]
 Rasaei M.R., Sahimi M. (2009) Upscaling of the permeability by multiscale wavelet transformations and simulation of multiphase flows in heterogeneous porous media, Comput. Geosci. 13, 187–214. [Google Scholar]
 Ebrahimi F., Sahimi M. (2004) Multiresolution wavelet upscaling of unstable miscible displacements in flow through heterogeneous porous media, Transp. Porous Med. 57, 75–102. [Google Scholar]
 Ebrahimi F., Sahimi M. (2006) Grid coarsening, simulation of transport processes in, and scaleup of heterogeneous media: application of multiresolution wavelet transformations, Mech. Mater. 38, 772–785. [Google Scholar]
 Marquez R.Y., Araujo M. (2006) Wavelet analysis for upscaling of two dimensional permeability fields, Física del Petróleo 52, 76–80. [Google Scholar]
 Pancaldi V., King P.R., Christensen K. (2007) Waveletbased upscaling of advection equations, Phys. A. 387, 4760–4770. [Google Scholar]
 Pancaldi V. (2008) Coarse graining equations for flow in porous media: a Haar wavelets and renormalization approach, PhD Thesis, Imperial College London, 189 p. [Google Scholar]
 Moslehi M., de Barros F.P.J., Ebrahimi F., Sahimi M. (2016) Upscaling of solute transport in disordered porous media by wavelet transformations, Adv. Water Res. 96, 180–189. [Google Scholar]
 Bouard L., Duval L., Preux C., Payan F., Antonii M. (2019) Refinableprecision in mesh compression for upscaling and upgridding in reservoir simulation with HexaShrink, in: 81st Conference and Exhibition EAGE Annual, London, 2019. [Google Scholar]
 Lie K.A., Kedia K., Skaflestad B., Wang X., Yang Y., Wu X.H., Hoda N. (2017) A general nonuniform coarsening and upscaling framework for reducedorder modeling, in: SPE Reservoir Simulation Conference, Texas, USA, pp. 1–15. [Google Scholar]
 Azad M.R., KamkarRouhani A., Tokhmechi B., Arashi M. (2019) Determination of optimal bandwidth in upscaling process of reservoir data using kernel function bandwidth, J. Min. Environ. 10, 3, 613–621. [Google Scholar]
 Azad M.R., KamkarRouhani A., Tokhmechi B., Arashi M. (2019) Simultaneous upscaling of two properties of reservoirs in one dimension using adaptive bandwidth in kernel function method, Pet. Explor. Dev. 46, 4, 756–762. [Google Scholar]
 Azad M.R., KamkarRouhani A., Tokhmechi B., Arashi M. (2019) Multiresolution simultaneous upscaling of two features of the reservoir using the bandwidth of the kernel function and wavelet transformation, J. Pet. Sci. Eng. https://doi.org/10.1016/j.petrol.2019.106697. [Google Scholar]
 Christie M.A., Blunt M.J. (2001) Tenth SPE comparative solution project: a comparison of upscaling techniques, Soc. Pet. Eng. 4, 1–10. [Google Scholar]
 Chen T., Gewecke N., Li Z., Aubiano A., Shuttleworth R., Yang B., Zhong X. (2009) Fast computational methods for reservoir flow models, 1–20. [Google Scholar]
 Silverman B.W. (1986) Density estimation for statistics and data analysis, Chapman & Hall/CRC. [Google Scholar]
 Simonoff J.S. (1996) Smoothing methods in statistics, in: Springer Series in Statistics, Springer, 340 p. [Google Scholar]
 Wang Sh, Li A., Wen K., Wu X. (2020) Robust kernels for kernel density estimation, Econ. Lett. 191, 109138. [Google Scholar]
 Raykar V.C., Duraiswami R., Zhao L.H. (2010) Fast computation of kernel estimators, J. Comput. Graph. Stat. 19, 205–220. [Google Scholar]
 Heidarinasab A., Dabir B., Sahimi M. (2004) Multiresolution waveletbased simulation of transport and photochemical reactions in the atmosphere, Atmos. Environ. 38, 37, 6381–6397. [Google Scholar]
All Tables
Comparison of the results of bandwidthbased and waveletbased upscaling methods.
All Figures
Fig. 1 (a) Porosity and (b) permeability fields for the SPE10 Top Layer [25]. 

In the text 
Fig. 2 Twodimensional upscaling algorithm based on the kernel bandwidth. 

In the text 
Fig. 3 (a) Wavelet decomposition and (b) 1D decomposition. 

In the text 
Fig. 4 (a) Finescale porosity model, (b) upscaled model, (c) enlarged section of fine model and (d) enlarged section of upscaled model. 

In the text 
Fig. 5 Waveletbased upscaling: (a) Fine scale of porosity model, (b) upscaled model, (c) enlarged section of the upscaled model. 

In the text 
Fig. 6 (a) SUTF model of porosity, (b) SUTF model of permeability, (c) enlarged section of porosity model, (d) enlarged section of permeability model, (e) data variability of section c, and (f) data variability of section d. 

In the text 
Fig. 7 (a) SUTF model of porosity, (b) SUTF model of permeability, (c) enlarged section of porosity model and (d) enlarged section of permeability model. 

In the text 
Fig. 8 Waveletbased SUTF, (a) upscaled model of porosity, (b) upscaled model of permeability, (c) SUTF model of porosity and (d) SUTF model of permeability. 

In the text 
Fig. 9 (a) Bandwidthbased upscaling model, (b) waveletbased upscaling model, (c) SUTF model based on bandwidth method, (d) SUTF model based on wavelet method, (e) data variability of SUTF model based on bandwidth method and (f) data variability of SUTF model based on wavelet method. 

In the text 
Fig. 10 Mean reservoir pressure variation based on the time in geological model and upscaled models. 

In the text 
Fig. 11 Injected water flow variation based on the time in geological model and upscaled models. 

In the text 
Fig. 12 Comparison of the reservoir oil production between upscaled models and fine model. 

In the text 
Fig. 13 Comparison of the reservoir pressure between upscaled models and fine model. 

In the text 