The parabolic equation (PE) is one of the most popular methods for the modeling of low-frequency sound in the atmosphere. Over the past decades, considerable efforts have been dedicated to the inclusion of medium inhomogeneities, such as wind and turbulence, extending the PE to more realistic atmospheric conditions. On the other hand, few models take topography into account, while its effects on infrasound propagation are important in some cases. A new development has enabled the inclusion of irregular terrain in the narrow-angle 3DPE for a moving inhomogeneous medium [Khodr et al., JASA (2020)]. A coordinate transform similar to the Beilis-Tappert mapping used in 2D [Parakkal et al., JASA (2012)] is introduced to account for the variation in topography. The numerical solution relies on a Crank–Nicolson implicit scheme along the propagating direction. The pressure field at every step is computed with an iterative fixed-point algorithm which consists of a succession of tridiagonal systems that can be efficiently solved. A validation is performed against 3D Boundary Element simulations for propagation above a Gaussian hill in a homogeneous atmosphere. The existence of transversal scattering in the shadow zone, not represented by 2D models, is highlighted. Further developments include the extension of the method to a split-step wide-angle formulation. [This work was conducted at the University of Bristol in partnership with AWE Blacknest.]