Spatial-Temporal Separation Based on the Dynamic Recurrent Wavelet Neural Network Modelling for ASP Flooding

In this paper, a three-dimensional spatial-temporal decomposition modelling method is proposed to build the alkali-surfactant-polymer (ASP) flooding model, in which a new dynamic recurrent wavelet neural network (DRWNN) is presented to identify the temporal coefficients. At first, the detailed mathematical model of ASP flooding is described which is a complex distributed parameter system. Then a three-dimensional spatial-temporal modelling method is inferred based on Karhunen-Loeve (K-L) decomposition to decompose the water saturation of reservoir into a series of spatial basis functions and corresponding temporal coefficients. Furthermore, the recurrent wavelet neural network is used to acquire the identification model, in which the injection concentrations of ASP flooding and temporal coefficients are taken as the input and output information. In order to improve the capability of dynamic modelling, DRWNN is proposed through adding feedback layers and setting the different weights with time to achieve dynamic memory of the past information. Considering the gradient descent method for the neural networks training easily leads to local minimum and slow convergence speed, the spectral conjugate gradient method is introduced to optimize the weights of DRWNN. At last, DRWNN is used to build the relation between the moisture content of production wells and the water saturation of the corresponding grids. Thus, the final approximate model of ASP flooding is finished. The accuracy is proved by model with four injection wells and nine production wells through data from the mechanism model.


Introduction
With the old oil fields entering the later period of development, moisture content of reservoir is increasing, and the oil production is reducing [1]. How to update the technical means to ensure oil recovery is one of the most important measures to stabilize oil production. ASP flooding is an important tertiary oil recovery technique which is widely studied [2]. It can enhance the oil production evidently by use of the interaction among alkali, surfactant and polymer to improve the physicochemical property of reservoir. However, there is lack of a uniform mathematical model description because of the uncertainty of alkali reaction. It is important to build an accurate and easily to be applied model in the research of ASP flooding.
In addition, since the mechanism model of ASP flooding is a complexed distributed parameter system, it is hard to carry out optimal control strategies for the oil production, because of the features of infinite dimensions, spatial-temporal coupling and complex nonlinear behavior. The inputs of ASP flooding are the injection concentrations of alkali, surfactant and polymer; the output is the water cut of production wells. The state variables of ASP flooding contain water saturation, pressure and grid concentrations. In applications, different methods are used to search the injection concentrations of ASP flooding to get the best performance index. One of common ways is to solve the mechanism equations of ASP flooding, such as fluid equations and seepage equations, etc. [3] However, the method based on mechanism model usually involves a lot of math operations, and the process of mathematic treatment always needs a fair amount of calculation. Though many people have done researches on modelling of ASP flooding, nearly all the works are about improving or enriching the primary model. The main work they have done is to consider more influencing factors or inductive researches; the model is still complex and difficult to be applied. So it is important to consider using model identification method to approximate the ASP flooding system so as to simplify the mathematical operation process, and it will be helpful to improve the computational efficiency and reduce the computational complexity of the algorithm [4].
Traditional system identification methods do not consider spatial information [5]. To model the spatial-temporal dynamics, spatial information is a very important part of system. Spatial-temporal decomposition technology [6,7], which comes from Fourier series expansion, is a very useful modelling method for distributed parameter systems. It can reflect the information in time domain and space domain of the system by decomposition. A spatiotemporal variable can be expanded into an infinite number of spatial basis functions and corresponding temporal coefficients. Generally speaking, the first few primary spatial basis functions can reflect the maximal inner information of the system, which provides a good approximation because of their separation properties [8]. In this way, the spatial-temporal method can obtain a finite-dimensional model.
An important condition of modelling is to guarantee the accuracy of model which is highly dependent on the choice of spatial basis functions. In particular, Karhunen-Loeve (K-L) decomposition, which is called proper orthogonal decomposition and principal component analysis [9,10], is a popular spatial-temporal decomposition approach to find principal spatial structures and reduce the dimension of the data. Among all linear expansions, K-L expansion is the most efficient in the sense that, for a given approximation error, the number of K-L bases required is minimal. As a result, K-L decomposition can help to reduce the model dimension and the number of estimated parameters.
Once the spatial basis functions are designed properly, the corresponding states can be determined by projecting the spatiotemporal data onto these spatial basis functions. To model the input-state dynamics, the relationship between inputs and temporal coefficients should be identified. In the process of modelling, the dynamic recurrent neural network is a good way for the system with unknown nonlinearity. The recurrent wavelet neural network (RWNN) [11,22] combining wavelet transform with dynamic neural network can get better modelling capability. It is widely used in nonlinear dynamical modelling problems because of the advantages of neural networks such as reliable theory basis, explicit practical sense, simple algorithm realization and strong adaptability, etc.
In order to improve the dynamic modelling capability of RWNN, a new dynamic recurrent wavelet neural network (DRWNN) is proposed, in which the feedback layers are added for the network and different weights are set for feedback layers which is decreasing with time to make the network have the dynamic memory function. However, the gradient descent method is commonly adopted to train the parameters of network, which may cause the local minimum solutions and slow convergence speed. Especially in large search space for multimodal function, the function usually can't find out the global optimal value [12]. Many researchers have studied the optimization methods, which can be divided into five classes in general: gradient descent methods, Newton methods, conjugate gradient methods, heuristic optimization methods and Lagrange multiplier methods [13]. But the gradient descent method has slow converge speed, the Newton method needs to storage and compute the Hessen matrix which has huge computing complexity, the results of heuristic optimization method has a certain uncertainty, and the Lagrange multiplier method is often used to deal with multi-objective optimization problems. While the conjugate gradient method, which balances the convergence speed and computing complexity, has good stability and with no need of any extraneous parameter. It is widely used in nonlinear optimization. PRP method is viewed as one of the most effective conjugate gradient (CG) methods, but the convergence property is not so good [14]. Reference [15] gave an improved PRP method which can also be named as NPRP method. By modifying the step-length parameter, the sufficient descent property and global convergence are guaranteed on the condition of strong wolf line search. To improve the performance of numerical methods further, on the basis of [15] a new spectral NPRP (SNPRP) conjugate gradient method was presented through adding the spectral parameter and revising the step-length [16]. The author also gave the proof about global convergence and fast descent. On account of these advantages, the SNPRP conjugate gradient method is adopted to train the weights of DRWNN.
In this paper, a new model of ASP flooding is obtained in terms of spatial-temporal decomposition. K-L decomposition is used to model the state parameters of reservoir (water saturation) into a series of temporal coefficients and spatial basis functions. To get the relation between the injection concentrations of ASP flooding and temporal coefficients, the DRWNN is proposed in which different weights decreasing with time are added for feedback layers to improve the dynamic modelling capability and the weights are trained by the SNPRP conjugate gradient method. Furthermore, build the model between the grid water saturation and the moisture content of production wells with DRWNN. Then the final model of ASP flooding is finished with the injection concentrations as the inputs and the moisture content as the outputs.
This article is organized as follows: The mechanism model description of ASP flooding is described in section 2. In section 3, the Karhunen-Loeve decomposition of three-dimensional state of ASP flooding is introduced. In section 4, a new dynamic recurrent wavelet neural network, named DRWNN, is proposed to identify the relationships in the modelling. Besides that, in order to improve the ability of DRWNN, the SNPRP conjugate gradient method is adopted to train the weights. Section 5 gives an illustrative simulation example. Finally, a few conclusions are presented in section 6.

Mechanism Model Description of ASP Flooding
Suppose that the region of the reservoir is ( , , ) x y z ∈ Ω , 3 R Ω ∈ , the main five components of ASP flooding are oil, water, polymer, surfactant, alkali. The mathematical model of ASP flooding can be described as: [19,20] The flow equation for oil phase: The flow equation for water phase: The flow equation for polymer phase: The flow equation for alkali phase: The flow equation for surfactant phase: The initial conditions: The boundary conditions: Polymer can improve the viscosity of pore water, so it can be represented as: where 0 w µ is the viscosity of pure water; 1 ap , 2 ap , 3 ap , sp are empirical constants; sep C is the salinity of reservoir.
Alkali can react with substances in the reservoir to produce new substances and cause alkali consumption. This process also influences oil-water interfacial tension, and the concrete mathematical description is where ( ) wo s C σ is oil-water interfacial tension when oil displacement agent is only surfactant.
The different oil displacement agents can affect the strength of the rock adsorption and can be described as: where ll a , b are adsorption parameters.
Because the alkali will lower adsorption of polymer and surfactant on the surface of the pores of rock, the adsorption should be recalculated when alkali is injected to reservoir, and the formula for calculation can be described as: The other parameters and model explanations can be found in the appendix.

The Three-Dimensional Spatial-Temporal Decomposition Modelling Method for ASP Flooding
ASP flooding system is a complex distributed parameter system which is very difficult to be modeled in general method. In this paper, a new three-dimensional spatial-temporal model method is proposed based on Karhunen-Loeve decomposition method [8]. During the modelling process for ASP flooding, the main problem is to identify an appropriate spatial-temporal model according to the input, output and system state data. The whole model for space/time separation of ASP flooding can be divided into three stages: (1) Decompose the states of ASP flooding into spatial basis functions and temporal coefficients; (2) Identify the relationship between the injection concentrations and the temporal coefficients; (3) Model the moisture content of production wells and corresponding grid states.

Three-Dimensional Spatial-Temporal Method Based on KL Decomposition
For the system of ASP flooding, suppose that the output of system is According to the theory of Fourier series, the spatial-temporal output ( ) , , , Y x y z t can be decomposed into a series of orthonormal spatial basis functions and temporal coefficients. So the decomposition form can be expressed as: In engineering applications, ensuring the accuracy of the approximation, the model output is usually truncated to the finite form: where n denotes the order of model, it mainly depends on the system spatial basis functions which reflect the dynamic characteristic of system. As the description above, it is clear that the main problem of modelling is to get the dominant spatial basis functions The problem can be converted into a minimization problem as the following objective function: Because the spatial basis functions are unit orthogonal, they are submitted to the below equation: (15) Thus, the temporal coefficients can be described as: Considering the orthonormal constraint ( ) through Lagrangian function construction method, the constrained optimization problem in equation (14) can be structured in the following form: Combine the known conditions and do the transformation to equation (17), then take the variation with respect to ( ) The necessary condition of extreme value for this functional problem is 0 J can be an arbitrary function, the below equation can be developed.
the correlation function between two points in the space.
Considering the mutual independence of value for 1, , , to make sure that (20) is always equal to zero, if and only if (21) is true. This is the necessary condition of extremum for (14).

Snapshots Method
In this way, the original optimization problem is transformed to the solution of integral equation (21). However, in normal conditions, this equation is hard to solve. When the number of time samples is less than space samples, the snapshots method [17] is a good way to solve (21). This method can significantly reduce calculation quantities. The snapshot is obtained from sampling for each position point in discrete time. For an example, ( ) When time nodes are less than space nodes, suppose the spatial basis function ( ) , , i x y z ϕ can be expressed as a linear combination of snapshots as follows: For eigenvalue problem of equation (21), below equation can be got. , , , ,  , ,  , ,   1  , , ,1 ,  , , , 2 , ,  , , , , Substitute (22) into (23), at the same time, and let In this way, N N × eigenvalue problem in equation (21)  The dimension of the model largely determines the modelling accuracy and complexity. When deciding the number of the dimension n , the eigenvalues are arranged in descending order. The eigenfunctions corresponding to the first several eigenvalues are the most advanced representative in the system output. So the proper number of spatial basis functions can be selected. Generally speaking, the larger the number of dimension n is, the more accurate the accuracy of modelling will be. However, the modelling complexity will increase with the increasing model order. As long as the proper order of model is chose, the modelling complexity and accuracy can be guaranteed, and the ideal mathematical model would be achieved.
Define that the total energy of system is equal to the sum of all eigenvalues. Then the bigger the value is, the more energy the basis function reflects. For every eigenvalue i λ , the proportion that the energy it reflects account for the total energy is 1 .
Generally speaking, when the proportion is more than 99% , we think these bases can reflect most energy, the corresponding n is the needed dimension. Then the spatial basis functions and temporal coefficients are obtained.

Dynamic Recurrent Wavelet Neural Network Modelling Based on SNPRP Conjugate Gradient Method
In order to establish the relationship between the temporal coefficients ( ) T t and the input ( ) u t , suppose that the dynamic between ( ) T t and ( ) u t can be expressed by a nonlinear autoregressive with exogenous input (NARX) model [18]: where u n and T n are the maximum input and output lags, respectively. The modelling process of NARX model which is black-box model doesn't need to know its internal mechanism, and the mathematical model can be established based on its input and output data. Because the neural network has good nonlinear approximation capability, powerful operation ability, strong fault tolerance and robustness, so it is widely used in the field of nonlinear system identification.
Dynamic neural networks can directly reflect the system dynamic characteristic, and is widely used in the dynamic system modelling. To improve the capability of dynamic modelling, in this paper, the feedback layers for dynamic recurrent wavelet neural network is added and different weights for them in chronological order are set. A new algorithm which called dynamic memory feedback wavelet neural network (DRWNN) based on the above-mentioned idea is proposed. Considering the weakness of the traditional gradient descent method for the neural networks training, the SNPRP conjugate gradient method is used to optimize the weights of DRWNN.

The Network Structure
Considering the conventional dynamic neural networks can only memory the output of last hidden layers, by adding the multi-layer of feedback to DRWNN, it can memory several dynamic feedback of hidden layers. In the structure of DRWNN, the number of feedback layers reflects the memory of the historical data, as shown in Figure 1. ( ) n c x t R ∈ is the output of feedback layers.
Through adding Q feedback layers to the hidden layers to dynamically memory the output of hidden layers. Suppose the data of the qth feedback layer is ( ) H k q − , and its weight is . In order to make the network achieve dynamic oblivion function, iq wb is set as q i wb . Because of the weights ranging from 0 to 1, q i wb will be smaller with time which can gradually forget further information. Suppose α is the feedback gain, the output of feedback layer can be computed by the following formula: The dynamic equations of DRWNN can be expressed as:

The SNPRP Conjugate Gradient Method
The SNPRP conjugate gradient method is a new conjugate gradient method which is presented in [16]. It has the sufficient descent property and global convergence property under the condition of strong line search. Furthermore, it needs less stored information which ensures the high computational efficiency. Then a brief introduction to this method is given as follows.
Considering a general unconstrained optimization problem is a one order continuously differentiable nonlinear function, and its gradient vector is ( ) ( ) For a general spectral conjugate gradient method [23], below equation is satisfied. 1 , where d k is the search direction, k β is the parameter, and k η is the step-length factor. The spectral coefficient is defined as In SNPRP conjugate gradient method, the parameter k β is defined as The corresponding convergence and stability has been proved in [16].
The algorithm process is as follows. Then the SNPRP conjugate gradient method is used to train the weights of DRWNN.  The training goal of DRWNN is to adjust the parameters to decrease E gradually. From the chain rule of neural network on the basis of SNPRP conjugate method, the weights can be inferred as follows,

The Training Algorithm of DRWNN
is the learning coefficients of 1 2 , , , , For the equations above, The initial training value is With the method above, the final DRWNN is determined.

Experimental Results and Analysis
In order to verify the modelling ability of the proposed DRWNN, A single input single output dynamical nonlinear system is used for the testing. The dynamic equation is described as: In this section, the neural network used for testing contains 6 hidden layers, 1 input layer, 1 output layer, and the number of feedback layers of DRWNN is 3. ( ) u t and ( ) y t of the first 200 time steps are taken as the training samples. The training signal is shown in Figure 2. In order to compare the performance, the RWNN is introduced to model with the same training way [11]. The training number of both RWNN and DRWNN is 200 times.
where * i y is the real value, i y is the calculated value, and N is the number of samples. The training results are shown as Table 1. It can be shown that no matter for the maximal absolute error, the mean error or the root-mean-square error, the results of DRWNN is always less than that of RWNN, and the model output curve of DRWNN is closer than that of RWNN. So DRWNN has better modelling ability. In order to verify the generalization ability of DRWNN, the input signal ( ) u t from the time step 201 to 400 is served as the testing signal ( ) ɶ u t , and the output ( ) ɶ y t is obtained. The training signal is shown in Figure 3 and the corresponding errors are shown in Table 2. These results demonstrate that the output of DRWNN is much closer to the real output, which verifies the DRWNN has better model generalization ability. Using the SNPRP conjugate gradient method to optimize the weights of DRWNN is effective.

Modelling for ASP Flooding with
Proposed Spatial-Temporal Separation Method

Reservoir Description
Suppose that the reservoir of ASP flooding consists of four injections and nine production wells. All wells distribute uniformly; there is one injection well at the center of every four production wells. The distribution of wells is shown in Figure 4. As for this reservoir, the length is 630m; the width is 630m, and the thickness is 19.990m. There are 7 layers in all; the thickness of each layer is 2.857m; and the net thickness is 1.4286m; the depth of upper surface is 2420m; the porosity of every layer is 0.3, and the pore volume is 1.1097 × 10 6 m 3 . The initial grid concentration of ASP is 0 g/L. The initial water saturation is shown in Figure 5. The grids of reservoir are divided in three directions , , x y z . The grids in x and y are divided into 21 units, respectively, and the grids in z are divided into 7 units. The total number is 21×21×7. The water injection rate of each injection well is 83 m 3 /d. The rate of production wells is as defined in Table 3. Using reservoir simulation software CMG, the whole oil production process can be obtained by simulation. The whole production time lasts for 96 months, and sampled output can be got by the snapshots method. The injecting time of ASP is served as the initial time, so the total time nodes are =97 L and the total space nodes are 21 21 7 N = × × . Table 3. The liquid volume quantity of production wells. wells S1-1 S1-2 S1-3 S2-1 S2-2 S2-3 S3-1 S3

Modelling and Verification for ASP Flooding
Given the ASP flooding above, the approximate model is built with the method proposed in this paper. In order to sufficiently motivate system, the simulation is run on the software CMG for 50 times. The software randomly generates different injection strategies for every month, and the process of injection totally lasts for 48 months. Then the three-dimensional spatial-temporal modelling method is used to get the spatial basis functions and the temporal coefficients of the 50 sets of grid saturation. So as to get the best structure to represent the ASP flooding system, use one group of the spatial basis functions to reconstitute the system with other 49 groups of temporal coefficients and calculate the average mean square errors of all sampling points for these 49 times. After executing this process for 50 times, a group of spatial basis functions with the smallest mean square error is served as the spatial basis functions for modelling. The relationship between injection concentrations of ASP and temporal coefficients is get by DRWNN, then combine the spatial basis functions to build the model between grid water saturation and injection concentrations of ASP. At the same time, build the model between the grid water saturation and moisture content of production wells with DRWNN. The mathematical form can be described as follows, , , where w n is the lag of moisture content, s n is the lag of grid water saturation.
Then the model between the injection concentration of ASP and moisture content of production wells can be built after integrating the models above, which can be expressed as equations (13), (28) and (50).
In order to compare the influence of the number of the spatial basis functions, the performance indictor is defined as follows: Generally speaking, the more the number of the spatial basis function is, the more information of system the model reflects, and the higher accuracy the model has. But at the same time, the model will be sensitive to the external disturbance and the generalization ability will decrease. Besides, the dimension of the model increases. On the other hand, the less quantity the spatial basis function is; the model will reflect less information and the error can be very big. It is important to choose the proper amounts of the spatial basis functions. In order to test the influence of different numbers spatial basis function on the accuracy of modelling, different numbers of spatial basis functions is chose to reconstruct the system with the temporal coefficients and calculate the average RMSE which is shown in Table 4. From Table 4, the accuracy of model increase with the number of spatial basis functions. However, when the number is over 4, the RMSE of model will not increase which indicate that the spatial basis functions can reflect the whole reservoir, so the number of spatial basis functions adopted in this paper is 4. Besides that, the number of feedback layers of the DRWNN networks for injection concentrations and temporal coefficients is 3; the corresponding hidden layers is 6; the lags of temporal coefficients is 4 T n = ; the lags of injection concentrations is 3 u n = . The number of feedback layers of DRWNN networks for moisture content of production wells and the corresponding grid water saturation is 2; the hidden layers is 8; the lags of moisture content is 3 w n = ; the lags of grid water saturation is 2 s n = . Model for the ASP flooding with the spatial-temporal separation method based on DRWNN which is proposed in this paper. The errors for this model are that the mean RMSE for grid water saturation is 0.0285, and the mean RMSE for moisture content of production wells is 1.1378%. This demonstrates the good modelling ability.
In order to verify the generalization ability, the injection concentrations of ASP flooding are given randomly:

( )
2.9,1.3, 0.6 S u = kg/m 3 . The displacing agent injection totally lasts for 48 months which are divided into 3 slugs uniformly, and the rest time is water flooding. The specific injection can be found in Figure 6-a. Because the whole system is five-dimensional considering the time and value of water saturation, it can't be plotted in the picture directly. In order to illustrate clearly, the 21 grids ranging from ( ) The injection mass concentration The output of real system The output of modelling To analyze the generalization ability further, we compute the error with equation (48). The mean absolute errors for the moisture content of production wells and for the grid water saturation are 1.2970% and 0.0318% , respectively. Figure 8 compares the moisture content of modelling and that of the Water Saturation simulation software. It can be known that the error is very small, which verifies that the whole model has better generalization ability.
From the above, the three-dimensional spatial-temporal separation modelling method based on DRWNN can model for ASP flooding well. It has good modelling accuracy and generalization ability.

Conclusions
In this paper, a three-dimensional spatial-temporal separation modelling approach is proposed to establish the identification model for ASP flooding. At first, spatial state (the water saturation) is expanded onto a series of dominant spatial basis functions and temporal coefficients. Then a new dynamic recurrent wavelet neural network, in which feedback layers are added and different weights are set with time to achieve dynamic memory of the past information, is proposed to build the relationship between the temporal coefficients and the injection concentrations of ASP flooding. To avoid the local minimum and low convergence, the SNPRP conjugate gradient method is adopted to train the weights of DRWNN. Besides that, DRWNN is also used to establish the model of the moisture content of production wells and the corresponding grid water saturation. The simulation of ASP flooding is carried out to show the effectiveness of this spatial-temporal modelling method. The result shows that the method proposed in this paper has good accuracy and generalization ability. It is suitable to be used for modelling for complex distributed parameter systems like ASP flooding.

D i w o j s p
∈ ∈ is the diffusion coefficient of component j in solution i .