It turns out that the statement in question appeared in full beauty in
There in fact is given the more general statement that the full Eilenberg-Moore category of the jet comonad is equivalent to that of partial differential equations in variables in the given base space.
Of course the Kleisli category in question is a full subcategory of that, and the statement that I asked for is the special case of the above made fully explicitly in
Now since this statement about differential operators is true, of course every other true statement about differential operators that one finds elsewhere in the literature will be compatible with it, but Marvan is the only author I am aware of who explicitly states the neat general abstract statement that I asked for. It's one of those abstract trivialities, if you wish, that are worth making explicit.
On the other hand, Marvan checks that $\mathrm{Jet}$ is a comonad directly, without realizing it as the base change comonad $\mathrm{Jet}_X \simeq i^\ast_X (i_X)_\ast$ along the unit of the de Rham stack monad $i_X \colon X \to X_{dR}$. For that statement of course one needs to be in a model of synthetic differential geometry in order for $X_{dR}$ to exist, such as (but not necessarily) algebraic geometry.
In the algebraic context, the reference stating explicitly $\mathrm{Jet} \simeq i^\ast_X (i_X)_\ast $ that I am aware of is
Here it is maybe noteworthy that from this it follows immediately that the Jet-comonad is right adjoint to the "infinitesimal disk bundle" monad $T_{\mathrm{inf}} := i_X^\ast (i_X)_!$.
This adjunction $(T_{\mathrm{inf}} \dashv \mathrm{Jet})$ itself, without however its origin in the base change adjoint triple of $i_X$, was observed in the context of synthetic differential geometry in
I'll just add that one immediate neat implication of making the general abstract comonad theory behind this explicit is that it gives in full generality that for any topos (or $\infty$-topos) $\mathbf{H}$ equipped with an "infinitesimal shape modality" $X\mapsto \Im X = X_{dR}$, then since $\mathrm{Jet} := (i_X)^\ast (i_X)_\ast \colon \mathbf{H}_{/X} \to \mathbf{H}_{/X}$ is a right adjoint, a standard fact (here for toposes, here maybe for $\infty$-toposes) gives that its EM-category, hence the category $\mathrm{PDE}(X)$ of PDEs in $X$ is itself a topos, sitting by a geometric morphism
$$\mathrm{PDE}(X) \stackrel{{\longleftarrow}}{\longrightarrow} \mathbf{H}_{/X}$$
over $\mathbf{H}$, whose direct image is (non-fully) that co-Kleisli category of bundles over $X$ with differential operators between them.
Now moreover, using the infinitesimal shape modality $\Im$ it also follows that on $\mathbf{H}_{/X}$ there is an étalification coprojection $Et$ (here). This is such that when $E \in \mathbf{H}$ is a coefficient for a differential cohomology theory (e.g. any stable object in the case that $\mathbf{H}$ is cohesive (here)), then $Et(X^\ast E)$ has the interpretation of being the "bundle of $E$-connections" over $X$. So with the above geometric morphism we may transfer all this to the topos $\mathrm{PDE}(X)$ along the composite
$$
\iota \colon \mathbf{H} \stackrel{X^\ast}{\longrightarrow} \mathbf{H}_{/X}\stackrel{Et}{\longrightarrow} \mathbf{H}_{/X} \stackrel{}{\longrightarrow} \mathrm{PDE}(X)
\,.
$$
Unwinding what this all means, in view of Marvan's insight above, it follows, I think, that for any bundle $F \in \mathbf{H}_{/X} \to \mathrm{PDE}(X)$ then maps in $\mathrm{PDE}(X)$ of the form
$$
F \longrightarrow \iota E
$$
are equivalently horizontal "$E$-differential forms" on the jet bundle of $F$, in the sense of variational bicomplex theory.
That might count as some non-trivial mileage gained out of making explicit the statement that I was asking for.