A lattice Boltzmann bounce-back scheme that can simulate the first-order slip boundary condition at a fluid-solid interface is derived and Verified with analytical and asymptotic solutions for shear and pressure driven flows in channels, through square arrays of cylinders, and in a Couette cell. We then simulated Stokes flows with surface slip through simple cubic (SC), body-centered cubic (BCC), and face-centered cubic (FCC) arrays of spheres and established the slip (apparent) permeability of these arrays as a function of the no-slip (absolute) permeability, the slip length and the porosity. It is found that k(s)/k(ns), the ratio between the slip permeability and the no-slip permeability, increases with increasing slip length and increases with decreasing porosity owning to enhancement of flow through the gaps between particles. Interestingly, k(s)/k(ns), is nearly independent of sphere configuration when the dimensionless slip length is small; with increasing dimensionless slip length, however, different sphere configurations start to have different k(s)/k(ns). The dependence of k(s)/k(ns) on the slip length is slightly nonlinear owning to the curvature of the surface and this dependence can be correlated by quadratic equations. (C) 2017 Elsevier B.V. All rights reserved.