A numerical modeling effort of flow and Biochemical Oxygen Demand (BOD) fate in horizontal subsurface flow (HSF) constructed wetlands (CWs), under Mediterranean conditions, is presented. The Visual MODFLOW family computer code, based on the finite difference method, was used for the numerical simulation of flow in five pilot-scale HSF CW units containing various vegetation and porous material types, and operating under various temperature and hydraulic residence time (HRT) conditions. BUD fate and transport was simulated using the MT3DMS computer code. Experimental data from these CW units were used in estimating the BUD removal coefficient, required in the above modeling effort. Values of this removal coefficient were found to depend on temperature and HRT. Relations are presented for the prediction of this coefficient. The model was calibrated and verified using the measured BUD concentrations, and was used in test runs comparing performance under various vegetation, porous media size, temperature and HRT conditions. Model performance was found acceptable, indicating its usefulness in the simulation of BUD fate in CWs and in the design of these facilities. (C) 2012 Elsevier B.V. All rights reserved.