In what follows, all manifolds are smooth, Hausdorff, paracompact, connected and oriented, and all maps between any two of them are assumed to be smooth. Let $\pi:E\rightarrow M$ be a fiber bundle over the base manifold $M$, which we assume to have dimension $n$. Given a bundle map $$\mathscr{L}:J^kE\rightarrow\wedge^nT^*M$$ over $M$, where $J^kE$ is (the total space of) the jet bundle of order $k$ of $\pi$ - i.e. $\mathscr{L}$ is a Lagrangian density of order $k$ - we define the action functional $S_{\mathscr{L}}:C^\infty_c(M)\times\Gamma(\pi)\rightarrow\mathbb{R}$ associated to $L$ (here $C^\infty_c(M)$ is the space of smooth real-valued functions on $M$ with compact support, and $\Gamma(\pi)=\{\phi\in C^\infty(M,E)\ |\ \pi\circ\phi=\text{id}_M\}$ is the space of smooth sections of $\pi$) as $$S_{\mathscr{L}}(f,\phi)=\int_M f(j^k\phi)^*\mathscr{L}\ .$$ The Euler-Lagrange operator $E(\mathscr{L})$ associated to $\mathscr{L}$ is the (usually nonlinear) partial differential operator of order $\leq 2k$ defined by the formula $$\int_M E(\mathscr{L})[\phi](\vec{\phi})=\int_M\left.\frac{\partial}{\partial t}\right|_{t=0}f(j^k\phi_t)^*\mathscr{L}\ ,\quad\vec{\phi}\in\Gamma_c(\phi^*VE\rightarrow M)\ ,f|_{\text{supp}\vec{\phi}}\equiv 1\ ,$$ where $\phi_t=\Phi(t,\cdot)$ for any $\Phi\in C^\infty(I\times M,E)$ with $0\in I\subset\mathbb{R}$ an open interval, $\phi_0=\phi$, $\left.\frac{\partial\Phi}{\partial t}\right|_{t=0}=\vec{\phi}$ and $\phi_t\in\Gamma(\pi)$, $\phi_t|_{M\smallsetminus\text{supp}\vec{\phi}}=\phi|_{M\smallsetminus\text{supp}\vec{\phi}}$ for all $t\in I$. The above definition of $E(\mathscr{L})$ can be shown to be independent of the particular choice of $f$ and $\Phi$ under the above conditions. There are other possible ways to define $E(\mathscr{L})$ in the literature, all leading to the same object.
One can see from the above definition that $E(\mathscr{L})[\phi]\in\Gamma(\phi^*V^\circledast E\rightarrow M)$, where $\pi_{V^\circledast E}:V^\circledast E=\pi^*(\wedge^nT^*M)\otimes V^*E\rightarrow E$ is the so-called twisted dual of the vertical bundle $VE=\ker T\pi\rightarrow E$ of $E$. Since $$\phi^*V^\circledast E=\{(p,q)\in M\times V^\circledast E\ |\ \phi(p)=\pi_{V^\circledast E}(q)\}\subset M\times_M V^\circledast E$$ for all $\phi\in\Gamma(\pi)$, it is clear that $\phi^*V^\circledast E$ is a sub-bundle of $M\times_M V^\circledast E$ over $M$. As such, we can write $$E(\mathscr{L})=\rho(E(\mathscr{L}))\circ j^{2k}\ ,$$ where $$\rho(E(\mathscr{L})):J^{2k}(E)\rightarrow M\times_M V^\circledast E$$ is a bundle map covering $\text{id}_M$ ($M$ is seen above as a fiber bundle over itself with singleton fibers and $\text{id}_M$ as the projection map).
When one writes a local formula for $E(\mathscr{L})[\phi]$ using a local trivialization of $E$ and a local chart for the typical fiber $Q$ of $\pi$, one sees that $E(\mathscr{L})[\phi]$ is always an affine function of the highest-order derivatives of $\phi$ if the order of $E(\mathscr{L})$ happens to be nonzero (which we assume to be the case), so it would make sense to say that $E(\mathscr{L})$ is a "quasi-linear" partial differential operator. However, $M\times_M V^\circledast E$ is generally not a vector or affine bundle over $M$ (unless $E$ itself is), despite $\phi^*V^\circledast E$ being a (different) vector bundle over $M$ for each $\phi\in\Gamma(\pi)$, so $E(\mathscr{L})$ does not fit into the usual way to globally define quasi-linear partial differential operators (see e.g. the discussion in Section IX.2, pp. 393ff of the book Exterior Differential Systems by R.L. Bryant et al.), which requires the target bundle of the operator to be a vector (or at least an affine) bundle.
Question: in which global (i.e. coordinate-independent) way one may define quasi-linear partial differential operators in order to encompass Euler-Lagrange operators $E(\mathscr{L})$ acting on smooth sections of general fiber bundles?