As common physical phenomena in porous media,capillarity behaviors exist in many engineering applications and natural science fields. The experimental,theoretical and numerical research on capillarity in porous media has lasted for more than a century,and the research results have been widely used in various fields,such as the development of conventional and unconventional resources. However,although the research has made great progress,the complex imbibition mechanism poses new challenges to us. The 1st National Conference on Imbibition Theory and Application in Porous Media was held in Beijing from April 22 to 24,2023, to gather researchers who are interested in imbibition research,exchange the latest progress and achievements in the field of imbibition in porous media,and discuss research hotspots and difficulties.

Acidization is a widely used stimulation technique for carbonate reservoirs aimed at removing formation damage, and if successful, can result in the creation of wormholes of specific lengths and conductivities around the wellbore. The formation of wormholes depends on the injection rate for a particular acid-mineral system and can be predicted through numerical simulations of the reactive phenomenon during acidization. In this paper, the commonly used two-scale continuum model is enhanced to encompass fractured-vuggy porous media. The fractures are characterized by a pseudo-fracture model, while vugs are represented by a cluster of anomalous matrices with high porosity. Moreover, a method for generating random pore-fracture-vuggy models is proposed. The governing equations are discretized by the finite volume method and are solved under three-dimensional linear and radial conditions. Sensitivity analysis of dissolution dynamics with respect to fracture and vug parameters is performed. The simulation results indicate that both fractures and vugs significantly impact wormhole development. Except for fractures perpendicular to the acid flow direction, fractures in other directions play a crucial role in determining the direction of wormhole growth.

In the last decades, shale gas development has relieved the global energy crisis and slowed global warming problems. The water bridge plays an important role in the process of shale gas diffusion, but the stability of the water bridge in the shale nanochannel has not been revealed. In this work, the molecular dynamics method is applied to study the interaction between shale gas and water bridge, and the stability can be tested accordingly. CO₂ can diffuse into the liquid H₂O phase, but CH₄ only diffuses at the boundary of the H₂O phase. Due to the polarity of H₂O molecules, the water bridge presents the wetting condition according to model snapshots and one-dimensional analyses, but the main body of the water bridge in the two-dimensional contour shows the non-wetting condition, which is reasonable. Due to the effect of the molecular polarity, CO₂ prefers to diffuse into kerogen matrixes and the bulk phase of water bridge. In the bulk of the water bridge, where the interaction is weaker, CO₂ has a lower energy state, implies that it has a good solubility in the liquid H₂O phase. Higher temperature does not facilitate the diffusion of CO₂ molecules, and higher pressure brings more CO₂ molecules and enhances the solubility of CO₂ in the H₂O phase, in addition, a larger ratio of CO₂ increases its content, which does the same effects with higher pressures. The stability of the water bridge is disturbed by diffused CO₂, and its waist is the weakest position by the potential energy distribution.

In order to make full use of the information provided in the physical reservoirs, including the production history and environmental conditions, the whole life cycle of reservoir discovery and recovery should be considered when mapping in the virtual space. A new concept of reservoir digital twin and the exploratory multi-scale framework is proposed in this paper, covering a wide range of engineering processes related with the reservoirs, including the drainage, sorption and phase change in the reservoirs, as well as extended processes like injection, transportation and on-field processing. The mathematical tool package for constructing the numerical description in the digital space for various engineering processes in the physical space is equipped with certain advanced models and algorithms developed by ourselves. For a macroscopic flow problem, we can model it either in the Navier-Stokes scheme, suitable for the injection, transportation and oil-water separation processes, or in the Darcy scheme, suitable for the drainage and sorption processes. Lattice Boltzmann method can also be developed as a special discretization of the Navier-Stokes scheme, which is easy to be coupled with multiple distributions, for example, temperature field, and a rigorous Chapman-Enskog expansion is performed to show the equivalence between the lattice Bhatnagar-Gross-Krook formulation and the corresponding Navier-Stokes equations and other macroscopic models. Based on the mathematical toolpackage, for various practical applications in petroleum engineering related with reservoirs, we can always find the suitable numerical tools to construct a digital twin to simulate the operations, design the facilities and optimize the processes.

In this work the problem of displacing a ganglion of a fluid by another immiscible one in capillaries is investigated. A modeling approach is developed to predict the location of the ganglion with time. The model describes two patterns; namely, when the ganglion totally exists inside the tube, and when the advancing interface of the ganglion has broken through the exit of the tube. The model is valid for the case in which the ganglion is wetting as well as when it is nonwetting to the wall of the tube. It also considers the situation in which both the advancing and the receding interfaces assume, generally, different contact angles. For the special case when the displacement process is quasistatic, both the receding and the advancing contact angles may be considered the same. Under these conditions, interfacial tension force plays no role and the ganglion moves as a plug inside the tube with a constant velocity. When the viscosity ratio between the invading fluid and the ganglion is one (i.e., both phases are having the same viscosity) the motion reduces to the Hagen-Poiseuille flow in pipes. Once the advancing interface breaks through the exit of the tube, interfacial tension starts to take part in the displacement process and the ganglion starts to accelerate or decelerate according to the viscosity ratio. When the ganglion is nonwetting, interfacial tension becomes in the direction of the flow and is opposite to the flow otherwise. The model accounts for external forces such as pressure and gravity in addition to capillarity. A computational fluid dynamics analysis of this system is conducted for both types of wettability scenarios and shows very good match with the results of the developed model during both the two modes of flow patterns. This builds confidence in the developed modeling approach. Other cases have also been explored to highlight the effects of other scenarios.

Water uptake in rock fractures caused by rainfall plays a significant role in slope stability analysis. Since the fracture network system has complicated structures and multiple scales, the models based on the averaged system cannot account for these properties. On the other hand, a model describing a single fracture with fractal characteristics and surface roughness fails to deal with the case of multiple fractures at spatial scales. In this study, a fracture-network model is established to account for the complex structures and multiple scales of fractures. By considering the connectivity between fractures and the limited area of aquifer, capillary pressure formulations in different fractures are derived based on the Young-Laplace equation, and the final water level under specific rainfall conditions is also obtained. The cross-section shapes and exhaust conditions of rainwater infiltration have important influences on the final water level. The results indicate that the final water level is proportional to the ratio of perimeter to cross-section area when the fracture is a cylinder, and a circular pipe can reduce water level elevation in the fracture system.

Compositional multiphase flow in subsurface porous media is becoming increasingly attractive due to issues related with enhanced oil recovery, CO ${}_{\mathrm{2}}$ sequestration and the urgent need for development in unconventional oil/gas reservoirs. One key effort to construct the mathematical model governing the compositional flow is to determine the phase compositions of the fluid mixture, and then calculate other related physical properties. In this paper, recent progress on phase equilibrium calculations in unconventional reservoirs has been reviewed and concluded with authors’ own analysis, especially focusing on the special mechanisms involved. Phase equilibrium calculation is the main approach to investigate phase behaviors, which could be conducted using different variable specifications, such as the NPT flash and NVT flash. Recently, diffuse interface models, which have been proved to possess a high consistency with thermodynamic laws, have been introduced in the phase equilibrium calculation, incorporating the realistic equation of state (EOS), e.g. Peng-Robinson EOS. In the NVT flash, the Helmholtz free energy is minimized instead of the Gibbs free energy used in NPT flash, and this thermodynamic state function is decomposed into two terms using the convex-concave splitting technique. A semi-implicit numerical scheme is applied to the dynamic model, which ensures the thermodynamic stability and then preserves the fast convergence property. A positive definite coefficient matrix is designed to meet the Onsager reciprocal principle so as to keep the entropy increasing property in the presence of capillary pressure, which is required by the second law of thermodynamics. The robustness of the proposed algorithm is demonstrated by using two numerical examples, one of which has up to seven components. In the complex fluid mixture, special phenomena could be captured from the global minimum of tangent plane distance functions and the phase envelope. It can be found that the boundary between the single-phase and vapor-liquid two phase regions shifts in the presence of capillary pressure, and then the area of each region changes accordingly. Furthermore, the effect of the nanopore size distribution on the phase behavior has been analyzed and a multi-scale scheme is presented based on literature reviews. Fluid properties including swelling factor, criticality, bubble point and volumetrics have been investigated thoroughly by comparing with the bulk fluid flow in a free channel.