We study numerically two-dimensional creeping viscoelastic flow past a biperiodic square array of cylinders within the Oldroyd B, fene-CR and fene-P constitutive models of dilute polymer solutions. Our results capture the initial mild decrease then dramatic upturn ('thickening') seen experimentally in the drag coefficient as a function of increasing Weissenberg number. By systematically varying the porosity of the flow geometry, we demonstrate two qualitatively different mechanisms underpinning this thickening effect: one that operates in the highly porous case of widely spaced obstacles, and another for more densely packed obstacles, with a crossover between these two mechanisms at intermediate porosities. We also briefly consider 2D creeping viscoelastic flow past a linear array of cylinders confined to a channel. Our results here suggest that recent numerical reports of viscoelastic turbulence in this 2D geometry may need to be checked for convergence with respect to grid resolution, with a possible return to steady flow for finer grids. ; Comment: 16 pages, 16 figures