Artificial neural network modeling of scale-dependent dynamic capillary pressure effects in two-phase flow in porous media

Characterizing single and two-phase flow in porous media is of particular interest in issues relating to remediation of contaminant, oil recovery, flow in pulp and paper systems, geological sequestration of CO₂, and others (Ingham & Pop 2005; Hill et al. 2007; Kobayashi et al. 2008; Bear 2013; Abidoye & Das 2015; Das et al. 2014). The topic of the paper, namely, the two-phase flow behavior in porous medium can be macroscopically described by the equations for conservation of fluids’ mass and momentum (Kalaydjian 1987); however, these should be coupled with constitutive relationships among capillary pressure (P^c) saturation (S) and relative permeability (K_r) (Abidoye & Das 2015; Khudaida & Das 2014). The traditional approaches for determining the P^c–S–K_r relationships assume that the flow parameters are functions of steady-state saturation, i.e., when dS/dt = 0. For example, the definition of P^c according to the Laplace law is valid under static conditions (Kalaydjian 1992) where is the interfacial tension between the two fluids and is the contact angle. As evident, the law defines P^c to be a function of contact angle which, in turn, depends on the physical properties of the fluids and the porous medium (e.g., wettability (Dullien et al. 1990) and viscosity ratio (Danish & Jacquin 1983)) that are in contact with one another. Despite these notions, the saturation-rate (dS/dt) dependency of the P^c–S relationships has been observed in many studies during dynamic two-phase flow in porous media which is referred to as the dynamic capillary pressure effect (Hassanizadeh et al. 2002; Mirzaei & Das 2007).

As explained by Das et al. (2007), the magnitude of relates to how close or far from capillary equilibrium the two-phase flow system is. However, its magnitude is reported to be dependent on the size of the domain (Dahle et al. 2005; Camps-Roach et al. 2010; Bottero et al. 2011a, 2011b; Das & Mirzaei 2012; Das & Mirzaei 2013) apart from other factors, e.g., fluid viscosity and density ratios (Gielen et al. 2005; Das et al. 2007; Goel & O'Carroll 2011; Joekar-Niasar & Hassanizadeh 2011), permeability of the medium (Camps-Roach et al. 2010; Tian et al. 2012; Hanspal et al. 2013) and heterogeneities (Manthey et al. 2005; Das et al. 2006; Das & Mirzaei 2013; Mirzaei & Das 2013). This paper is specifically concerned with the domain scale effects on the dynamic coefficient as discussed below.

Understanding the influence of domain scale on the magnitude of is important as the two-phase flow can occur in porous domains at pore, core and/or larger field scales. Previous experimental (Bottero et al. 2011b) and numerical (Dahle et al. 2005) studies showed an increasing magnitude of as the domain size increases. However, some authors have expressed different conclusions about the effect of domain size on ⁠. Camps-Roach et al. (2010) and Das & Mirzaei (2012, 2013) do not find significant difference between the locally measured and the upscaled value of ⁠. This indicates the possible inconsistency surrounding the exact trend of the magnitude of with the domain scale.

Meanwhile, investigating the dynamic capillary pressure effects on two-phase flow systems in porous media often requires time-consuming experiments and/or cost-intensive modeling and simulations, which generally involve complex procedures to set up and run the simulations (Spalding 1981; Hanspal et al. 2013; Khudaida & Das 2014). Furthermore, determining the effects of domain scale on also imposes further challenges with the implications of various averaging techniques which are proposed and applied in the literature. For example, Bottero et al. (2011a) proposed the centroid-corrected averaging method as the most appropriate for two-phase systems while authors such as Das & Mirzaei (2012, 2013), Camps-Roach et al. (2010) and Manthey et al. (2005) employed different averaging techniques to determine the scale dependency of ⁠. While the scale dependency of continues to be a topic of discussion, there is an obvious lack of industrially relevant, easy to use, tools that can easily determine these behaviors. To address these challenges/issues, we wish to investigate the domain scale effect on using alternative platforms offering less complex and fast implementation procedures to achieve the same end. As such, a computationally economical and simpler platform for investigating the effects of domain scale on is presented using an artificial neural network (ANN).

Indeed, the ANN is a powerful modeling tool with the ability to learn and generalize functions from rounds of training as well as extract essential information from data (Wang & Fu 2008; Khashei & Bijari 2013). It provides a novel, elegant and valuable class of computational tools for data analysis and prediction (White 1989; Deka & Quddus 2014). Its building blocks or the elements are the ‘neurons’ which are grouped into input, hidden and output layers with respective biases, weights and transfer functions (Mueller & Hemond 2013; Yurdakul & Akdas 2013). The network manipulates the values of the biases and weights in a sequence of training processes and uses the transfer functions to establish the relationships between the inputs and the outputs. It has found applications in wide areas of science and engineering problems including medical fields to illustrate medical diagnosis (Amato et al. 2013), renewable energy systems, economics, psychology and many more (Kalogirou 2000).

Although the modeling parameters involved in the two-phase flow in porous media are interrelated in a complex manner, it has been shown that an ANN can approximate the relevant functions to the desired accuracy (Zhang et al. 1998; Hanspal et al. 2013; Das et al. 2015). This quality of the ANN can, thus, be harnessed to investigate the complex behavior of two-phase flow in porous media. For example, please see the works on the application of the ANN to study two-phase flow pattern (Mehta et al. 2013), oil flow rate (Ahmadi et al. 2013), groundwater contamination and pollutant infiltration forecasting (El Tabach et al. 2007), optimization of groundwater remediation problems (Rogers & Dowla 1994; Johnson & Rogers 2000), large-scale water resource management (Yan & Minsker 2006; Mounce et al. 2014), and permeability modeling in petroleum reservoir management (Karimpouli et al. 2010). Recently, Hanspal et al. (2013) demonstrated the effectiveness of an ANN in the determination of dynamic effects in a two-phase flow system, though they did not investigate the effect of domain scale on ⁠. They concluded that a well-trained and validated ANN structure can give reliable prediction of in a two-phase flow system. In addition to being inexpensive, an ANN offers a faster alternative to modeling of a complex system with the freedom from excessive imposition of constraints on the complex relationships between the input and output variables. Thus, our work explores the above qualities of the ANN to investigate the domain scale dependency of in the two-phase flow system in porous media.

ANN and multivariate regression (MVR) techniques were used to investigate the domain scale dependency of the in the two-phase flow system. The MVR was chosen because of the ease of implementation and to provide an alternative approach to compare different ANN configurations against.

For successful modeling of dynamic two-phase flow behavior using an ANN, the impacts of the network configuration, training and evaluation procedures cannot be overemphasized. In this work, different network configurations were investigated using a feed forward network. This is the commonest network in engineering application (Kalogirou 2000; Deka & Quddus 2014). For the purpose of training, a back-propagation algorithm was employed. Details of the configurations, training and data processing are discussed below.

In this paper, the literature data were obtained from the results by Das et al. (2007), Mirzaei & Das (2007), Goel & O'Carroll (2011), Das & Mirzaei (2012) and Hanspal & Das (2012) and additional data (Abidoye & Das 2015) from in-house laboratory experiments using methodology described by Das & Mirzaei (2012). Nine independent variables that have been identified as important in the literature were used as input variables in the simulations. These include, water saturation (S), media permeability (k), capillary entry pressure (P_d), fluid viscosity ratio (μ_r) defined as the ratio of the non-wetting phase viscosity (μ_nw) to that of the wetting phase viscosity (μ_w), fluid density ratio (D_r) defined as the ratio of the non-wetting phase density (D_nw) to that of the wetting phase density (D_w), temperature (T), pore size distribution index (λ), porosity (ɸ), and domain volume (V). The output variable is the corresponding ⁠. The number of data points under each variable is 307. In selecting these data sources, efforts were made to ensure that the data contain varying experimental or simulation parameters and conditions needed to fulfill the objective of this paper. For example, the works of Mirzaei & Das (2007), Goel & O'Carroll (2011), Das & Mirzaei (2012) and Hanspal & Das (2012) were conducted using different domain volumes. Furthermore, our laboratory experiments (Abidoye & Das 2015) were conducted for different heights of domain (4, 8 and 12 cm height) but with the same diameter for 500 cSt viscosity ratio of silicone–oil water system. These features enhance the training of the ANN to easily capture the nonlinear relationships between the domain size and ⁠. Important statistics of the variables are listed in Table 1.

Compared to the number of data points commonly required to plot a complete curve (typically less than 10 data points, see, e.g., Bottero et al. (2011b)), the amount of data used in this work (>300 data points) is over 30 times more than what is typically required. Also, this work employs a simple ANN structure. These features ensure that artificial over-fitting of the data is avoided. Complex ANN structures and few data can lead to artificial over-fitting in ANN modeling (Hanspal et al. 2013).

Various configurations of ANN were developed and tested to determine the most suitable network. The configuration approach followed that demonstrated in Hanspal et al. (2013) as proposed by Srinivasulu & Jain (2006). The ANNs include single and double hidden layers. A program file with lines of code was written and implemented in MATLAB to create, train, validate and test the networks as well as to generate the goodness of fit of the parameters, e.g. correlation coefficients and slope for the predicted output ⁠. The developed networks consist of different layers comprising the input, hidden and the output layers. The input layer is occupied by the independent variables while the output layer is for the dependent variable. The hidden layer is occupied by the neurons which are the constitutive units that receive the input and operate on them to produce the output. The code divides the data set randomly into 60, 20 and 20% corresponding to the data for training, validation and testing, respectively. As stated before, the training was performed with the Levenberg–Marquardt function (Marquardt 1963) using a back-propagation algorithm. The Levenberg–Marquardt function is a curve-fitting function applied in nonlinear least squares problems. It optimizes the parameter of the model curve by minimizing the sum of the squares of the deviation from the empirical dependent variable. The back-propagation learning algorithm operates by iterative adjustment of the weights and biases in response to the error value between the predicted and the desired outputs. ‘Tansig’ and ‘Purelin’ transfer functions were used in this work. These transfer functions calculate a layer's output from its net input. While ‘Tansig’ is nonlinear, ‘Purelin’ is linear. For a network with a single hidden layer, ‘Tansig’ was used between the input and the hidden layers while ‘Purelin’ was used between the hidden and the output layer. For a network with a double hidden layer, ‘Tansig’ was used between the input and the hidden layers as well as between the first and the second hidden layers while ‘Purelin’ was used between the second hidden layer and the output.

Mean square error (MSE) was employed as the network default performance criterion relating the calculated outputs from the ANN to the actual target (dependent variable) in the training, validation and testing processes. In the simulation, pre-processing was performed using lines of code in the script with function ‘mapminmax’. This function scales the inputs so that they fall into the range of −1 to 1.

In the training process, the epochs and goals serve as the stopping criteria of the number of iterations and the error tolerance, respectively. Epoch is the maximum number of times all of the training sets presented to the network while goal refers to the maximum error tolerance to be met by the developed network. Thus, the training stops if the error goal is met or the maximum number of epochs is attained. In this work, an epoch of 200 and a goal of zero were used. Different network configurations were constructed and each configuration differs in the number of hidden layers or neurons. The number of neurons was gradually increased for either single or double hidden layers. In this work, the representation of the layers in the ANN configurations is ANN [X-H1-Y] and ANN [X-H1-H2-Y] for single and double hidden layers, respectively. ‘X’ represents the input layer and its number refers to the number of independent variables, ‘H1’ and ‘H2’ represent the first and the second hidden layers, respectively, and their numbers represent the number of neurons in that layer. ‘Y’ is the output layer and its number represents the number of the dependent variable.

For the purpose of comparisons with the performances of the different ANN configurations, multiple linear (LR) and nonlinear (NLR) regression models were investigated with the aid of MATLAB. Both regression models utilized the entire data set.

The performance of all ANNs were weighed with different statistical evaluations as demonstrated in Hanspal et al. (2013) using the following statistical analyses.

Following the rigorous statistical evaluation of the models developed and described above, prediction of the effect of domain size on the τ−S relationships was performed using the best-performing model. Separate data from independent experiments are predicted separately. To do this, the actual domain volume (V) of the experiment was increased by 10 or 20% and the corresponding was predicted as a function of saturation using the best-performing ANN.

A computationally cost-effective and reliable ANN structure that predicts the domain scale dependency of relationships will serve a useful purpose in determining the significance of dynamic capillary pressure for two-phase flow systems. Our results aim to demonstrate this possibility. So, the results of training, validation and testing of the ANNs are discussed below. Also, the performances of the different ANN configurations together with the MVR models are compared on the bases of the different performance criteria. Scale dependency of relationships were then predicted for two-phase flow systems using the best-performing configuration.

The training, validation and testing as well as the post-training regression analyses for different ANN configurations are shown in Figures 1 and 2 for ANN [9-13-15-1] and ANN [9-15-17-1], respectively. Figure 1(a) shows how the MSE reduces during training, validation and testing as the number of epoch increases. This eventually culminates in the optimal performance during validation at 108 epochs having approximately zero MSE value, i.e., 0.0023. This behavior shows that the network learns well as the number of epochs increase. The testing session shows acceptable MSE that is very close to zero as well. The post-training regression analysis (Figure 1(b)) shows how the LR line fits the data points. This regression line has a correlation coefficient (c) and slope (m) of 0.99 and 0.97, respectively, which are very close to 1. These show the reliability of the fit. In the figure, it can be observed that the target data cluster around the regression line in a way that shows reliable prediction. Similarly, the behavior of the network for ANN [9-15-17-1] is shown in Figure 2. The network exhibits similar behavior as discussed above. The learning improves with the number of epochs for the training, validation and testing (Figure 2(a)). This is indicated by the reduction in the MSE values and the optimal MSE with validation occurs before the number of epochs reaches 51. In comparison to the behavior of ANN [9-13-15-1], shown in Figure 1(a), the testing and validation errors are larger in Figure 2(a), under the corresponding condition. The regression line for the ANN structure [9-15-17-1] is shown in Figure 2(b). The fit shows good c and m values of approximately 0.95 and 1.00, respectively. The regression analysis in Figure 2(b) shows that the cluster of the target data around the regression fit line is more scattered unlike that shown for ANN [9-13-15-1] (Figure 1(b)). This gives indication that ANN [9-13-15-1] may be more reliable than ANN [9-15-17-1]. However, the latter learns faster than the former.

The performances of various ANN configurations, in terms of the values of c and m of the regression line of fits to the target data are listed in Table 2. As is well known, values of the slope m and correlation coefficient c closer to 1 indicate reliable prediction. From the table, it can be observed that the majority of the ANN configurations perform well as indicated by the high values of c and m. However, the two-hidden-layer models perform better than the single-hidden-layer models as they have the slope and correlation coefficient closer to 1 than the single-hidden-layer structure. In all, ANN [9-13-15-1] and ANN [9-15-17-1] appear to be leading in performance. However, the criteria listed in the ‘ANN performance testing criteria’ subsection are further employed in weighing the performance of all the models including the MVR models. In Table 2, the slope of the regression line obtained from ANN [9-15-17-1] is shown to be slightly greater than 1 (i.e., 1.02). This can be explained to mean slight overprediction of the target data by the model. This is also visible in Figure 2(b) where the line of fit is slightly above the best line of fit (i.e., Y = T). In contrast, Figure 1 shows that the regression line obtained from ANN [9-13-15-1] lies slightly below the best line of fit (i.e., Y = T) and hence, the slope is slightly below 1 (i.e., 0.97).

Figure 3 shows the plots of AARE for all ANN configurations, LR and NLR models. In comparison, AARE is generally low for the different ANN configurations while the LR and NLR generally have high AARE. The lower the AARE the better the performance (Hanspal et al. 2013). Thus, it seems that the ANNs perform better than the LR and NLR models. Among the ANN configurations, ANN [9-13-15-1] has the least AARE followed by ANN [9-11-13-1].

In Figure 4, the plots of the SSE similarly show that the SSE is generally higher for LR and NLR models. For the ANN structures, the ANN [9-13-15-1] configuration has the least SSE followed by ANN [9-11-13-1]. Comparison of the model output in relation to the target is described in terms of Nash-Sutcliffe efficiency coefficient (E) depicted in Figure 5 for all the models. Again, it is visible that ANN [9-13-15-1] has the highest efficiency followed by ANN [9-11-13-1]. Also, TS for all the models show that ANN [9-13-15-1], ANN [9-11-13-1] and ANN [9-15-17-1] have the leading percentages. The plots are shown in Figure 6 for TS 5, TS 10, TS 25, TS 50 and TS 100. High values of TS imply good model performance.

From the above discussions, the performance criteria show that ANN structures have better reliability in predicting the two-phase flow parameters than MVR models (both linear and nonlinear). From the results, ANN [9-13-15-1] has shown the best performance.

Results of the statistical analyses discussed in the ‘ANN configurations’ section compare the predicted output of the models to the actual target output ⁠. The discussions in this section focus on the comparison of the actual relationships (target) to that predicted by some selected ANN, which include the best-performing model: ANN [9-13-15-1]. Figure 7 shows the plot of relationships using the entire data set in comparison with the prediction by ANN [9-13-15-1]. One can observe that the ANN structure has a good predictive ability of the relationships at both low and high saturation. Almost the entire data set are overlaid by the predicted values. In the figure, is shown to increase as the saturation reduces. From the start of the displacement of the wetting phase by the non-wetting phase, there exists only minimal change in as the saturation reduces. The trend, however, changes around the irreducible saturation where the value rises very steeply. This trend is widely reported in literature (Camps-Roach et al. 2010; Sakaki et al. 2010; Goel & O'Carroll 2011; Das & Mirzaei 2013). According to Das et al. (2007), increase in the magnitude of indicates increased deviation of the P^c − S relationships from the equilibrium condition. Since the magnitude of becomes spectacularly large toward irreducible saturation, one can infer that the system properties at this region exhibit wider deviations from equilibrium.

Predictions of the relationships by ANN [9-15-17-1], ANN [9-10-10-1] and ANN [9-11-13-1] are shown in Figures 8–10, respectively. While these plots show good prediction of the relationships for the entire data set, they show more mismatches in comparison with Figure 7, especially at low water saturation or at high values of ⁠. However, ANN [9-11-13-1] and ANN [9-10-10-1] appear to predict better than ANN [9-15-17-1], having a lower number of mismatches especially at low water saturation.

Plots of the predictions of relationships by LR and NLR models are shown in Figures 11 and 12. The predictions by these models are less reliable, especially at low water saturation where the values are higher. The performances of these regression models are much less satisfactory in comparison to any of the ANN configurations with plots shown in Figures 7,⁸⁹–10. This reveals the limitations of the regression models in the prediction of relationships for two-phase flow systems.

Judging from the results of the statistical analyses on the prediction of as well as the above models’ performances on the prediction of relationships, one can conclude that ANN [9-13-15-1] is the best structure among the models tested in this work. This conclusion is similar to that of Hanspal et al. (2013). They found that the regression models performed poorly in the prediction of the relationships and concluded that generally the regression models have much less predictive ability than ANN structures for two-phase flow system characteristics. In their work, the display of NLR models seems better than shown in Figure 12 of this work even though similar functions were used. This can be explained by the fact that they utilize only five independent variables in their work as different from the nine used in this work. Therefore, one can infer that the performance of the regression models becomes less reliable as the number of independent variables increases for two-phase flow systems.

In the previous analyses and discussions, the ANN [9-13-15-1] structure is shown to be the best-performing network in the context of this work. In this section, the network is used to predict the domain scale dependency of relationships on the basis that the network is trained and validated. Separate data from different experiments are independently predicted. These include data from the literature as well as our in-house laboratory experiments.

Figures 13,¹⁴–15 display the results of the predictions. In the figures, for the original domain size is represented by two plots. One plot was obtained from the experimental data and the other was obtained from the ANN prediction of the original experimental data at the original domain scale (volume). These are labeled original volume (experiment) and the original volume (ANN) for the actual experimental data and the ANN prediction, respectively. The other two plots in each of the figures are obtained when the domain scale (volume) is increased by 10 and 20%, respectively, and the corresponding relationships are predicted using ANN [9-13-15-1]. The results show that increasing the volume of the domain increases the magnitude of ⁠. This is shown in Figure 13. At 10% increase in original domain volume, curve lies higher than at the original domain size. This effect becomes greater at 20% increase in domain volume.

This trend points to an important observation in the literature about the dependency on the domain scale. Several authors have reported the same phenomenon (Dahle et al. 2005; Bottero et al. 2011a, b). Bottero et al. (2011a, b) found that there is a shift to higher values in the relationships as the scale goes from local measurements to higher averaging windows. Also, Dahle et al. (2005) using a bundle of tubes model reports this phenomenon with a greater effect of scale on relationships which is said to be proportional to the square of the length of the domain.

Furthermore, using the experimental results by Goel & O'Carroll (2011), increasing the domain volume shows an increase in the relationships. This is shown in Figure 14. The relationship for 10% increases in the domain size lies above the curve from the original experimental domain size. A 20% increase in the domain size also shows further rise in the relationships.

Finally, results from our laboratory for a higher viscosity ratio silicone oil-water system (500) were also tested using the ANN [9-13-15-1]. This is shown in Figure 15. It can be seen from the figure that the curve rises as the domain size increases. But the rise in this case is rather sluggish, especially at higher water saturation. This can be attributed to the high viscosity ratio (500) in this case. The original curve of Das & Mirzaei (2013) shown in Figure 13 uses a viscosity ratio of 200 and the work of Goel & O'Carroll (2011) uses a viscosity ratio of 5. Viscosity ratio refers to the ratio of the viscosity of the oil to that of the water.

As demonstrated above, the trend in as the domain size increases, which indicates the dependency of the dynamic effects in the system properties of two-phase flow systems on the media characteristics. Since increase in indicates increased deviation from equilibrium (Das et al. 2007), the domain size certainly impacts the dynamic P^c– S relationships in two-phase flow in the porous media. Judging from Equation (1), increasing the magnitude of implies two things. Provided the value of P^c,static remains unaffected by domain scale, then the increase in the value of is influenced by increasing value of P^c,dyn and/or decreasing value of the ⁠. Thus, it may be rightly considered that P^c,dyn increases as the domain scale increases. Similarly, it is plausible to consider that decreases as the domain scale increases which may be caused by the decreasing pressure gradient as the domain height or size increases. While the simultaneous impact of the changes, mentioned above, can cause the observed effects on ⁠, the ratio of their contribution may differ. Bottero et al. (2011a, 2011b) observe that the marginal change in pressure difference (P^c,dyn − P^c,static) with upscaled windows of observation is less significant. They attribute the change of with domain scale to the ⁠, which decreases significantly as the domain size or length increases. Thus, plays significant role in the domain scale dependency of ⁠.

Application of an ANN for the prediction of the scale dependence of the dynamic capillary pressure effects in two-phase flow in porous media has been elaborately demonstrated.

Statistical analyses of the models tested showed that ANN configurations with double hidden layers outperformed those with single layers. Further comparison of the ANNs with LR and NLR models showed that ANNs have better prediction ability of the two-phase flow system characteristics.

Using the best-performing ANN structure (ANN [9-13-15-1]) in this work, the prediction of the domain size dependency for relationships reveals that the curve rises as the domain size increases in all the viscosity ratios tested. It was pointed out that the rate of change of saturation plays a more significant role in the domain scale dependency of ⁠.

Our findings showed the reliability and applicability of the ANN in characterizing and predicting the complex relationships for two-phase flow in porous media. Since the ANN system can be readily accessed and conveniently set up, it offers savings in cost and computational time in comparison to the flow-physics-based simulators.

The work in this paper was conducted in the framework of the Engineering and Physical Sciences Research Council (EPSRC), UK, project GR/S94315/01, ‘micro-heterogeneity and temperature effects on dynamic capillary pressure–saturation relationships for two-phase flow in porous media’. A PhD studentship awarded to L. K. A. by the Petroleum Technology Development Fund (PTDF), Nigeria, is gratefully acknowledged. Comments of three anonymous referees which have helped to improve the paper are appreciated.