Upscaling the porosity–permeability relationship of a microporous carbonate for Darcy-scale flow with machine learning

The permeability of a pore structure is typically described by stochastic representations of its geometrical attributes (e.g. pore-size distribution, porosity, coordination number). Database-driven numerical solvers for large model domains can only accurately predict large-scale flow behavior when they incorporate upscaled descriptions of that structure. The upscaling is particularly challenging for rocks with multimodal porosity structures such as carbonates, where several different type of structures (e.g. micro-porosity, cavities, fractures) are interacting. It is the connectivity both within and between these fundamentally different structures that ultimately controls the porosity–permeability relationship at the larger length scales. Recent advances in machine learning techniques combined with both numerical modelling and informed structural analysis have allowed us to probe the relationship between structure and permeability much more deeply. We have used this integrated approach to tackle the challenge of upscaling multimodal and multiscale porous media. We present a novel method for upscaling multimodal porosity–permeability relationships using machine learning based multivariate structural regression. A micro-CT image of Estaillades limestone was divided into small 603 and 1203 sub-volumes and permeability was computed using the Darcy–Brinkman–Stokes (DBS) model. The microporosity–porosity–permeability relationship from Menke et al. (Earth Arxiv, https://doi.org/10.31223/osf.io/ubg6p, 2019) was used to assign permeability values to the cells containing microporosity. Structural attributes (porosity, phase connectivity, volume fraction, etc.) of each sub-volume were extracted using image analysis tools and then regressed against the solved DBS permeability using an Extra-Trees regression model to derive an upscaled porosity–permeability relationship. Ten test cases of 3603 voxels were then modeled using Darcy-scale flow with this machine learning predicted upscaled porosity–permeability relationship and benchmarked against full DBS simulations, a numerically upscaled Darcy flow model, and a Kozeny–Carman model. All numerical simulations were performed using GeoChemFoam, our in-house open source pore-scale simulator based on OpenFOAM. We found good agreement between the full DBS simulations and both the numerical and machine learning upscaled models, with the machine learning model being 80 times less computationally expensive. The Kozeny–Carman model was a poor predictor of upscaled permeability in all cases.


Methods
Creating the training dataset. The porosity-permeability curve for the micro porosity is provided by the previous work in Menke et al. 32 . A core of Estaillades limestone was scanned in a micro-CT both dry and saturated with a high contrast brine. The differential image was used to estimate the porosity of the connected micro porosity. Microporous subsections of the core were then scanned in a nano-CT and the grain size distribution of the micro porosity was modelled numerically and the results used to generate a synthetic porosity-permeability for the micro porosity. This curve was then used in a Stokes-Brinkman simulation of the whole core and benchmarked against experimental permeability measurements of the core with high accuracy. A detailed discussion of this multi-scale imaging and benchmark modelling can be found in Menke et al. 32 . The image of Estaillades that is used in this study is the same as the image originally taken in Menke et al. 32 and due to the extensive characterization performed in this study for the purposes of this work we will thus assume the permeability results and porosity-permeability curve derived in Menke et al. 32 to be ground truth. The raw and processed micro and nano CT images are all publicly available open access on the image archive of the British Geological Survey 44 .
We used the 15-phase (pore, solid grain, viton sleeve and 12 phases of microporosity) segmented micro-CT image of Estaillades which is 1200 × 1200 × 6000 voxels with a resolution of 3.9 μm. The top 10% of the image was split into two training datasets: one with 30,000 sub-volumes of 60 3 voxels with a 50% overlap between sequential sub-volumes to increase the number of training images and one with 30,000 sub-volumes of 120 3 voxels with a 75% overlap between sequential sub-volumes. Overlapping the datasets allowed us to keep the two training datasets directly comparable. www.nature.com/scientificreports/ Structural features were then extracted from each of the sub-volumes using the image analysis toolbox in python's SciKitImage 45 . This analysis included the volume fractions of each phase in the sub-volume as well as the cumulative connectivity of the phases in each orthogonal direction expressed mathematically as the first connected phase. An example of the calculation of the connectivity of the primary porosity of a sub-volume is shown in Fig. 1. This image analysis created a robust feature vector set of 18 features: 15 volume fractions and 3 connectivity features.
Calculation of permeability using GeoChemFoam. The training datasets (1) and (2) were solved for permeability in each direction using the DBS approach 40 in which one equation is used to model the flow within the fully resolved pores (i.e. voxel porosity equal to 1.0) and the micropores (i.e. voxel porosity lower than 1.0).
where u [m·s −1 ] is the fluid velocity, p [Pa] is the pressure, ϕ [−] is the porosity, µ kg m −1 s −1 is the viscosity and K [m 2 ] is the permeability. Porosity and permeability of each microporous voxel is assigned using the segmented phases and values from Menke et al. 32 (Table 1).
For the upscaling models, Eq. (1) Can be simplified into the Darcy equations.
Here K is an anisotropic diagonal tensor that represents the permeability of each sub-volume (60 3 or 120 3 voxels). Each diagonal coefficient ( K x , K y and K z ) has been computed using one of the three upscaling models presented in "Upscaling to the Darcy-scale" section. Each sub-volume is modelled using a 4 × 4 × 4 grid with constant coefficients, so that the total grid for the Darcy simulation is 24 × 24 × 24 when using the 60 3 subvolumes and 12 × 12 × 12 when using the 120 3 sub-volumes.

Results and discussion
The results and discussion are organized into four parts. (1) In "Sub-volume permeability, Kozeny-Carman model, and multivariate regression" section we look at the relationship between porosity and permeability of the sub-volumes as solved by DBS and compare it to both, a power law fitted Kozeny-Carman model and the predicted permeabilities of our machine learning multivariate regression model. (2) In "Analysis of feature importance" section we examine the features used in the machine learning regression model and their importance relative to the number of features used in the regression as well as the error associated with the number of training images used to train the regression model. The feature vectors obtained from 29,000 of the 30,000 sub-volumes in training dataset 1 were then used to train an extra randomized trees ensemble machine learning regression model using SciKitLearn 41 . The remaining 1000 sub-volumes were then used to test the regression model performance. To avoid model overfitting, bootstrapping was used during tree creation where a randomly selected subset of models was used to train each tree, and the number of trees was limited to 50. The average out-of-bag error calculated on each tree with unused training images was always within 10% of the average R 2 score of the model calculated with the images used in each tree, indicating that the models were not overfitted.
The porosity values from the training sub-volumes were then plugged into the Kozeny-Carman model (Fig. 2B). The machine learning regression model significantly outperformed the Kozeny-Carman model where the RMSE of the Kozeny-Carman permeability predictions was 29.7% while the machine learning regression model had a RMSE of 4.3%. Furthermore, the Kozeny-Carman model's predicted values clustered towards the mean and the model was unable to predict the highest and lowest permeability values. This clustering indicates that the Kozeny-Carman model lacks enough complexity to incorporate the multiscale heterogeneities inherent in carbonate rocks with accuracy.
Analysis of feature importance. Three test cases of feature vectors, comprising 3, 15, and 18 features respectively, were assembled. In the '3 features' test set, only the connectivity feature vectors were included. In www.nature.com/scientificreports/ the '15 features' set, only the porosity volume fraction features vectors were considered, while in the '18 features' set, both the connectivity information and the porosity volume fraction features were included in the model. Figure 3A shows the RMSE of the different feature sets with the number of training sub-volumes used to train  www.nature.com/scientificreports/ the model. All models showed a sharp drop in RMSE when increasing the number of sub-volumes used in the regression followed by diminishing returns as the number of sub-volumes increased past 500. This trend indicates that for this particular size of the sub-volume in this rock, very few sub-volumes are required to train an accurate model. However, this assertion is rock dependent, and we expect this number to change significantly with other rock types. Additionally, we found that the RMSE of the '18 features' case was the best with a minimum RMSE of 4.3%. However, this is only a slight improvement over the '3 features' set with a minimum RMSE of 5.6%, but a large improvement over the '15 features' set that has a minimum RMSE of 15.3%. This trend indicates that, as expected, the connectivity features are more important in predicting permeability than porosity alone. Figure 3B shows the relative model weighting of each feature. Here we see that the connectivity features are weighted an order of magnitude more highly than the porosity features in all cases, confirming that when predicting permeability, the connectivity is a better predictor than porosity. This observation also highlights the importance of extracting the most informative features during the image processing and analysis. It would be interesting to also see if the connectivity could be quantified more robustly in the primary porosity (phase 1) to increase the accuracy of prediction in highly connected porosity sub-volumes. However, this analysis is out of the scope of this study and is a target of future work.
Analysis of the size of the training data. As the ultimate goal of this work is to upscale to the Darcyscale, the choice of sub-volume size must be carefully considered. The sub-volumes must both be of a size relevant to Darcy-scale imaging and modelling, but also be below the REV scale for porosity and permeability so that the spread of permeabilities is sufficient to capture structural heterogeneity while still containing enough  www.nature.com/scientificreports/ tractable information on both the nano and micron-scale porosity structures to make the feature selection in a multivariate regression representative of the characteristic properties influencing flow. The image resolution of a medical CT scanner typically ranges between 250 and 500 microns, which corresponds approximately to the sizes of our sub-volumes chosen for investigation, i.e. 60 3 voxels and 120 3 voxels. To keep the total volume of rock in the training set constant (so as not to introduce any additional heterogeneity), the 120 3 sub-volumes overlapped by 75% instead of the 50% in the 60 3 sub-volumes. The same feature vector extraction workflow was used for both training datasets with between 100 and 29,000 training images and a further 1,000 test sub-volumes. Figure 4A shows the RMSE of the regression model predictions against the number of sub-volumes used in the training for both the 60 3 and 120 3 cases. Both cases converge to approximately the same RMSE of ~ 4% with 5000 training images. However, from Fig. 4B it is apparent that the spread of permeabilities is much greater for the 60 3 case, indicating that the 60 3 sub-volume size might be a better size to characterize the rock more accurately. This assertion is investigated further in "Upscaling to the Darcy-scale" section.
Upscaling to the Darcy-scale. In this section we compare the ground truth DBS solution in ten 360 3 test cases (Model 1) to the permeability predictions of three different upscaling models (Models 2, 3 & 4). First, we divided the remaining 90% of the micro-CT image not used in regression model training or benchmarking into ten 360 3 voxel test cases which were solved with DBS (Model 1). The test cases were then further subdivided into 6 × 6 × 6 and 3 × 3 × 3 matrices of sub-volumes of 60 3 and 120 3 voxels respectively for the Darcy upscaling tests (Fig. 5) These sub-volumes were then solved with DBS and the output permeability used in each block of the Darcy model (Model 2). Next, the feature vectors were extracted for each of the sub-volumes and the trained regression model was used to predict the permeability of each sub-volume using the feature set. Each subvolume was then assigned the permeability for the Darcy-scale flow model with its calculated total porosity from the feature vector set (Model 3). Finally, the porosity of each sub-volume was used to predict permeability using the Kozeny-Carman model fit described in "Sub-volume permeability, Kozeny-Carman model, and multivariate regression" section and that was used for each grid block of the Darcy-scale model (Model 4).  The permeability results of all these simulations are shown in Fig. 6 and Table 2. We found that in the case of the 60 3 sub-volumes the numerical upscaling and the machine learning upscaling both did equally well in predicting permeability with RMSE errors of 0.10 and 0.14, respectively, while the Kozeny-Carman did very poorly with a RMSE of 0.44. This poor performance is especially apparent in the case of test volume 2 which has the lowest porosity but the highest overall permeability in the ground truth, which is not something the Kozeny-Carman fit can predict. The upscaling results were not as accurate for the 120 3 sub-volume cases with the machine learning upscaling slightly outperforming the numerical upscaling with a RMSE of 0.16 over the numerical upscaling's 0.11. The Kozeny-Carman predictions were the worst, with a RMSE of 0.58. The overall increase in error for the 120 3 volumes can be explained by the smaller spread of sub-volume permeabilities shown in Fig. 4B as compared to the 60 3 sub-volumes. The 120 3 sub-volumes are too big to characterize the range of flow heterogeneities in the rock at this resolution and thus the output permeabilities all trend towards the mean permeability. It is also important to note that in all cases the total CPU time for the machine learning model after model training was as low as the Kozeny-Carman model but had similar prediction performance to both, the full DBS case and numerical upscaling.

Conclusion
In this study we have used the DBS model in GeoChemFoam in combination with decision tree based multivariate regression to upscale a microporous carbonate from the pore scale to the Darcy-scale. We found that multivariate regression can be used to upscale and predict permeability with very few training images and performed equally well to numerical upscaling with a fraction of the computational cost. Additionally, increased sub-volume size had little effect on model predictions. However, a bigger sub-volume meant increased model CPU cost and decreased the accuracy of the upscaling models. Furthermore, we found that appropriate choice of feature vectors for extraction is important for regression model performance and connectivity information  www.nature.com/scientificreports/ is the most important feature to include in the models. Machine learning based multivariate regression is thus an effective way of increasing prediction speed and accuracy during upscaling, but this method requires both precise multiscale imaging and an in depth understanding of connectivity in multiscale porosity structures. It is important to note that for the purposes of this study, we have defined the Darcy-scale to be any scale where the porosity-permeability relationship is explicitly defined for every grid block in the model. However, all samples are derived from a single core plug that is reasonably homogenous for a carbonate, and thus do not include larger features such as vugs and fractures that would likely be present in core plugs from a highlyheterogenous carbonate reservoir. It would be relatively straightforward to train a model on a different rock sample of similar size and complexity using the same imaging and modelling procedures. However, larger porosity structures would require an additional step in the workflow. They could be incorporated into the modelling framework using another upscaling step with the same D-B-S approach, where Darcy flow is computed in the porous matrix and Stokes in the fractures and vugs.
Reservoir characterization and modelling as well as history matching and uncertainty quantification at scales larger than a core plug require many core samples, which are expensive to drill and image. One advantage of this upscaling technique is that realistic synthetic geometries could be created using already established machine learning tools 50 with a small number of samples to ultimately create a large image library with accurate porosity-permeability relationships for use in uncertainty quantification without the need for large numbers of drilled samples or laboratory experiments. The addition of these larger structural scales is out of the scope of this work, but is a target of future research.
Furthermore, using this workflow to calculate permeability is just a steppingstone toward investigating more complex physics, in particular multiphase and reactive transport. Future work will also include investigation into using features that do not require high resolution imaging such as Darcy-scale porosity and tracer transport curves.