Variational calculations of ground-state properties of $^4$He, $^{16}$O, and $^{40}$Ca are carried out employing realistic phenomenological two- and three-nucleon potentials. The trial wave function includes two- and three-body correlations acting on a product of single-particle determinants. Expectation values are evaluated with a cluster expansion for the spin-isospin dependent correlations considering up to five-body cluster terms. The optimal wave function is obtained by minimizing the energy expectation value over a set of up to 20 parameters by means of a nonlinear optimization library. We present results for the binding energy, charge radius, one- and two-body densities, single-nucleon momentum distribution, charge form factor, and Coulomb sum rule. We find that the employed three-nucleon interaction becomes repulsive for $A\geq16$. In $^{16}$O the inclusion of such a force provides a better description of the properties of the nucleus. In $^{40}$Ca instead, the repulsive behavior of the three-body interaction fails to reproduce experimental data for the charge radius and the charge form factor. We find that the high-momentum region of the momentum distributions, determined by the short-range terms of nuclear correlations, exhibit a universal behavior independent of the particular nucleus. The comparison of the Coulomb sum rules for $^4$He, $^{16}$O, and $^{40}$Ca reported in this work will help elucidate in-medium modifications of the nucleon form factors.