Evidence is presented for the electroweak (EW) production of two jets (jj) in association with two Z bosons and constraints on anomalous quartic gauge couplings are set. The analysis is based on a data sample of proton-proton collisions at (\sqrt{s}=13\phantom{\rule{0.2em}{0ex}}\text{TeV}) collected with the CMS detector in 2016–2018, and corresponding to an integrated luminosity of 137(\phantom{\rule{0.2em}{0ex}}{\text{fb}}^{-1}). The search is performed in the fully leptonic final state (ZZ\to \ell \ell {\ell }^{\prime }{\ell }^{\prime }), where (\ell,{\ell }^{\prime }=e,\mu ). The EW production of two jets in association with two Z bosons is measured with an observed (expected) significance of 4.0 (3.5) standard deviations. The cross sections for the EW production are measured in three fiducial volumes and the result is ({\sigma }{\mathrm{EW}}\left(pp\to ZZ\mathrm{jj}\to \ell \ell {\ell }^{\prime }{\ell }^{\prime }\mathrm{jj}\right)={0.33}{-0.10}^{+0.11}\phantom{\rule{0.2em}{0ex}}{\text{(stat)}}_{-0.03}^{+0.04}\phantom{\rule{0.2em}{0ex}}\text{(syst)}\phantom{\rule{0.20em}{0ex}}\text{fb}) in the most inclusive volume, in agreement with the standard model prediction of (0.275±0.021\phantom{\rule{0.20em}{0ex}}\text{fb}). Measurements of total cross sections for jj production in association with two Z bosons are also reported. Limits on anomalous quartic gauge couplings are derived in terms of the effective field theory operators T0, T1, T2, T8, and T9.