To understand the mechanism of the turbulence transport in porous media, large eddy simulations have been performed by a multiple-relaxation time lattice Boltzmann method. The porous media presently considered are square rod arrays. Changing the square rod size, five different cases of the porous media are analysed. Using the simulation data, the budget terms of the volume averaged Reynolds stress equations are investigated in detail. It is found that the behaviour of the production, redistribution and dissipation terms drastically changes when the porosity of the porous media becomes higher than 0.75. When the porosity becomes higher, it is observed that a perturbation by the Kármán vortex shedding propagates in the whole domain. This perturbation is the source of the the cross-streamwise component of the Reynolds stress. It is confirmed that there is a relation between the porous characteristic scales and the turbulent length scale and one of the production terms of the transport equation of the volume averaged Reynolds stress can be modelled by the Darcy-Forchheimer form.