A general numerical model has been developed to predict the probability density function (pdf) of the magnitudes of complex pulse-compressed broadband echoes (broadband echo pdf) due to arbitrary aggregations of scatterers that are detected with a single-beam echosounder. The model is based on physics principles and rigorously accounts for the broadband frequency-dependent characteristics of the system, signal, scatterers, and the beampattern modulation effects of the sonar/radar transceiver. A key aspect to modeling the statistics of broadband echoes is accounting for the scenarios where the pulse-compressed echoes may only partially overlap. The echo statistics under those conditions will be significantly different than those associated with a narrowband system with the same center frequency but whose echoes will completely overlap given the same density of scatterers. As a result, pulse-compressed broadband echoes will generally deviate more from the Rayleigh distribution (i.e., be more “heavy-tailed”) than narrowband echoes-a feature that is critical to analyzing real-world data. As in the case of narrowband signals, the shape of the broadband echo pdf is shown to vary from strongly non-Rayleigh to Rayleigh as the number of dominant scatterers in the beam increases. The model is applied to sonar in the ocean in which the numerical density of fish is inferred using broadband echoes (30-70 kHz). The results are compared with those from conventional echo energy methods. As with narrowband systems, statistics of broadband echoes can be used to estimate the numerical density of scatterers without the need for absolute calibration of the system.