  We develop numerical methods to simulate the fluid-mechanical erosion of many
bodies in two-dimensional Stokes flow. The broad aim is to simulate the erosion
of a porous medium (e.g. groundwater flow) with grain-scale resolution. Our
fluid solver is based on a second-kind boundary integral formulation of the
Stokes equations that is discretized with a spectrally-accurate Nystrom method
and solved with fast-multipole-accelerated GMRES. The fluid solver provides the
surface shear stress which is used to advance solid boundaries. We regularize
interface evolution via curvature penalization using the $\theta$-$L$
formulation, which affords numerically stable treatment of stiff terms and
therefore permits large time steps. The overall accuracy of our method is
spectral in space and second-order in time. The method is computationally
efficient, with the fluid solver requiring $\mathcal{O}(N)$ operations per
GMRES iteration, a mesh-independent number of GMRES iterations, and a one-time
$\mathcal{O}(N^2)$ computation to compute the shear stress. We benchmark
single-body results against analytical predictions for the limiting morphology
and vanishing rate. Multibody simulations reveal the spontaneous formation of
channels between bodies of close initial proximity. The channelization is
associated with a dramatic reduction in the resistance of the porous medium,
much more than would be expected from the reduction in grain size alone.
