In the present work, a k-epsilon model, based on the work of Lee and Howell (Proceedings of the ASME-JSME Thermal Engineering Hawaii, 1987), is rigorously derived based on time average of spatially averaged Navier-Stokes equations. The model is then employed to solve for a flow in a backward-facing step channel with a porous insert. The numerical solver is modified from the STREAM code (Lien and Leschziner, Comput. Meth. Appl. Mech. Eng. 114 (1994a) 123-148), and it has been validated against the experimental data of Seegmiller and Driver (AIAA Journal 23 (1985) 163-171). The code is then used to perform simulation for cases with a porous insert. The resistance of the porous insert can be altered by changing its permeability (k), Forchheimer's constant (F), or thickness (b). The goal is to examine the influence of each parameter on the resulting flow and turbulent kinetic energy (k) distributions. It is discovered that, by increasing the resistance of the insert, flow eventually enters a transitional regime towards relaminarization. This is due to the contribution of Darcy's and Forchheimer's terms in the governing equations, and modifying these two terms changes the levels of P-k and, hence, k and . Generally speaking, lowering k or raising F results in a greater suppression of P-k than epsilon, causing the flow to relaminarize. Meanwhile, if the pore size is reasonably large to sustain turbulence within the porous media, increasing b reduces but does not eliminate the turbulent activity in the porous insert.