  We introduce a new adaptive and fully Bayesian grid-based method to model
strong gravitational lenses with extended images. The primary goal of this
method is to quantify the level of luminous and dark-mass substructure in
massive galaxies, through their effect on highly-magnified arcs and Einstein
rings. The method is adaptive on the source plane, where a Delaunay
tessellation is defined according to the lens mapping of a regular grid onto
the source plane. The Bayesian penalty function allows us to recover the best
non-linear potential-model parameters and/or a grid-based potential correction
and to objectively quantify the level of regularization for both the source and
the potential. In addition, we implement a Nested-Sampling technique to
quantify the errors on all non-linear mass model parameters -- ... -- and allow
an objective ranking of different potential models in terms of the marginalized
evidence. In particular, we are interested in comparing very smooth lens mass
models with ones that contain mass-substructures. The algorithm has been tested
on a range of simulated data sets, created from a model of a realistic lens
system. One of the lens systems is characterized by a smooth potential with a
power-law density profile, twelve include a NFW dark-matter substructure of
different masses and at different positions and one contains two NFW dark
substructures with the same mass but with different positions. Reconstruction
of the source and of the lens potential for all of these systems shows the
method is able, in a realistic scenario, to identify perturbations with masses
>=10^7 solar mass when located on the Einstein ring. For positions both inside
and outside of the ring, masses of at least 10^9 solar mass are required (...).
