The conventional approach for analysis of liquid seepage in partially saturated fracture media usually employs macroscale continuum concepts
but notable unevenness of local velocity and flux in fractured rocks has been observed in the field test. Unsaturated flow in naturally fractured rocks is generally a multi-phase
non-isothermal flow that occurs in a three-dimensional fracture network. Stochastic simulation has become a common tool for characterizing and visualizing heterogeneity on incomplete information. 3D fracture systems are approximated as 2D planar fractures that are conceptualized as 2D heterogeneous porous media. Simulated annealing algorithm is applied to incorporate data from various sources into their models by formulating a suitable objective function; and then spatial fractures are characterized.