"Let $P : \CI(M; E) \to \CI(M; F)$ be an order $\mu$ differential operator with coefficients $a$ and $P_k := P : H^{s_0 + k +\mu}(M; E) \to H^{s_0 + k}(M; F)$. We prove polynomial norm estimates for the solution $P_0^{-1}f$ of the form $$\|P_0^{-1}f\|_{H^{s_0 + k + \mu}(M; E)} \le C \sum_{q=0}^{k} \, \| P_0^{-1} \|^{q+1} \,\|a \|_{W^{|s_0|+k}}^{q} \, \| f \|_{H^{s_0 + k - q}},$$ (thus in higher order Sobolev spaces, which amounts also to a parametric regularity result). The assumptions are that $E, F \to M$ are Hermitian vector bundles and that $M$ is a complete manifold satisfying the Fr\'echet Finiteness Condition (FFC), which was introduced in (Kohr and Nistor, Annals of Global Analysis and Geometry, 2022). These estimates are useful for uncertainty quantification, since the coefficient $a$ can be regarded as a vector valued random variable. We use these results to prove integrability of the norm $\|P_k^{-1}f\|$ of the solution of $P_k u = f$ with respect to suitable Gaussian measures."