  We consider associated $Zh$ production with $Z \to \ell^+ \ell^-$ and $h \to
b \bar b$ decays in hadronic collisions. In the framework of the Standard Model
effective field theory (SMEFT) we calculate the QCD corrections to this process
and achieve next-to-next-to-leading order plus parton shower (NNLO$+$PS)
accuracy using the MiNNLO$_{\rm PS}$ method. This precision is obtained for a
subset of six SMEFT operators, including the corrections from effective Yukawa-
and chromomagnetic dipole-type interactions. Missing higher-order QCD effects
associated with the considered dimension-six operators are estimated to have a
relative numerical impact of less than a percent on the total rate once
existing experimental limits on the relevant Wilson coefficients are taken into
account. We provide a dedicated Monte Carlo (MC) code that evaluates the NNLO
SMEFT corrections on-the-fly in the event generation. This MC generator is used
to study the numerical impact of NNLO$+$PS corrections on the kinematic
distributions in $pp \to Zh \to \ell^+ \ell^- b \bar b$ production employing
simple SMEFT benchmark scenarios. We identify the invariant mass $m_{b \bar b}$
of the two $b$-tagged jets as well as the three-invariant jet mass $m_{b \bar b
j}$ as particularly interesting observables to study SMEFT effects. These
distributions receive contributions that change both their normalisation and
shape with the latter modifications depending on the exact jet definition. To
our knowledge SMEFT effects of this type have so far not been discussed in the
literature. The presented MC generator can also serve as a starting point to
obtain NNLO$+$PS accuracy for a suitable enlarged set of effective operators in
the future.
