Fluid flow in unpropped and natural fractures is critical in many geophysical processes and engineering applications. The flow conductivity in these fractures depends on their closure under stress, which is a complicated mechanical process that is challenging to model. The challenges come from the deformation interaction and the close coupling among the fracture geometry, pressure, and deformation, making the closure computationally expensive to describe. Hence, most of the previous models either use a small grid system or disregard deformation interaction or plastic deformation.

In this study, a numerical model is developed to simulate the stress-driven closure and the conductivity for fractures with rough surfaces. The model integrates elastoplastic deformation and deformation interaction, and can handle contact between heterogeneous surfaces. Computation is optimized and accelerated by use of an algorithm that combines the conjugate-gradient (CG) method and the fast-Fourier-transform (FFT) technique. Computation time is significantly reduced compared with traditional methods. For example, a speedup of five orders of magnitude is obtained for a grid size of 512 × 512. The model is validated against analytical problems and experiments, for both elastic-only and elastoplastic scenarios.

It is shown that interaction between asperities and plastic deformation cannot be ignored when modeling fracture closure. By applying our model, roughness and yield stress are found to have a larger effect on fracture closure and compliance than Young’s modulus. Plastic deformation is a dominant contributor to closure and can make up more than 70% of the total closure in some shales. The plastic deformation also significantly alters the relationship between fracture stiffness and conductivity. Surfaces with reduced correlation length produce greater conductivity because of their larger apertures, despite more fracture closure. They have a similar fraction of area in contact as compared with surfaces with longer fracture length, but the pattern of area in contact is more scattered. Contact between heterogeneous surfaces with more soft minerals leads to increased plastic deformation and fracture closure, and results in lower fracture conductivity. Fracture compliance appears not to be as sensitive to the distribution pattern of hard and soft minerals. Our model compares well with experimental data for fracture closure, and can be applied to unpropped or natural fractures. These results are obtained for a wide range of conditions: surface profile following Gaussian distribution with correlation length of 50 µm and roughness of 4 to 50 µm, yield stress of 100 to 1500 MPa, and Young’s modulus of 20 to 60 GPa. The results may be different for situations outside this range of parameters.