Using a Cauchy integral formulation of the boundary integral equations, we simulate the erosion a porous medium comprised of up to 100 solid bodies embedded in a Stokes flow. The grains of the medium are resolved individually and erode under the action of surface shear stress. Through nonlinear feedback with the surrounding flow fields, microscopic changes in grain morphology give way to larger-scale features in the medium such as channelization. The Cauchy-integral formulation and associated quadrature formulas enable us to resolve dense configurations of nearly contacting bodies. We observe substantial anisotropy to develop over the course of erosion; that is, the configurations that result from erosion generally permit flow in the longitudinal direction more easily than in the transverse direction by up to a factor of six. These results suggest that the erosion of solid material from groundwater flows may contribute to previously observed anisotropy of natural porous media.