Abstract The Wall-Modeled-Large Eddy Simulation (WMLES) technique was used to analyze turbulent flow in a porous medium composed of a staggered arrangement of square cylinders. The study focused on a porosity range of 0.3 to 0.8 at a pore Reynolds number of 10×103. The research provides novel information on pressure and shear force distributions around the cylinders and time-averaged values of particle drag and fluctuating lift coefficients. Results indicate that the flow physics are mainly driven by an expansion region in front of the central particle and another in the rear, where the flow is contracted and strongly accelerated. In the first region, velocities, turbulence kinetic energy (k), and its dissipation rate (e) are attenuated by a sudden pressure increase, while in the contraction region, the effect is the opposite. As porosity decreases, flow gradients and the overall levels of k, e, and Reynolds stresses become more significant. Turbulent anisotropy increases as porosity decreases below 0.6. Hydrodynamic stresses in the last interval present relatively uniform levels, which is a unique feature that may be considered in designing dedicated engineering devices with controlled stresses.