This work concerns modeling of very high frequency (>100 kHz) sonar images obtained from a sandy seabed. The seabed is divided into a discrete number of 1D height profiles. For each height profile the backscattered pressure is computed by an integral equation method for interface scattering between two homogeneous media as formulated by Chan [IEEE Trans. Antennas Propag. 46, 142-149 (1998)]. However, the seabed is inhomogeneous, and volume scattering is a major contributor to backscattering. The SAX99 experiments revealed that the density in the unconsolidated sediment within the first 5 mm exhibits a high spatial variation. For that reason, additional roughness is introduced: For each surface point a stochastic realization of the density along the vertical is generated, and the sediment depth at which the density has its maximum value will constitute the new height field value. The matrix of the full integral equation is reduced to a band matrix as the interaction between the point sources on the seabed is neglected from a certain range; this allows computations on long height profiles with lengths up to approximately 25 m (at 300 kHz). The equivalent roughness approach, combined with the band-matrix approach, agrees with SAX99 data at 300 kHz. © 2007 Acoustical Society of America.