  The few determinant (FED) methodology, introduced in our previous works for
1D lattices, is here adapted for the repulsive two-dimensional Hubbard model at
half-filling and with finite doping fractions. Within this configuration mixing
scheme, a given ground state with well defined spin and space group quantum
numbers, is expanded in terms of a nonorthogonal symmetry-projected basis
determined through chains of variation-after projection calculations. The
results obtained for the ground state and correlation energies of half-filled
and doped 4$\times$4, 6$\times$6, 8$\times$8, and 10$\times$10 lattices, as
well as momentum distributions and spin-spin correlation functions in small
lattices, compare well with those obtained using other state-of-the-art
approximations. The structure of the intrinsic determinants resulting from the
variational strategy is interpreted in terms of defects that encode information
on the basic units of quantum fluctuations in the considered 2D systems. The
varying nature of the underlying quantum fluctuations, reflected in a
transition to a stripe regime for increasing onsite repulsions, is discussed
using the intrinsic determinants belonging to a 16$\times$4 lattice with 56
electrons. Such a transition is further illustrated by computing spin-spin and
charge-charge correlation functions with the corresponding multireference FED
wave functions. In good agreement with previous studies, the analysis of the
pairing correlation functions reveals a weak enhancement of the extended
$s$-wave and d$_{x^2-y^2}$ pairing modes. Given the quality of results here
reported together with those previously obtained for 1D lattices and the
parallelization properties of the FED scheme, we believe that symmetry
projection techniques are very well suited for building ground state wave
functions of correlated electronic systems, regardless of their dimensionality.
