We derive a model to describe the interaction of an rf-SQUID (radio f requency Superconducting QUantum Interference Device) based metasurface with free space electromagnetic waves. The electromagnetic fields are described on the base of Maxwell's equations. For the rf-SQUID metasurface we rely on an equivalent circuit model. After a detailed derivation, we show that the problem that is described by a system of coupled differential equations is wellposed and, therefore, has a unique solution. In the small amplitude limit, we provide analytical expressions for reflection, transmission, and absorption depending on the frequency. To investigate the nonlinear regime, we numerically solve the system of coupled differential equations using a finite element scheme with transparent boundary conditions and the Crank-Nicolson method. We also provide a rigorous error analysis that shows convergence of the scheme at the expected rates. The simulation results for the adiabatic increase of either the field's amplitude or its frequency show that the metasurface's response in the nonlinear interaction regime exhibits bistable behavior both in transmission and reflection.