Obtaining information of the spatial distribution of subsurface natural and induced fractures is critical in improving the production of geothermal or hydrocarbon fluids. Traditional seismic characterization methods for subsurface fractures are usually based on the effective anisotropy medium theory, which may not be true in reality where the fracture distribution is non-uniform. In this abstract, we propose to test the double-beam method to characterize non-uniformly distributed fractures that are commonly observed in the unconventional reservoirs. We built a 3D layered reservoir model and the reservoir layer is geometrically irregular and it contains a set of randomly spaced fractures with spatially varying fracture compliances. We used an elastic full-wave finite-difference method to model the wavefield where we treat the fractures as linear-slip boundaries and the recorded data include all elastic multiple scattering. Taking the surface seismic data as input, the double-beam method forms a focusing source beam and a focusing receiver beam toward the fracture target. The fracture information is derived from the interference pattern of these two beams, which gives fracture orientation, fracture spacing, and fracture compliance as a function of spatial location. The fracture orientation parameter is the most readily determined parameter. The beam interference amplitude depends on both fracture spacing and compliance in a local average sense for random fractures. The beam interference amplitude is large when there are dense fractures or the compliance value is large, which is important in the interpretation of the fluid transport properties of a reservoir.