A treatment wetland (TW) is a nature-based solution for wastewater treatment that utilizes engineered systems to enhance processes that naturally occur in the environment. Computational Fluid Dynamics (CFD) is a consolidated technique used to design and perform the optimization of wastewater treatment systems. However, few CFD studies have addressed TW as this requires additional assumptions to properly include the porous flow modeling. This study aims to investigate the hydraulic patterns in horizontal subsurface flow treatment wetlands by using CFD coupled with the discrete element method (DEM), which provides a detailed representation of flow movement and porous media. The simulations were carried out in the software CFD-DEM. In terms of model settings, laminar flow and unsteady simulation were assumed. A case study using a 2 m long pilot-scale treatment wetland was used, as this model was previously validated for this setup. In order to support the design of a system with higher flow rates, three different flow rates were applied to check the differences in the system from the point of view of hydraulic behavior, in terms of flow direction, preferential flows, and hydraulic retention time. The results of the three scenarios were compared and a significant change was noted in the hydraulic behavior from the lower to the higher flow rate applied. The CFD-DEM simulations were effective to model the hydraulic patterns of the case under study and shown to be a good approach to accurately investigate hydraulics in TW, enabling operation optimization.