Obtaining information on the spatial distribution of subsurface natural and induced fractures is critical in the production of geothermal or hydrocarbon fluids. Traditional seismic characterization methods for subsurface fractures are based on the assumption of effective anisotropy medium theory, which may not be true in reality when the fracture distribution is random. We have tested the recently proposed double-beam method to characterize nonuniformly distributed fractures. We built a 3D layered reservoir model; the reservoir layer was geometrically irregular, and it contained 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 data include all elastic multiple scattering. Taking the surface seismic data as input, the doublebeam method forms a focusing source beam and a focusing receiver beam toward the fracture target. The fracture information is derived fromthe interference pattern of these two beams, which includes fracture orientation, fracture spacing, and fracture compliance as a function of spatial location. The fracture orientation parameter is the most readily determined parameter even for multiple nonorthogonal coexistent fracture sets. The beaminterference amplitude depends on the fracture spacing and compliance in a local average sense for random fractures. The beaminterference amplitude is large when there are many fractures or the compliance value is large, which is important in the interpretation of the fluid-transport properties of a reservoir.