Natural phenomena indicate there exist cases when fracture characteristics (e.g., spacing, trace length, dip direction, and dip) have apparent spatial patterns rather than being independent of each other. This study proposes a geostatistical approach to incorporate spatial variability of fracture characteristics into fracture network modeling. Results indicate the proposed approach is valid, robust, flexible, and able to generate fracture networks with given statistics (i.e., the mean
μ, the standard deviation
σ, and the correlation range
λ) of spatial variability of fracture characteristics. Meanwhile, not only the autocorrelation of fracture characteristics, but also the cross-correlation of different fracture sets is incorporated into the fracture networks. Furthermore, our approach is successfully applied to generate fracture networks with certain rock mass structures such as the conjugate structure, the blocky structure, the finely laminated structure, the interlayer structure, and the cataclastic structure, even independent, nonperiodic or purely random cases. These results demonstrate the proposed approach facilitate the fracture network modeling to better meet needs from different engineering practices by generating fracture networks with given characteristics.