This study aimed at optimizing horizontal subsurface flow constructed wetlands (CWs) to improve hydraulic performance and pollutant removal efficiency. A groundwater modeling package (MODFLOW) was used to optimize three design parameters (length-to-width ratio, inlet/outlet-to-length ratio, and substrate size configuration). Using the optimized parameters, three pilot-scale CWs were built to treat actual wastewater. For model validation, we used a tracer test to evaluate hydraulic performance, and investigated the pollutant spatial distributions and removal efficiencies. We conclude that MODFLOW is suitable for designing CWs, accurately predicting that increasing hydraulic conductivity from surface to bottom layers could improve performance. However, the effect of vegetation, which decreased the hydraulic conductivity of the surface layer, should be considered to improve simulation results. Multilayer substrate configuration, with increasing hydraulic conductivity from the surface to bottom layers, significantly increased pollutant removal compared with monolayer configuration. The spatial variation in pollutant transport and degradation through the filling substrate showed that the multilayer configuration was able to increase use of the available space and moderately reduced short-circuiting and dead zones. Thus, multilayer CWs had higher experimental retention times, effective volume fractions and hydraulic efficiencies, and lower short-circuiting compared with monolayer CWs operating under similar conditions.