Radar sounders (RS) provide information on subsurface targets for planetary investigations. Several simulation techniques have been developed to support the RS design and the data interpretation. Each technique has different properties and modeling capabilities, achieving different trade-offs between accuracy and computational requirements. The state-of-the-art RS simulation techniques include: i) numerical methods, such as the Finite- Difference Time-Domain (FDTD) technique, which allows the modelling of small-scale scattering phenomena at the cost of high computational requirements; ii) facet modeling and ray-tracing based methods, such as the multi-layered coherent RS (MCS) technique, which requires less computational resources than FDTD, allowing the modeling of large-scale scattering phenomena. Recently an integrated simulation methodology has been presented, for simulating small-scale scattering phenomena in large scenarios. However, this methodology was designed for modeling only surface scattering. In this paper, we propose a method that extends the capabilities of the integrated methodology to model both large and small-scale roughness in a multi-layer scenario. The proposed method uses the FDTD technique to evaluate the effects associated with small-scale roughness in terms of i) scattering phenomena associated with the layers and ii) power losses associated with the signal transmitted through a rough layer. To recursively apply scattering and transmission to multiple layers of the subsurface, a coherent ray-tracing method is used. We experimentally assessed the effectiveness of the proposed methodology on three-layer models by integrating the effect of roughness imposed on the layers and in transmission through them.