Abstract: We analyze theoretically and experimentally vortex configurations in mesoscopic superconducting squares. Our theoretical approach is based on the analytical solution of the London equation using Green's-function method. The potential-energy landscape found for each vortex configuration is then used in Langevin-type molecular-dynamics simulations to obtain stable vortex configurations. Metastable states and transitions between them and the ground state are analyzed. We present our results of the first direct visualization of vortex patterns in micrometer-sized Nb squares, using the Bitter decoration technique. 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 formation of vortex “shells”.