A Gaussian beam summation (GBS) formulation is introduced for a doubly rough boundary waveguide, wherein the coherent and incoherent scattered fields are decomposed into a discrete phase-space summation of Gaussian beams (GB) that emanate from the rough surfaces in all directions. The scheme involves deterministic GB propagators and stochastic GB-to-GB (GB2GB) scattering matrices for the coherent and incoherent fields, where each scattered beam is propagated inside the waveguide and is scattered again from the rough boundaries. The GB2GB matrices are calculated from the statistical moments of the scattering amplitude, which are given either analytically or empirically. An analytical and numerical example for a waveguide with weak boundary roughness is presented and discussed. The formulation reveals explicitly the phase-space footprint of the stochastic multiple scattering process at the rough boundaries, thus providing a cogent physical interpretation and an effective mathematical representation to the field. The formulation also accommodates the receiver's pattern in the same phase-space format. Bistatic reverberations inside a rough surface waveguide as a function of the range and of the source and the receiver directions are thus examined as an implementation example.