$\newcommand{\id}{\operatorname{Id}}$ $\newcommand{\TM}{\operatorname{TM}}$ $\newcommand{\Hom}{\operatorname{Hom}}$ $\newcommand{\M}{\mathcal{M}}$ $\newcommand{\N}{\mathcal{N}}$ $\newcommand{\tr}{\operatorname{tr}}$ $\newcommand{\TM}{\operatorname{T\M}}$ $\newcommand{\TN}{\operatorname{T\N}}$ $\newcommand{\TstarM}{T^*\M}$
Edit: I narrowed the focus of the question.
Summary: I suspect $\delta(df \wedge df)=0$ is not an E-L equation because it pus too many constraints on $f$. Can this heuristic be formalized?
Let $\M,\N$ be $d$-dimensional oriented Riemannian manifolds, and let $f:\M \to \N$ be smooth; $df \in \Omega^1(\M,f^*{\TN})$, and for $1 < k \le d$ let $\bigwedge^k df\in \Omega^k\big(\M,\Lambda_k(f^*{\TN})\big)$ be the induced map.
Notation:
$\nabla^{\Lambda_k(f^*{\TN})}$ is the induced connection on $\Lambda_k(f^*{\TN})$, $\delta_{\nabla^{\Lambda_k(f^*{\TN})}} $ is the adjoint of $d_{\nabla^{\Lambda_k(f^*{\TN})}}$.
Question: Is $\delta_{\nabla^{\Lambda_k(f^*{\TN})}} \big( \bigwedge^k df \big) =0$ an Euler-Lagrange equation of some functional $E(f)$?
I suspect the answer is negative.
Edit: Can this be proved using known results on inverse problems in calculus of variations? e.g if we can show the equation in the Euclidean setting ($M=N=\mathbb{R}^d$) is not an E-L eq, then we are done.
(I am not aware of much work on such inverse problems in general Riemannian settings, but if the Euclidean can be decided, it's enough).
Heuristic:
$\delta_{\nabla^{\Lambda_k(f^*{\TN})}} \big( \bigwedge^k df \big) =0$ puts too many constraints on $f$:
$\delta_{\nabla^{\Lambda_k(f^*{\TN})}} \big( \bigwedge^k df \big) \in \Omega^{k-1}(\M,\Lambda_k(f^*{\TN}))$, so locally the "equation" is a system of ${d \choose k-1}{d \choose k}$ scalar equations.
However, it seems to me an Euler-Lagrange eq of a functional cannot consist of more than $\dim(f^*\TN)=d$ equations, since (roughly) this is the number of degrees of freedom we have in choosing the variation field $V\in\Gamma(f^*\TN)$.
More explicitly, since $dE(V)$ is linear in $V$, it should always be in the form of $dE(V)=\langle A(f),V \rangle$, where $A(f) \in \Gamma(f^*\TN)$ so the E-L eq should be in the form of $A(f)=0$.
Of course, this argument gives only an upper bound to the number of "independent" constraints - degeneracies can occur (e.g Null-Lagrangians).
So, suppose ${d \choose k-1}{d \choose k}>d$. Is it possible $\delta_{\nabla^{\Lambda_k(f^*{\TN})}} \big( \bigwedge^k df \big) =0$ is degenerate and reduces to no more than $d$ independent eqs (locally)?
Can we prove this is not happening? (If we can, than this shows our equation is not an E-L eq of any functional).
Analysis of the borderline cases:
We now turn to see what happens when ${d \choose k-1}{d \choose k}=d$. This happens (for $1 \le k \le d$) iff $k=1$ or $k=d$. The former corresponds to harmonicity. Suppose that $k=d$. We will try to understand if our equation is an E-L equation in this case.
For concreteness, let's work with $k=d=2$.
I proved in my answer that the E-L equation of the functional $$ E_2(f)=\frac{1}{2}\int_{M} \| \bigwedge^2 df\|^2 \text{Vol}_{M}, \, \text{is}$$ $$h_{f^*\TN} \big( \tr_{\TM}\big( df \otimes \delta_{\nabla^{\Lambda_2(f^*{\TN})}}(\bigwedge^2 df)\big) \bigg)=0.$$ Since $h_{f^*\TN}$ is injective in this case, the E-L eq is equivalent to $$\tr_{\TM}\big( df \otimes \delta_{\nabla^{\Lambda_2(f^*{\TN})}}(\bigwedge^2 df)\big)=0. \tag{1}$$
When restricting the discussion to immersions, this reduces to $ \delta_{\nabla^{\Lambda_2(f^*{\TN})}}(\bigwedge^2 df)=0$.
It is not clear what happens in the general case; If we don't assume $f$ is an immersion, are the equations equivalent? (If not, perhaps there is another way to realize $\delta_{\nabla^{\Lambda_2(f^*{\TN})}}(\bigwedge^2 df)=0$ as an E-L eq).
The problem is that we cannot conclude (at least not immediately) that at every point where $df$ is not invertible, $\delta_{\nabla^{\Lambda_2(f^*{\TN})}}(\bigwedge^2 df)=0$. (We know $\bigwedge^2 df=0$ at the point, but $\delta$ is a differential operator, it sees beyond the pointwise behaviour).