Because of its low expense and high effectiveness, waterflooding is the most popular enhanced oil recovery technique applied in Chinese oil fields, such as Daqing oil field and Shengli oil field. During the waterflooding process, water is injected into the formation via the injection wells continually and the injected water will make great impacts on the formation parameters, such as permeability, porosity, content of clay minerals and contact pattern of skeleton particles. Consequently, these impacts could be of great importance for the distribution of the remaining oil and development of the oil fields (Han, 2010; Shokri & Babadagli, 2016) .
Usually, researches about the formation parameters variation are carried out from the macroscopic perspective and the main techniques are the indoor displacement experiment and the core analysis (You et al., 2007; Li, 2005; Wu et al., 2002; He & Xu, 2010; Shi et al., 2013) . For indoor displacement experiment, rock cores obtained from the field are used and the water displacement process is applied to them. During the displacement process, experimental data of the permeability and porosity at different injected pore volume are recorded and curves based on the experiment data will be made for the analysis of parameters variation. For core analysis, coring is conducted and cores at various water-cut stages are gained. Then the parameters of these cores could be measured and analyzed. However, results of these methods are just the reflection of the variation rules of formation parameters and the reasons and mechanisms of these variations are still unknown.
For its advantages of low expense and dynamic prediction, reservoir simulation is quite a significant technique used during the long period of reservoir development. Besides, it is also quite useful for the study of the remaining oil distribution (Rege & Fogler, 1987; Jiang et al., 2005; Cui et al., 2012; Feng et al., 2009; Crandell et al., 2012; Feng & Bai, 2011) . However, formation permeability in the traditional reservoir simulators is set as a constant which is quite different from that of the actual conditions and this will result in the inaccuracy of predictive development results.
Based on the analysis above, several works are conducted here. A network model for formation parameter variation is introduced by taking various particle mechanisms into consideration, such as particle detachment, particle deposition. Variation of the formation permeability and porosity is studied and a mathematical model is proposed based on the simulation results. Then the coupled mathematical model for formation parameters variation is established and simulations under different conditions are conducted based on a novel reservoir numerical simulator. The reservoir permeability variation and its influence on the remaining oil distribution are analyzed.
2. Pore Network Modeling for Formation Parameters Variation
Pore network modeling was firstly proposed by Fatt in 1950s to study microscale multiphase flow. In recent years, there has been increasingly interest in pore-scale modeling for researchers from various fields and great progress has been made. Pore network modeling is no longer limited to the simple two-phase flow and the computation of relative permeability. For two-phase flow, effects of wettability, wetting hysteresis and mass transfer between phases could also be studied. Besides, some more applications of pore network modeling have been developed, such as modeling of three-phase flow, non-Newtonian flow and formation impairment (Hou et al., 2005; Blunt et al., 2013; Feng et al., 2015; Li et al., 2017; Watson et al., 2017) .
In this part, a pore network model of formation parameters variation is proposed by taking into account different mechanisms about particles. The main mechanisms considered in this model are detachment and entrapment. In addition, pores and throats of the rock are represented by cylinders in this model and have circular cross-section.
2.1. Detachment of Particles from Pore/Throat
As analyzed above, flow of injected water will bring drag force on the particles attached on the surface of pore/throat. And if the fluid velocity goes higher, particles will be detached from the surface and flow out of the reservoir with injected fluid. Then the pore/throat size will increase, which will be in favor of fluid flowing through the reservoir. The mathematical model presented by Jalel & Jean-Francois (1999) is employed to calculate the detachment rate of particles from the pore surface,
where, is the detachment rate of particles from the pore/throat surface per unit area, ; is the rate of fluid flowing in pore space, m/s; is the critical flowing rate and the particles just begin to detach when the fluid flowing rate is higher than the critical rate, m/s; is the release coefficient, and equals to zero when ; is the volumetric concentration of particles on the pore/throat surface, dimensionless.
2.2. Entrapment of Particles to Pore/Throat Surface
Particles flowing with injected water may also be precipitated and attached on the surface of pore/throat again or strained at the pore/throat entrance which will result in the blockage of pore/throat. The main entrapment mechanisms include surface deposition, direct blockage and bridging.
Surface deposition means that detached particles with smaller size may be reattached to the pore/throat surface under the effect of gravity or electric forces and result in the decrease of pore/throat radius. The equation for describing particle deposition in the pore network simulation was proposed by Jalel & Jean-Francois (1999) ,
where, is the deposition rate per unit pore surface area, ; Cfi is the volumetric particle concentration of the fluid in pore space, dimensionless; is the radius of pore/throat, m; is the fluid flowing rate through the pore/throat, m/s; is the particle radius, m; is the length of throat, m.
Particle blockage happens when particles flow through a pore or throat whose radius is smaller than the particles. Then the particle will be entrapped at the entrance and make the pore/throat blocked.
Particle bridging means that several particles whose radius is smaller than that of their passing pore/throat could plug the pore/throat by bridging at its entrance. Generally, scholars consider that when the size of the particle is bigger than 1/3 of its passing pore/throat, bridging would take place and result in the blockage of the pore/throat (Faruk, 2010) .
Computation algorithm for the pore network modeling is shown in Figure 1. Simulations are conducted to study formation parameters variation under various pressure gradients. The simulation results of permeability and porosity variation are shown in Figure 2 and Figure 3.
As shown in Figure 1, firstly, set the initial conditions, including model parameters, boundary conditions, and fluid parameters. Secondly, construct the pore network for formation parameters variation. Thirdly, based on the pore network model, calculate the formation permeability and porosity with the iterative method.
As shown in Figure 2, for a certain pressure gradient, permeability of the network increases gradually and the curve goes stable after long-term water flooding. For simulations at different pressure gradient, the permeability increases more quickly at higher pressure gradient. And Figure 3 reveals that the porosity variation principle is similar to that of permeability.
Figure 1. Computation algorithm for the pore network modeling.
Figure 2. Simulation results of permeability variation.
Figure 3. Simulation results of porosity variation.
3. The Coupled Mathematical Model for Formation Parameters Variation
3.1. The Coupled Mathematical Model
Based on the pore network modeling for formation parameters variation, the variation curves of reservoir permeability under different conditions are shown in Figure 4. The mathematical model for formation permeability variation is obtained by multiple linear regressions.
Figure 4. Variation curves of reservoir permeability.
where, k, k0, kmax mean the absolute permeability at every moment, the initial absolute permeability, and the maximum absolute permeability, respectively. μm2; Q, Qmax mean the cumulative flowing rate per unit pore/throat area, the maximum cumulative flowing rate per unit pore/throat area, respectively, m/s; uc means the cementation strength among sand particles, which is classified into three groups (uc = 1, 2, 3), dimensionless, and uc = 3 indicates the minimum cementation strength. Vsh means the clay content, dimensionless.
The basic mathematical model for fluids flowing in porous media is,
where, so, sg, sw mean the saturation of oil phase, gas phase, water phase, respectively; ρo, ρg, ρw mean the density of oil phase, gas phase, water phase, respectively; Rso, Rsw mean the dissolved gas-oil ratio, the dissolved gas-water ratio, respectively; Bo, Bg, Bw mean the formation volume factor of oil phase, gas phase, water phase, respectively; kro, krg, krw mean the relative permeability of oil phase, gas phase, water phase, respectively; μo, μg, μw mean the viscosity of oil phase, gas phase, water phase, respectively; Po, Pg, Pw mean the pressure of oil phase, gas phase, water phase, respectively; qvo, qvg, qvw mean the injection/production rate of oil, gas, water, respectively.
Combined Equations (3)-(5) with Equation (6), the formation permeability variation will be taken into consideration in the mathematical model for fluids flowing in porous media, i.e., the coupled mathematical model for formation parameters variation is established.
3.2. Numerical Solution
The coupled mathematical model for formation parameters variation consists of a set of nonlinear equations, which mainly includes the continuity equations of oil, water, and gas phase, the equations of formation permeability variation, and a series of auxiliary equations. The finite difference method is used to solve the nonlinear equation system since the analytical solutions of the system are intractable. In this work, the implicit pressure and explicit saturation method (IMPES) and Gauss-Seidel iterative technique are used to solve the coupled mathematical model. The procedures of the solution are listed as follows:
1) The pressure and saturation of oil, water, and gas are obtained first by IMPES.
2) The cumulative flowing rate per unit pore/throat is calculated by Darcy’s law.
3) The new absolute permeability (k) is calculated by formation permeability variation model (Equations (3)-(5)).
4) The next iteration is computed if maximum time is not reached.
4. Numerical Results and Discussions
Based on the coupled mathematical model for formation parameters variation, a novel three dimensional (3D) field-scale reservoir numerical simulator is developed in Fortran 90. A conceptual model is simplified from a typical well group in a China offshore oilfield. The geological model parameters are given in Table 1 and Table 2, and the relative permeability curves are shown in Figure 5. On the basis of the novel 3D field-scale reservoir numerical simulator, the permeability and remaining oil distribution is studied in reverse rhythm reservoir with inverted nine-spot pattern.
4.1. Comparisons of the Permeability and Remaining Oil Distribution with and without Considering Formation Parameters Variation
Take the conceptual model with permeability ratio of 5 for example. When the water cut reaches 90%, the permeability distribution and remaining oil distribution of the 10th layer is shown in Figure 6. However, when the variation of
Table 1. The geological model description.
Table 2. Reservoir initial horizontal permeability.
Figure 5. The relative permeability curves.
Figure 6. The permeability and remaining oil distribution of the 10th layer considering the formation parameters variation. (a) Permeability distribution; (b) Remaining oil distribution.
formation parameters is not taken into consideration in the reservoir numerical simulator, the permeability distribution and remaining oil distribution of the 10th layer is shown in Figure 7.
Compared Figure 6 with Figure 7, when the variation of formation parameters is taken into consideration, the reservoir permeability changes obviously around the injection well and the main streamlines between the injection and production wells, and the remaining oil saturation, as is shown in Figure 6, is higher than that in Figure 7. The results can be validated by the experiment data carried out by Chen et al. (2016) . In the experiment, the initial formation permeability is 5661.9 × 10−3 μm2. After long-term water flooding, the formation permeability increases to 7574.5 × 10−3 μm2. For the reverse rhythm reservoir with large permeability ratio, the injection water flow to production wells along the top layers preferentially. Therefore, the high-permeability channels are more likely to develop in the top layers and the remaining oil mainly accumulates in the bottom layers.
4.2. The Sensitivity Analysis for Permeability Ratio
Take the conceptual model with permeability ratio of 1 and 5 for example. When the water cut reaches 90%, the permeability distribution and remaining oil distribution of the 2nd and the 12th layer are shown in Figure 8 and Figure 9 respectively.
Compared Figure 8 with Figure 9, the reservoir permeability of the 2nd layer changes greater and the remaining oil saturation of the 12th layer is higher when the permeability ratio is 5. For the reverse rhythm reservoir with larger permeability ratio, the injection water flows to the production well along the top layers more easily, and the erosion effect of the injection water on the top layers is stronger. Therefore, the high-permeability channels evolve seriously in the top layers and the remaining oil mainly accumulates in the bottom layers. The results can be validated by the experiment data carried out by Liu et al. (2014) . In
Figure 7. The permeability and remaining oil distribution of the 10th layer without considering the formation parameters variation. (a) Permeability distribution; (b) Remaining oil distribution.
Figure 8. The permeability and remaining oil distribution when permeability ratio is 1. (a) Permeability distribution of the 2nd layer; (b) Remaining oil distribution of the 12th layer.
Figure 9. The permeability and remaining oil distribution when permeability ratio is 5. (a) Permeability distribution of the 2nd layer; (b) Remaining oil distribution of the 12th layer.
the experiment, the permeability of high permeability ratio formation increases due to the injected water during the later stage of production. The high-permeability channels are more obvious influenced by water flooding and later physical and chemical reformation.
The novel simulation framework combining microcosmic and macroscopic mechanism provides considerable guiding significance for predicting the formation parameters variation. What’s more, the study also demonstrates that the novel simulation framework provides sufficient information for remaining oil description.
But the formation parameters variation is difficult to determine. Therefore, it can only be approximately simulated through the novel simulation framework. To obtain more precise results, a series of experiments on formation parameters variation are necessary to be carried out.
 Chen, D. Q., Li, J. Y., Zhu, W. S., & Xin, Z. (2016). Experimental Re-search on Reservoir Parameters Variation after Water Flooding for Offshore Unconsolidated Sandstone Heavy Oil Reservoirs. China Offshore Oil and Gas, 28, 64-60.
 Crandell, L. E., Peters, C. A., Um, W., Jones, K. W., & Lindquist, W. B. (2012). Changes in the Pore Network Structure of Hanford Sediment after Reaction with Caustic Tank Wastes. Journal of Contaminant Hydrology, 131, 89-99.
 Cui, C. Z., Geng, Z. L., Wang, Y. Z., Huang, Y., & Liu, H. (2012). Calculation Model of Dynamic Permeability Distribution and its Application to Water Drive Reservoir at High Water Cut Stage. Journal of China University of Petroleum (Edition of Natural Science), 36, 118-122.
 Faruk, C. (2010). Non-Isothermal Permeability Im-pairment by Fines Migration and Deposition in Porous Media Including Dispersive Transport. Transport in Porous Media, 85, 223-258.
 Feng, Q. H., Li, S., Han, X. D. et al. (2015). Network Simulation for Formation Impairment due to Suspended Particles in Injected Water. Journal of Petroleum Science and Engineering, 133, 384-391.
 Feng, Q. H., Qi, J. L., Yin, X. M. et al. (2009). Simulation of Fluid-Solid Cou-pling during Formation and Evolution of High-Permeability Channels. Petroleum Exploration and Development, 36, 498-502.
 Jalel, O., & Jean-Francois, V. (1999). A Two-Dimensional Network Model to Simulate Permeability Decrease under Hydrodynamic Effect of Particle Release and Capture. Transport in Porous Media, 37, 303-325.
 Jiang, H. Q., Gu, J. W., Chen, M. F., & Sun, M. (2005). Reservoir Simulation of Remaining Oil Distribution Based on Time-Variant Reservoir Model. Petroleum Explo-ration and Development, 32, 91-93.
 Li, J., McDougall, S. R., & Sorbie, K. S. (2017). Dynamic Pore-Scale Network Model (PNM) of Water Imbibition in Porous Media. Advances in Water Resources, 107, 191-211.
 Liu, H. Q., Li, B. Y., Wang, W. F. et al. (2014). Characteristics Analysis of the Reservoir and Predominate Channel in Block Jin-16. Journal of Southwest Petroleum University (Sci-ence & Technology Edition), 36, 60-68.
 Rege, S. D., & Fogler, H. S. (1987). Network Model for Straining Dominated Particle Entrapment in Porous Media. Chemical Engineering Science, 42, 1553-1564.
 Shi, C. L., Zhang, F. H., & Chen, P. (2013). Affection of Sim-ulating Water-flooding by Water Injection Tests on Reservoir Properties. Journal of Southwest Petroleum University (Science & Tech-nology Edition), 35, 87-93.
 Shokri, A. R., & Babadagli, T. (2016). Field Scale Modeling of CHOPS and Solvent/Thermal Based Post CHOPS EOR Applications Considering Non-equilibrium Foamy Oil Behavior and Realistic Representation of Wormholes. Jour-nal of Petroleum Science and Engineering, 137, 144-156.
 Watson, M. G., Bondino, I., Hamon, G., & McDougall, S. R. (2017). A Pore-Scale Investigation of Low-Salinity Waterflooding in Porous Media: Uniformly Wetted Systems. Transport in Porous Media, 118, 201-223.
 You, Q. D., Zhou, F. X., Zhang, J. L., & Chen, Y. (2007). Law and Mechanism of Parameters Change in Saliferous Reservoir. Journal of China University of Petroleum (Edition of Natural Science), 31, 79-82.