Predicting the Compressibility Factor of Natural Gas by Using Statistical Modeling and Neural Network

: Simulating the phase behavior of a reservoir ﬂuid requires the determination of many parameters, such as gas–oil ratio and formation volume factor. The determination of such parameters requires knowledge of the critical properties and compressibility factor (Z factor). There are many techniques to determine the compressibility factor, such as experimental pressure, volume, and temperature (PVT) tests, empirical correlations, and artiﬁcial intelligence approaches. In this work, two different models based on statistical regression and multi-layer-feedforward neural network (MLFN) were developed to predict the Z factor of natural gas by utilizing the experimental data of 1079 samples with a wide range of pseudo-reduced pressure (0.12–25.8) and pseudo reduced temperature (1.3–2.4). The statistical regression model was proposed and trained in R using the “rjags” package and Markov chain Monte Carlo simulation, while the multi-layer-feedforward neural network model was postulated and trained using the “neural net” package. The neural network consists of one input layer with two anodes, three hidden layers, and one output layer. The input parameters are the ratio of pseudo-reduced pressure and the pseudo-reduced temperature of the natural hydrocarbon gas, while the output is the Z factor. The proposed statistical and MLFN models showed a positive correlation between the actual and predicted values of the Z factor, with a correlation coefﬁcient of 0.967 and 0.979, respectively. The results from the present study show that the MLFN can lead to accurate and reliable prediction of the natural gas compressibility factor.


Introduction
The relationship between pressure, volume, and temperature can be summarized by the equation of state (EOS).EOS has been widely used in the petroleum industry, especially when dealing with natural gas.Natural gas reservoirs play a vital role in responding to the massive global energy market.Knowledge of the physical properties of gas or oil related to pressure, volume, and temperature (PVT) is of great importance to petroleum engineers and academics.The physical properties of natural gas, which are used in the estimation of hydrocarbon reserves in reservoirs, the study of gas flow in the reservoir and wellbore and its thermodynamic behavior are essential in planning, gas metering, the design of pipelines, and surface facilities [1].Most of these properties can be determined in a PVT laboratory.One of the most important properties is the compressibility factor of a gas, because most of the other hydrocarbon gas properties depend directly or indirectly on it.The compressibility factor is also called the gas deviation factor and is represented by the Z factor.The deviation in behavior between the real gas and the ideal gas can be represented by the Z factor from the thermodynamic perspective.
The molecules of ideal gas are assumed not to be affected by the concept of corresponding states to have the same Z factor in the same two dimensionless properties: reduced pressure (P r ) and reduced temperature (Tr) [2,3].Consequently, the z-factor of any pure gas at a given reduced pressure (P r ) and reduced temperature (Tr) could be determined from the P r and Tr using the equation of state.The high expense of laboratory equipment and the lack of these laboratories are the main reason for seeking alternative techniques to determine or predict the Z factor of natural gases.The most well-known ways to predict the Z factor comprise the EOS, empirical correlations, and artificial intelligence [1,[4][5][6].In the past, EOS was developed and adapted to determine many of the physical properties of oil and gas, especially when dealing with volumetric properties and Z factor calculations [6][7][8].However, although the EOS gives good results for these properties, it has difficulties during the application, such as dealing with many parameters; and its specific mixing rules cause complexity during the process [9][10][11][12].Moreover, it features difficulties in estimating the critical properties of a large number of components, as in oil and gas mixtures.Therefore, many researchers began to conduct empirical correlations based on specific data sets to predict oil and gas physical properties, such as compressibility factor values [13][14][15][16][17].
Neural networks are excellent at simulating equations and correlations.Simple models have the ability to fit any practical function successfully.The combination of experimental PVT data, the EOS, and neural network approaches could lead to the proposal of a precise neural network model for predicting the behavior of the compressibility factor (Z factor) of natural gases.The aim of this work is to calculate the compressibility factor of natural hydrocarbon gases by two different techniques: statistical regression model and multilayered feedforward network.It was noted that there was a nonlinear relationship between the compressibility factor and both the pseudo-reduced pressures and temperatures.Two models for the Z factor, based on statistical analysis and multilayered feedforward network (MLFN), were employed to conduct a precise model for predicting the Z-factor of gas reservoirs.Statistical regression analysis was utilized to test the data sets and recognize the hidden patterns in the data.In addition, the statistical regression model could be used as a pre-modeling step for the MLFN.It is important to mention that the MLFN developed in this work does not provide a model for natural hydrocarbon gases; rather it is simply a reliable, precise, accurate, and effective method to calculate the compressibility factor of natural gases.
A literature review of different methods for calculating the compressibility factor is offered in the next section.

Literature Review
Empirical correlations simply and easily predict physical properties, so they are widely utilized in the oil and gas field.Based on the experimental PVT data for gas reservoirs (non-ideal gases) and using the pseudo-reduced pressure and pseudo-reduced temperature of these gases, Standing and Katz presented their standard chart for the petroleum industry to estimate the Z-factor for natural gas.The following two equations (Equations ( 1) and ( 2)) are used to estimate the P r and Tr.P r = P /P c (1) where P c and T c are the critical pressure and critical temperature, respectively.Later, several [13] attempts were conducted to correlate the Z-factor mathematically and to fit with the Standing and Katz chart [16].Some correlations were conducted by several researchers, such as Biggs and Brill in 1973 [15], Hall and Yarbrough in 1974 [18], Dranchuk and Abou-Kassem in 1975 [17], Kumar in 2004 [19], Azizi et al. in 2010 [20], among others [13,14,21].
Hall-Yarbrough [18] presented their correlation for the first time in 1973 in the following form: Energies 2022, 15, 1807 3 of 15 where P pr is the pseudo-reduced pressure, t is the reciprocal of the reduced temperature, and Y can be determined by solving Equation (4).
Despite the high confidence level of this correlation, it cannot be utilized when the pseudo-reduced temperature is less than one.Biggs and Brill presented their correlation to evaluate the Z factor of natural gases via a correlation based on the Standing and Katz compressibility chart [15].The correlation was limited to a pseudo-reduced temperatures higher than 0.92.The Dranchuk and Abou-Kassem (DAK) correlation for the calculations of gas density and compressibility is currently considered a standard in the petroleum industry [1].The correlation was conducted based on Benedict's equation of state, as follows [17].
The correlation is limited to the range of P r (0.2-30) and T r (1)(2)(3).In addition, to proceed with the DAK correlation, assuming an initial value of Z is necessary.At the beginning of 2014, Kumar et al. introduced the Shell Oil Company correlation to estimate the Z factor of gases as follows: The This correlation covers a wider range of both pressures and temperatures than that of Biggs and Brill.Despite the wide spread of empirical correlations in the petroleum industry, they have some limitations, such as inadequate prediction at elevated parameters, including pressure and temperature.In addition, some correlations are convenient for specific compositions and not widely used.For all these reasons, the researchers looked for a new way to overcome the disadvantages of empirical correlations.Currently, artifi-cial intelligence systems are used in various fields, especially in the petroleum sector, to overcome the drawbacks of empirical correlations, with promising levels of success [16,22].An artificial neural network (ANN) model containing only one hidden layer was used by Normandin et al. as an early attempt to predict the compressibility factor of some pure gases [23,24].The model was used only for pure gases, with inaccurate results; moreover, it cannot be used for natural gas mixtures.Mohanty et al. predicted the vapor phase and the bubble pressure of carbon dioxide-difluoromethane using ANN with a hidden layer [25].Saemi et al. reported that the linked adaptive neural network (LANN) and genetic algorithm (GA) could be used successfully to estimate reservoir permeability [26].
Currently, most of the physical properties of hydrocarbon gas and oil can be predicted using ANN [27], and there was another attempt to use ANN to predict the Z factor for hydrocarbon gas mixtures by Kamyab et al., using data from the Standing and Katz chart [16].The ANN model was applied using two input parameters, pseudo-reduced pressure and pseudo-reduced temperature, with two hidden layers.The proposed model was more accurate than the widely used and most popular correlation conducted by Dranchuk and Abo-Kassem.Fayazi et al. used the least square support vector machine (LSSVM) model to predict the Z factor of natural gas [13].Saghafi et al. introduced a model for the Z factor calculation of gas condensate based on the ANN and the genetic programming framework [28].

Methodology
The study aimed at predicting the Z factor from the pseudo-reduced pressure (P pr ) and temperature (T pr ).The first step was data gathering, which aimed at gathering the Z factor values and calculating the (P pr ) and (T pr ) values that were used to predict the Z factor.Secondly, data were explored and visualized in order to determine the relationship between the model variables.Based on the first two steps, two model types were selected and postulated.The first model was statistical regression model and the second was MLFN.Finally, the models were assessed and compared based on the correlation coefficient and residual analysis.The graphical representation of the workflow is shown in Scheme 1.
Energies 2022, 15, x FOR PEER REVIEW 5 of 17 Scheme 1.The graphical representation of the statistical model.

Data Gathering
The measured compressibility factor, compositional analysis, and molecular weight, as well as pressures and temperatures of a wide range of natural lean, rich, sweet, and sour gases, were gathered from the literature in order to predict the gas compressibility factor by using statistical modeling and MLFN [29][30][31][32][33][34].The pseudo-reduced pressure (Ppr) and temperature (Tpr) of each sample were calculated according to Kay's rule, using the following equations:

Data Gathering
The measured compressibility factor, compositional analysis, and molecular weight, as well as pressures and temperatures of a wide range of natural lean, rich, sweet, and sour gases, were gathered from the literature in order to predict the gas compressibility factor by using statistical modeling and MLFN [29][30][31][32][33][34].The pseudo-reduced pressure (P pr ) and temperature (T pr ) of each sample were calculated according to Kay's rule, using the following equations: P pr = P P pc and T pr = T T pc P pc = ∑ n i=1 y i P ci and T pc = ∑ n i=1 y i T ci where: P pc is the pseudo critical pressure of the natural gas mixture; T pc is the pseudo critical temperature of the natural gas mixture; P ci is the critical pressure of component i in the natural gas mixture; T ci is the critical temperature of component i in the natural gas mixture; y i is the mole fraction of component i in the natural gas mixture.The collected data from the literature are provided in the supplementary material (Tables S1 and S2).

Data Exploration
The gathered data were explored by calculating the statistical properties, such as the mean, standard deviation, and range.This step is important to assess the data quality and check whether there are any outliers and how to handle them.The next step is to visualize the relationship between the Z-factor and each of the (P pr ) and (T pr ), which is essential to decide whether the exploratory variables can be used to directly predict the Z-factor or whether they need to be transformed into new variables to facilitate the model selection process.
Based on data exploration and visualization, the Z-factor was found to be difficult to predict directly from the (P pr ) and (T pr ).Accordingly, variable transformation is needed to recognize the hidden patterns in the relationship between the Z-factor and each of the (P pr ) and (T pr ).One of the common ways to transform variables is to combine different variables together.Accordingly, the ratio (X) was obtained by dividing the P pr by the T pr .
Another pre-modeling step is to normalize the variables by subtracting each variable's observation from the mean and dividing the result by the standard deviation.This process is called "Centering and Scaling".It can be performed by using the "scale" function of the "base" package in R [35].This step is vital to enhance the stability of the model to be postulated [36][37][38].

Model Selection
The relationship between the Z-factor and (X) was found to be curvilinear.Therefore, a quadratic regression model was thought to be more appropriate than a linear model.In order to prove this, the "lm" function in R was used to compare the linear and quadratic cases.The comparison criteria between the linear and quadratic models, were the R-squared value and standard error.Moreover, the residuals were compared quantitatively and plotted in order to determine which model would be more representative of the data.Finally, it was proven that a quadratic model should be selected.
As the linearity between the Z-factor and X is not particularly high, another model is needed to better estimate the Z-factor based on machine learning.Recently, neural networks have been widely used to solve non-linear problems due to their advantages over regression models [39,40].Neural networks are among the supervised-learning tools that provide a function based on a network of input(s) and output(s) [41].Therefore, a neural network is expected to be more accurate than the regression model.
In this study, the multi-layer-feedforward neural network (MLFN) was selected as it is widely used for general-purpose non-linear regression models and due to its ability to extract hidden patterns from data.Another advantage of the MLFN is that it can perform input-output mapping by modifying the weights of the coefficients until reaching the least difference between the desired output and the network's actual output [42].

Model Postulation Statistical Regression Modeling
Based on the Z-X relationship, the model is expected to follow the quadratic formula: where Z is the output (response) variable, X is the P r -over-T r ratio (explanatory variable), and the coefficients of the equations are β 1 , β 2 , and β 3 .Given the observations of the Z-factor and X, the models' coefficients can be obtained, based on the Bayes' theorem [43], from the joint distribution of the likelihood and priors, as shown below: P(β |y ) = P( y| β)P(β) P(y) where P(β|y) is the posterior probability of the coefficients' vector (β) given the likelihood (y), which represents the actual observations, the expression "P(y|β)P(β)" is the joint distribution of the likelihood and the priors of the coefficients, while P(y) is the marginal distribution of the likelihood (y).Therefore, the posterior probability is equivalent to the joint distribution of the likelihood and priors, as shown below: In some cases, the joint probability of the model is too complex to be integrated.That's why the Markov chain Monte Carlo (MCMC) simulation can be used to obtain the posterior distribution of such models [44].A Markov chain (MC) is a chain of numbers in which each number depends on the previous number in the sequence.If {x 1 , x 2 , x 3 , . . ., x t } is a Markov chain, where {1, 2, 3, . . ., t} are successive points in time, the probability of the variables can be expressed according to the chain role, as shown below: P(x 1 , x 2 , x 3 , . . . ,x t ) = P(x 1 ) P(x 2 |x 1 ) P(x 3 |x 2 , x 1 ), . . ., P(x t |x t−1 , x t−2 , . . . ,x 2 , x 1 ) Assuming that the transition probabilities are not time-dependent, the probability at the time (t) depends only on the value of {x t−1 }: P(x 1 , x 2 , x 3 , . . .x n ) = P(x 1 ) P(x 2 |x 1 ) P(x 3 |x 2 ), . . ., P(x t |x t−1 ) The stationary distribution of the MC is the initial distribution at which the transition probability does not change in any given state.The aim of the MCMC is to draw multiple random values from a proposed distribution so that the sequence of all simulations is a Markov chain, and the stationary distribution of that chain is the posterior distribution.According to the law of large numbers [45], the MC should converge to the true mean of the posterior distribution that can be assigned to the coefficient.Therefore, each coefficient in the model is obtained by calculating the posterior mean of all the realizations drawn from the proposed distribution.
The aim of the statistical regression model is to fit a linear relationship between the Z-factor, which is the response variable, and each of the squared-X and X, which are the explanatory variables of the model.Accordingly, the likelihood function of the model can be expressed as shown below: Given the coefficients' vector, the explanatory variables X 2 , and the variance σ 2 , the Z-factor follows a normal distribution.The mean of this distribution equals the linear expression of the covariates and coefficients.The symbol (i) is the observations' index.
The next layer of the model consists of the coefficients' vector β and the variance σ 2 , which are considered the priors of the model.In order to determine an idea as to the probability distribution function of each of the coefficients, the "lm" function of the "stats" package [39] was used, in R, to fit a preliminary linear model that resulted in the hyper-parameters (mean and standard deviation) of each coefficient.Accordingly, the prior functions of the three coefficients can be expressed as follows: where the coefficient β j is normally distributed with a mean µ j and standard deviation τ j .The variance of the model follows an inverse-gamma distribution, which is a type of continuous probability distribution [46] that is always positive and depends on two parameters: the shape parameter α and the scale parameter β.The inverse-gamma is used to draw the unknown variance of a normal distribution, where the priors are uninformative.Therefore, the variance (σ 2 ) is given by: where α and β are the priors of the variance σ 2 .The hyper-parameters of the variance's distribution were treated as non-informative priors.The values of the α and β were assumed to be greater than 1 to include high random-error values.
The aim of the model is to use the MCMC simulation to randomly draw multiple realizations of the coefficients from their common distribution and, then, to calculate the Z-factor by substituting the posterior means of the coefficients to Equation (8).The graphical representation of the statistical model is shown in Figure 1, where: Z i is the response variable of an observation (i); X 2 i and X i are the explanatory factors; i is the observations' index, which ranges from 1 to n; b j is the coefficients vector, which depends on the mean µ and standard deviation τ; and σ 2 is the variance of the model, which depends on the shape parameter α and scale parameter The statistical model was postulated and trained in R using the "rjags" package, which uses the MCMC to generate a sequence of dependent samples from the posterior distribution of the parameters [47,48].

MLFN Model
Recently, neural networks have been widely used to solve non-linear problems due to their advantages over regression models [39,40].One of the common neural networks is the multi-layer feedforward (MLFN) neural network, which consists of three main types of layers: an input layer, one or more hidden layers, and an output layer [49].The advantage of this network is that hidden layers enable it to estimate variables based on nonlinear relationships.
The MLFN layers are fully connected, so that all units in one layer are connected to all units in the next layer [49].Each connection is controlled by a random weight that should be adjusted by an algorithm, such as the backpropagation method [50], until the optimum set of weights is obtained, so that the predicted output best matches the actual output.The aim of this algorithm is to feed the hidden layer(s) with the input units associated with their weights and, then, to produce a function to calculate the output from the weighted sums of the inputs.
The structure of the MLFN in the current study consists of three layers: an input layer, three hidden layers, and an output layer.The input layer consists of two neurons, which are X and X 2 , while the output layer consists of one neuron, which is the Z-factor.The number of neurons in the hidden layer is adjusted until the best performance is achieved.The network was postulated and trained using the "neural net" package in R, which calculates the generalized weights [51] and provides customized settings for error and activation function.

MLFN Model
Recently, neural networks have been widely used to solve non-linear problems due to their advantages over regression models [39,40].One of the common neural networks is the multi-layer feedforward (MLFN) neural network, which consists of three main types of layers: an input layer, one or more hidden layers, and an output layer [49].The advantage of this network is that hidden layers enable it to estimate variables based on non-linear relationships.
The MLFN layers are fully connected, so that all units in one layer are connected to all units in the next layer [49].Each connection is controlled by a random weight that should be adjusted by an algorithm, such as the backpropagation method [50], until the optimum set of weights is obtained, so that the predicted output best matches the actual output.The aim of this algorithm is to feed the hidden layer(s) with the input units associated with their weights and, then, to produce a function to calculate the output from the weighted sums of the inputs.
The structure of the MLFN in the current study consists of three layers: an input layer, three hidden layers, and an output layer.The input layer consists of two neurons, which are X and X 2 , while the output layer consists of one neuron, which is the Z-factor.The number of neurons in the hidden layer is adjusted until the best performance is achieved.The network was postulated and trained using the "neural net" package in R, which calculates the generalized weights [51] and provides customized settings for error and activation function.
Eventually, the statistical and MLFN models were compared based on the correlation coefficient and residual analysis.The correlation coefficient was obtained using the "cor.test"function of the "stats" package in R based on Pearson's product moment correlation coefficient [52].The residual analysis was performed by plotting the residuals to show how far the residuals were from the zero line [53].

Data Gathering
The collected data from different sources in the literature included about 1079 data points for each of the Z-factor, P pr , and T pr to achieve the goal of this study.The statistical distribution of the data points, such as the minimum (Min), maximum (Max), mean and standard deviation (Std), are reported in Table 1.

Data Exploration
It is clear that the minimum and maximum Z values of the data points were far from the mean.In addition, the Z-P pr and Z-T pr scatter plots (Figure 2a,b, respectively) show some residuals.Accordingly, the residuals were removed to avoid high variance in the data.
Energies 2022, 15, x FOR PEER REVIEW 9 of

Data Gathering
The collected data from different sources in the literature included about 1079 da points for each of the Z-factor, Ppr, and Tpr to achieve the goal of this study.The statistic distribution of the data points, such as the minimum (Min), maximum (Max), mean an standard deviation (Std), are reported in Table 1.

Data Exploration
It is clear that the minimum and maximum Z values of the data points were far fro the mean.In addition, the Z-Ppr and Z-Tpr scatter plots (Figure 2a,b, respectively) sho some residuals.Accordingly, the residuals were removed to avoid high variance in t data.Another observation is that the Ppr values were far from those of the Z-factor an Tpr, as shown in the box plots in Figure 3.This is why each variable was normalized b centering and scaling.
The relationship between the Z-factor and each of the Ppr and Tpr seemed non-line and showed a high degree of randomness.In addition, there were some residuals th were far away from the data range.Therefore, the data should be cleaned and the vari bles transformed in order to be ready for modeling.
It is obvious how the Z-X relationship (Figure 4) was more meaningful and less noi than the Z-Ppr and Z-Tpr relationships (Figure 2a,b, respectively).The next step was postulate a model that could estimate the Z-factor from the X.Two methods were appli to do so: statistical regression and multi-layer feedforward neural network (MLFN).Another observation is that the P pr values were far from those of the Z-factor and T pr , as shown in the box plots in Figure 3.This is why each variable was normalized by centering and scaling.
The relationship between the Z-factor and each of the P pr and T pr seemed non-linear and showed a high degree of randomness.In addition, there were some residuals that were far away from the data range.Therefore, the data should be cleaned and the variables transformed in order to be ready for modeling.

Results of the Statistical Regression Model
The data points were divided into training and testing points with the percentages 75% and 25%, respectively.Figure 5 shows the trace-plots of the Markov chains, which aim to estimate the models' coefficients and precision (prec), which is the reciprocal of the variance.The number of iterations is shown on the X-axis and the value of each coefficient is shown on the Y-axis.It is obvious how the Z-X relationship (Figure 4) was more meaningful and less noisy than the Z-P pr and Z-T pr relationships (Figure 2a,b, respectively).The next step was to postulate a model that could estimate the Z-factor from the X.Two methods were applied to do so: statistical regression and multi-layer feedforward neural network (MLFN).

Results of the Statistical Regression Model
The data points were divided into training and testing points with the percentages 75% and 25%, respectively.Figure 5 shows the trace-plots of the Markov chains, which aim to estimate the models' coefficients and precision (prec), which is the reciprocal of the variance.The number of iterations is shown on the X-axis and the value of each coefficient is shown on the Y-axis.

Results of the Statistical Regression Model
The data points were divided into training and testing points with the percentages 75% and 25%, respectively.Figure 5 shows the trace-plots of the Markov chains, which aim to estimate the models' coefficients and precision (prec), which is the reciprocal of the variance.The number of iterations is shown on the X-axis and the value of each coefficient is shown on the Y-axis.Each parameter was simulated by using three chains, and was expected to converge at the mean value of its multiple realizations.The chains reached the stationary distribution at a number of iterations that ranged from 1500 to 2000.The mean values of the three coefficients (β1, β2, and β3) were 0.6, 0.45, and 0.08, respectively.When applying these values to Equation (8), the following equation is obtained: Z = 0.6 X 2 + 0.45 X + 0.08 After un-scaling the predicted and actual values of the Z-factor for both the training and testing data, the correlation coefficient was calculated using the "cor.test"function of the "stats" package in R. Figure 6a,b show the predicted-versus-actual plots for the train- Each parameter was simulated by using three chains, and was expected to converge at the mean value of its multiple realizations.The chains reached the stationary distribution at a number of iterations that ranged from 1500 to 2000.The mean values of the three coefficients (β 1 , β 2 , and β 3 ) were 0.6, 0.45, and 0.08, respectively.When applying these values to Equation (8), the following equation is obtained: Z = 0.6 X 2 + 0.45 X + 0.08 After un-scaling the predicted and actual values of the Z-factor for both the training and testing data, the correlation coefficient was calculated using the "cor.test"function of the "stats" package in R. Figure 6a,b show the predicted-versus-actual plots for the training and testing data, resulting in the correlation-coefficient values 0.969 and 0.967, respectively.
The plots show a positive correlation between the predicted and actual values.However, there are some residuals in the bottom-left of the plot.
ing and testing data, resulting in the correlation-coefficient values 0.969 and 0.967, respectively.The plots show a positive correlation between the predicted and actual values.However, there are some residuals in the bottom-left of the plot.

Results of the MLFN Model
The graphical representation of the MLFN model is shown in Figure 7, where the network consists of an input layer with two nodes (green circles), three hidden layers (black circles), and an output layer (blue circle).The number of hidden layers and the number of nodes in each of them was adjusted until reaching the least mean-square error (MSE).The best network performance was reached after 2011 steps by using 2, 5, and 2 nodes in the hidden layers, resulting in an MSE value equal to 0.461.The network was applied to the data, resulting in a strong positive correlation between the predicted and actual Z-factor, as shown in Figure 8, for the training and testing data, respectively.The correlation coefficient of the training and testing data were 0.982 and 0.979, respectively.

Results of the MLFN Model
The graphical representation of the MLFN model is shown in Figure 7, where the network consists of an input layer with two nodes (green circles), three hidden layers (black circles), and an output layer (blue circle).The number of hidden layers and the number of nodes in each of them was adjusted until reaching the least mean-square error (MSE).The best network performance was reached after 2011 steps by using 2, 5, and 2 nodes in the hidden layers, resulting in an MSE value equal to 0.461.The network was applied to the data, resulting in a strong positive correlation between the predicted and actual Z-factor, as shown in Figure 8, for the training and testing data, respectively.The correlation coefficient of the training and testing data were 0.982 and 0.979, respectively.
ing and testing data, resulting in the correlation-coefficient values 0.969 and 0.967, respectively.The plots show a positive correlation between the predicted and actual values.However, there are some residuals in the bottom-left of the plot.

Results of the MLFN Model
The graphical representation of the MLFN model is shown in Figure 7, where the network consists of an input layer with two nodes (green circles), three hidden layers (black circles), and an output layer (blue circle).The number of hidden layers and the number of nodes in each of them was adjusted until reaching the least mean-square error (MSE).The best network performance was reached after 2011 steps by using 2, 5, and 2 nodes in the hidden layers, resulting in an MSE value equal to 0.461.The network was applied to the data, resulting in a strong positive correlation between the predicted and actual Z-factor, as shown in Figure 8, for the training and testing data, respectively.The correlation coefficient of the training and testing data were 0.982 and 0.979, respectively.

Discussion
The purpose of engaging both the statistical and MLFN models in this study was to use the regression model as a pre-modeling step for the neural network.In other words, it is important to use statistical models to recognize the hidden patterns in data and test various sets of model variables in order to obtain the best group of predictors that could be later used as inputs to neural networks.Another benefit of combining statistical models and neural networks is that it shows how machine-learning could overcome the randomness inside data and recognize the patterns hidden in them.This is clear when comparing the actual-predicted Z-factor plots of the statistical and MLFN models, as shown in Figures 5 and 8, respectively, where the MLFN model shows stronger linearity and less residuals relative to the statistical model.This hypothesis is evidenced by the correlationcoefficient values, which were higher in MLFN compared to the regression model, as shown in Table 2.Both models showed a high degree of stability because the correlation-coefficient values of the training and testing data were so close to each other, which indicates that there was no over-fitting.However, a more precise comparison was drawn between the two models by using the residual analysis shown in Figure 9a,b, where the residuals of the statistical and MLFN models are plotted, respectively.Unlike the statistical model, the MLFN model showed residuals that were closer to the zero red line and more centered around it.This is how the MLFN model improved the Z-factor prediction by using the power of machine learning.

Discussion
The purpose of engaging both the statistical and MLFN models in this study was to use the regression model as a pre-modeling step for the neural network.In other words, it is important to use statistical models to recognize the hidden patterns in data and test various sets of model variables in order to obtain the best group of predictors that could be later used as inputs to neural networks.Another benefit of combining statistical models and neural networks is that it shows how machine-learning could overcome the randomness inside data and recognize the patterns hidden in them.This is clear when comparing the actual-predicted Z-factor plots of the statistical and MLFN models, as shown in Figures 5 and 8, respectively, where the MLFN model shows stronger linearity and less residuals relative to the statistical model.This hypothesis is evidenced by the correlation-coefficient values, which were higher in MLFN compared to the regression model, as shown in Table 2.Both models showed a high degree of stability because the correlation-coefficient values of the training and testing data were so close to each other, which indicates that there was no over-fitting.However, a more precise comparison was drawn between the two models by using the residual analysis shown in Figure 9a,b, where the residuals of the statistical and MLFN models are plotted, respectively.Unlike the statistical model, the MLFN model showed residuals that were closer to the zero red line and more centered around it.This is how the MLFN model improved the Z-factor prediction by using the power of machine learning.

Conclusions
The gas compressibility factors of natural gases were calculated using a multi-layerfeedforward neural network and statistical regression as a function of the ratio of the

Scheme 1 .
Scheme 1.The graphical representation of the statistical model.

Energies 2022 , 17 Figure 1 .
Figure 1.The graphical representation of the statistical model.

Figure 1 .
Figure 1.The graphical representation of the statistical model.

Figure 2 .
Figure 2. Cross plots of (a) Z vs. P pr and (b) Z vs. T pr .

Figure 3 .
Figure 3. Box plots of the range values of Ppr, Tpr, and Z factor.

Figure 4 .
Figure 4. Cross plots of Z vs. X.

Figure 3 .
Figure 3. Box plots of the range values of P pr , T pr , and Z factor.

Energies 2022 , 17 Figure 3 .
Figure 3. Box plots of the range values of Ppr, Tpr, and Z factor.

Figure 4 .
Figure 4. Cross plots of Z vs. X.

Figure 4 .
Figure 4. Cross plots of Z vs. X.

Figure 5 .
Figure 5.The trace-plots of the Markov chains that simulate the statistical model's coefficients (β1, β2, and β3) and the precision (prec).

Figure 6 .
Figure 6.Plots of predicted Z-factor vs. actual Z-factor for the statistical models of both the training and testing data.

Figure 6 .
Figure 6.Plots of predicted Z-factor vs. actual Z-factor for the statistical models of both the training and testing data.

Figure 6 .
Figure 6.Plots of predicted Z-factor vs. actual Z-factor for the statistical models of both the training and testing data.

Figure 7 .
Figure 7.The graphical representation of the MLFN model.

Figure 7 .
Figure 7.The graphical representation of the MLFN model.

Figure 8 .
Figure 8. Plots of predicted Z-factor vs. actual Z-factor for the MLFN model of both the training and testing data.

Figure 8 .
Figure 8. Plots of predicted Z-factor vs. actual Z-factor for the MLFN model of both the training and testing data.

Figure 9 .
Figure 9.The residual plots of the (a) statistical and (b) MLFN models.

Table 1 .
Data statistics before cleaning.

Table 1 .
Data statistics before cleaning.

Table 2 .
The correlation-coefficient values of the training and testing data for the statistical and MLFN models.

Table 2 .
The correlation-coefficient values of the training and testing data for the statistical and MLFN models.