Abstract: We study vortex configurations in mesoscopic superconducting squares. Our theoretical approach is based on the analytical solution of the London equation using the Green's function method. The potential energy landscape found is then used in Langevin-type molecular-dynamics simulations to obtain stable vortex configurations. We show that the filling rules for vortices in squares with increasing applied magnetic field can be formulated, although in a different manner than in disks, in terms of the formation of vortex 'shells'. We discuss metastable states and the stability of the vortex configurations found with respect to variations of the material parameters and deformations of the shape of the sample.