The straightforward way to do this is to note that we have, for some $u_1$ and $u_2$, $$ \mathrm{d}u = u_1\,\eta^1 + u_2\,\eta^2 $$ Taking the exterior derivative of this and using Cartan's Lemma then shows that, for some $u_{11}$, $u_{12}$, and $u_{22}$, $$ \begin{aligned} \mathrm{d}u_1 &= \phantom{-}u_2\,\omega + u_{11}\,\eta^1 + u_{12}\,\eta^2\\ \mathrm{d}u_2 &= -u_1\,\omega + u_{12}\,\eta^1 + u_{22}\,\eta^2 \end{aligned} $$ Then taking exterior derivatives of these equations and using Cartan's Lemma yields $$ \begin{aligned} \mathrm{d}u_{11} &= & 2u_{12}\,&\omega + u_{111}\,\eta^1 + u_{112}\,\eta^2\\ \mathrm{d}u_{12} &= &(u_{22}-u_{11})\,&\omega + u_{121}\,\eta^1 + u_{122}\,\eta^2\\ \mathrm{d}u_{22} &= &-2u_{12}\,&\omega + u_{221}\,\eta^1 + u_{222}\,\eta^2\\ \end{aligned} $$ for some $u_{ijk}$ satisfying $u_{112}-u_{121} = -K\,u_2$ and $u_{122}-u_{221} = K\,u_1$.
Now just plug these formulae into the left hand side of Bochner's formula and compute.
Comment: I understand that the OP would like to have some kind of exterior algebra 'magic' that just makes the calculation a consequence of $\mathrm{d}^2=0$, but that's an unreasonable request, it seems to me. The strength of the method of calculating by moving coframes is that it does show that all the local identities that you can prove are consequences of $\mathrm{d}^2=0$ and exterior algebra, though you may have to apply this principle multiple times if the identity you want involves higher derivatives.
For example, in this particular case, the general structure equations for an orthonormal coframe (i.e., $g = {\eta_1}^2 + \cdots + {\eta_n}^2$) are that there are unique $1$-forms $\eta_{ij}=-\eta_{ji}$ satisfying $$ \begin{aligned} \mathrm{d}\eta_i &= -\eta_{ij}\wedge\eta_j\\ \mathrm{d}\eta_{ij} &= -\eta_{ik}\wedge\eta_{kj} + \tfrac12 R_{ijkl}\,\eta_k\wedge\eta_l \end{aligned}\tag1 $$ where $R_{ijkl}=-R_{ijlk}$ and all terms are summed on repeated indices.
Now, if $u$ is a function in the domain of $\eta$, we have $$ \mathrm{d}u = u_i\,\eta_i\tag2 $$ for some functions $u_i$. Applying the exterior derivative and then Cartan's Lemma to (2) yields $$ \mathrm{d}u_i = -u_j\,\eta_{ij} + u_{ij}\,\eta_j\tag3 $$ for some functions $u_{ij}=u_{ji}$. Applying the exterior derivative and then Cartan's Lemma to (3) yields $$ \mathrm{d}u_{ij} = -u_{ik}\,\eta_{kj} -u_{jk}\,\eta_{ki}+ u_{ijk}\,\eta_k\tag4 $$ where $u_{ilk}-u_{ikl} = R_{ijkl}\,u_j$ (and, of course, $u_{ijk}=u_{jik}$).
Now note that $|\nabla u|^2 = u_iu_i$ and then plug in the above formulae for derivatives to get $$ \tfrac12\Delta|\nabla u|^2 = u_{ij}u_{ij} + u_iu_{ijj}. $$ Since, by the above formulae, $u_{ijj}=u_{jij} = u_{jji} + R_{jlji}\,u_l$, this yields $$ \tfrac12\Delta|\nabla u|^2 = u_iu_{jji} + u_{ij}u_{ij} + R_{jlji}\,u_iu_l\,, $$ and the rest is just interpreting the three terms, where, in general dimension, the final term involving curvature becomes $$ \mathrm{Ric}_g(\nabla u,\nabla u). $$
Whether you choose to use moving coframes or not is a matter of taste, of course, since you can equally well work with the definition of curvature in terms of the Lie brackets of the dual vector fields. (Of course, the $u_i$, $u_{ij}$ and $u_{ijk}$ are just the result of 'covariantly differentiating' $u$ with respect to these vector fields.)