Start with my 'wave equation in spacetime' worksheet for an introduction to a scalar wave propagating in curved spacetime: https://thenumbernine.github.io/lua/symmath/tests/output/wave%20equation%20in%20spacetime.html

unit conversion: $dx^0 = c dt, \partial_0 = \frac{1}{c} \partial_t$

Written as a linear system:

$\left[ \begin{matrix} \Pi \\ \Psi_i \end{matrix} \right]_{,0} = \left[ \begin{matrix} \Pi_{,k} \beta^k + \alpha \Psi_{j,k} \gamma^{jk} + \Psi^j \alpha_{,j} + K \Pi \alpha - \alpha \Psi_i \Gamma^i - f \alpha \\ \alpha \delta^k_i \Pi_{,k} + \beta^k \delta^j_i \Psi_{j,k} + \Psi_j {\beta^j}_{,i} + \alpha_{,i} \Pi \end{matrix} \right]$

$\left[ \begin{matrix} \Pi \\ \Psi_i \end{matrix} \right]_{,0} + \left[ \begin{matrix} - \Pi_{,k} \beta^k - \alpha \Psi_{j,k} \gamma^{jk} \\ - \alpha \delta^k_i \Pi_{,k} - \beta^k \delta^j_i \Psi_{j,k} \end{matrix} \right] = \left[ \begin{matrix} \Psi_i \alpha_{,j} \gamma^{ij} + K \Pi \alpha - \alpha \Psi_i \Gamma^i - f \alpha \\ \Psi_j {\beta^j}_{,i} + \alpha_{,i} \Pi \end{matrix} \right]$

$\left[ \begin{matrix} \Pi \\ \Psi_i \end{matrix} \right]_{,t} + c \left[ \begin{matrix} -\beta^k & -\alpha \gamma^{jk} \\ - \alpha \delta^k_i & - \beta^k \delta^j_i \end{matrix} \right] \left[ \begin{matrix} \Pi \\ \Psi_j \end{matrix} \right]_{,k} = \left[ \begin{matrix} \Psi_i \alpha_{,j} \gamma^{ij} + K \Pi \alpha - \alpha \Psi_i \Gamma^i - f \alpha \\ \Psi_j {\beta^j}_{,i} + \alpha_{,i} \Pi \end{matrix} \right]$

Let $n = n^i \partial_i$ be the vector normal to the flux surface.
Contravariant normals would have values $n^i = \gamma^{ij} n_j$.
The length of the normal is given by $|n| = \sqrt{n^i n_i}$.

Then normals along the cartesian $x,y,z$ would have values, i.e. in the k direction: $n_i = \delta_{ik}$ and $n^i = \gamma^{ik}$.
The length of a normal $|n| = \sqrt{n^i n_i}$ along any cartesian basis in the k direction would be $|n| = \sqrt{\underset{i}{\Sigma} (\gamma^{ik} \delta_{ik})} = \sqrt{\gamma^{kk}}$
Also let $n_2$ and $n_3$ be orthogonal to $n$ such that $n^i (n_2)_i = n^i (n_3)_i = 0$ and $|n_2| = |n_3| = |n|$

(TODO use a different letter or a denotation since $n^\mu$ is already the spatial hypersurface normal in ADM formalism. How about $(n_\Sigma)^i$, or $\mu^i$, or $m^i$ ... or sticking with the original ADM notation of ${}^4 g_{\mu\nu}$ being the 4-metric and ${}^4 {\Gamma^\alpha}_{\beta\gamma}$ being the 4-connection, how about ${}^4 n^i$ is the 4-normal and $n^i$ is the 3-normal?)

Flux Jacobian:

$F = \left[ \begin{matrix} -c \beta^k & -c \alpha \gamma^{jk} \\ -c \alpha \delta^k_i & -c \beta^k \delta^j_i \end{matrix} \right]$

Flux Jacobian in terms of flux surface normal:

$F = \left[ \begin{matrix} -c \beta_n & -c \alpha n^j \\ -c \alpha n_i & -c \beta_n \delta^j_i \end{matrix} \right]$

...where $\beta^n = \beta^i n_i$

Expand components into x,y,z:

$F = \left[ \begin{matrix} -c \beta_n & -c \alpha n^x & -c \alpha n^y & -c \alpha n^z \\ -c \alpha n_x & -c \beta_n & 0 & 0 \\ -c \alpha n_y & 0 & -c \beta_n & 0 \\ -c \alpha n_z & 0 & 0 & -c \beta_n \end{matrix} \right]$

---

Eigen decomposition:

Acoustic matrix:

$A = F + c \beta_n I$

$A = \left[ \begin{matrix} 0 & -c \alpha n^x & -c \alpha n^y & -c \alpha n^z \\ -c \alpha n_x & 0 & 0 & 0 \\ -c \alpha n_y & 0 & 0 & 0 \\ -c \alpha n_z & 0 & 0 & 0 \end{matrix} \right]$

Eigenvalues of A:
$A R_A = \lambda R_A$
$(A - I \lambda) R_A = 0$
$|A - I \lambda| = 0$

$\left| \begin{matrix} -\lambda & - c \alpha n^x & - c \alpha n^y & - c \alpha n^z \\ -c \alpha n_x & -\lambda & 0 & 0 \\ -c \alpha n_y & 0 & -\lambda & 0 \\ -c \alpha n_z & 0 & 0 & -\lambda \end{matrix} \right| = 0$

$c \alpha n_z \left| \begin{matrix} - c \alpha n^x & - c \alpha n^y & - c \alpha n^z \\ -\lambda & 0 & 0 \\ 0 & -\lambda & 0 \\ \end{matrix} \right| - \lambda \left| \begin{matrix} -\lambda & - c \alpha n^x & - c \alpha n^y \\ -c \alpha n_x & -\lambda & 0 \\ -c \alpha n_y & 0 & -\lambda \end{matrix} \right| = 0$

$-c^2 \alpha^2 n_z n^z \lambda^2 - \lambda \left( -c \alpha n_y \left| \begin{matrix} - c \alpha n^x & - c \alpha n^y \\ -\lambda & 0 \end{matrix} \right| - \lambda \left| \begin{matrix} -\lambda & - c \alpha n^x \\ -c \alpha n_x & -\lambda \end{matrix} \right| \right) = 0$

$-c^2 \alpha^2 n_z n^z \lambda^2 - \lambda \left( c^2 \alpha^2 n_y n^y \lambda - \lambda ( \lambda^2 - c^2 \alpha^2 n^x n_x ) \right) = 0$

$\lambda^2 ( \lambda^2 - c^2 \alpha^2 n_x n^x - c^2 \alpha^2 n_y n^y - c^2 \alpha^2 n_z n^z ) = 0$

$\lambda^2 (\lambda^2 - c^2 \alpha^2 |n|^2) = 0$
$\lambda^2 (\lambda + c \alpha |n|) (\lambda - c \alpha |n|) = 0$

...has eigenvalues:
$\lambda = \{-\alpha c |n|, 0, 0, \alpha c |n| \}$

For $\lambda = 0$:

$\left[ \begin{matrix} 0 & -c \alpha n^x & -c \alpha n^y & -c \alpha n^z \\ -c \alpha n_x & 0 & 0 & 0 \\ -c \alpha n_y & 0 & 0 & 0 \\ -c \alpha n_z & 0 & 0 & 0 \end{matrix} \right] \left[ \begin{matrix} X^\Pi \\ X^{\Psi_x} \\ X^{\Psi_y} \\ X^{\Psi_z} \end{matrix} \right] = 0$

$\left[ \begin{matrix} 0 & n^x & n^y & n^z \\ 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{matrix} \right] \left[ \begin{matrix} X^\Pi \\ X^{\Psi_x} \\ X^{\Psi_y} \\ X^{\Psi_z} \end{matrix} \right] = 0$

...we have the solution...

$\left[ \begin{matrix} X^\Pi \\ X^{\Psi_x} \\ X^{\Psi_y} \\ X^{\Psi_z} \end{matrix} \right] = \left[ \begin{matrix} 0 & 0 \\ (n_2)_x & (n_3)_x \\ (n_2)_y & (n_3)_y \\ (n_2)_z & (n_3)_z \end{matrix} \right]$

For $\lambda = \pm \alpha c |n|$:

$\left[ \begin{matrix} \mp c \alpha |n| & -c \alpha n^x & -c \alpha n^y & -c \alpha n^z \\ -c \alpha n_x & \mp c \alpha |n| & 0 & 0 \\ -c \alpha n_y & 0 & \mp c \alpha |n| & 0 \\ -c \alpha n_z & 0 & 0 & \mp c \alpha |n| \end{matrix} \right] \left[ \begin{matrix} X^\Pi \\ X^{\Psi_x} \\ X^{\Psi_y} \\ X^{\Psi_z} \end{matrix} \right] = 0$

...we have the solution...

$\left[ \begin{matrix} X^\Pi \\ X^{\Psi_x} \\ X^{\Psi_y} \\ X^{\Psi_z} \end{matrix} \right] = \left[ \begin{matrix} \mp |n| \\ n_x \\ n_y \\ n_z \end{matrix} \right]$

Combined:

$\Lambda_A = \left[ \begin{matrix} -c \alpha |n| & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & c \alpha |n| \end{matrix} \right]$

$R_A = \left[ \begin{matrix} |n| & 0 & 0 & -|n| \\ n_x & (n_2)_x & (n_3)_x & n_x \\ n_y & (n_2)_y & (n_3)_y & n_y \\ n_z & (n_2)_z & (n_3)_z & n_z \end{matrix} \right]$

$L_A = (R_A)^{-1} = \left[ \begin{matrix} \frac{1}{2|n|} & \frac{1}{2|n|^2} n^x & \frac{1}{2|n|^2} n^y & \frac{1}{2|n|^2} n^z \\ 0 & \frac{1}{|n|^2} (n_2)^x & \frac{1}{|n|^2} (n_2)^y & \frac{1}{|n|^2} (n_2)^z \\ 0 & \frac{1}{|n|^2} (n_3)^x & \frac{1}{|n|^2} (n_3)^y & \frac{1}{|n|^2} (n_3)^z \\ -\frac{1}{2|n|} & \frac{1}{2|n|^2} n^x & \frac{1}{2|n|^2} n^y & \frac{1}{2|n|^2} n^z \\ \end{matrix} \right]$

Now since we have no change of basis from F to A we can say that $R_F = R_A$ and $L_F = L_A$

Last, $\Lambda_F = \Lambda_A - c \beta_n I$

$\Lambda_A = \left[ \begin{matrix} -c \beta_n -c \alpha |n| & 0 & 0 & 0 \\ 0 & -c \beta_n & 0 & 0 \\ 0 & 0 & -c \beta_n & 0 \\ 0 & 0 & 0 & -c \beta_n + c \alpha |n| \end{matrix} \right]$

---

If we want to use use volume weighting of our flux by $\frac{1}{\sqrt\gamma}$ then it looks like we will only be adding new terms instead of removing.
Of course the only new term we would produce is associated with $\Psi_i$, where $\nabla_j \Psi_i = \Psi_{i,j} - {\Gamma^k}_{ij} \Psi_k$.
Removing extra terms is typically the case when we look at the covariant derivative form of most equations, where the volume form naturally pops out from the covariant derivative of a vector field: $\nabla_i v^i = {v^i}_{,i} + {\Gamma^i}_{ki} v^k = \frac{1}{\sqrt\gamma} (\sqrt\gamma v^i)_{,i}$

After re-inserting $\frac{1}{\sqrt\gamma} \cdot \sqrt\gamma \cdot ...$ next to the state variable spatial derivatives, we find...

$\left[ \begin{matrix} \Pi \\ \Psi_i \end{matrix} \right]_{,0} + \left[ \begin{matrix} - \beta^k \frac{1}{\sqrt\gamma} ((\sqrt\gamma \Pi)_{,k} - (\sqrt\gamma)_{,k} \Pi) - \alpha \gamma^{jk} \frac{1}{\sqrt\gamma} ((\sqrt\gamma \Psi_j)_{,k} - (\sqrt\gamma)_{,k} \Psi_j) \\ - \alpha \delta^k_i \frac{1}{\sqrt\gamma} ((\sqrt\gamma \Pi)_{,k} - (\sqrt\gamma)_{,k} \Pi) - \beta^k \delta^j_i \frac{1}{\sqrt\gamma} ((\sqrt\gamma \Psi_j)_{,k} - (\sqrt\gamma)_{,k} \Psi_j) \end{matrix} \right] = \left[ \begin{matrix} \Psi_i \alpha_{,j} \gamma^{ij} + K \Pi \alpha - \alpha \Psi_i \Gamma^i - f \alpha \\ \Psi_j {\beta^j}_{,i} + \alpha_{,i} \Pi \end{matrix} \right]$

$\left[ \begin{matrix} \Pi \\ \Psi_i \end{matrix} \right]_{,0} + \frac{1}{\sqrt\gamma} \left[ \begin{matrix} - \beta^k & - \alpha \gamma^{jk} \\ - \alpha \delta^k_i & - \beta^k \delta^j_i \end{matrix} \right] \left[ \begin{matrix} \sqrt\gamma \Pi \\ \sqrt\gamma \Psi_j \end{matrix} \right]_{,k} = \left[ \begin{matrix} K \Pi \alpha - \beta^k {\Gamma^j}_{kj} \Pi + \Psi^j \alpha_{,j} - \alpha \Psi_i \Gamma^i - \alpha {\Gamma^j}_{kj} \Psi^k - f \alpha \\ \alpha_{,i} \Pi - \alpha {\Gamma^j}_{ij} \Pi + \Psi_j {\beta^j}_{,i} - \beta^k {\Gamma^j}_{kj} \Psi_i \end{matrix} \right]$

This just looks more complex.
How about if we weight the time derivatives as well, to see what would happen if we replaced the variables altogether with $\tilde{U} = \sqrt\gamma U$?

$\frac{1}{\sqrt\gamma} \left[ \begin{matrix} \sqrt\gamma \Pi_{,0} + (\sqrt\gamma)_{,0} \Pi - (\sqrt\gamma)_{,0} \Pi \\ \sqrt\gamma \Psi_{i,0} + (\sqrt\gamma)_{,0} \Psi_i - (\sqrt\gamma)_{,0} \Psi_i \end{matrix} \right] + \frac{1}{\sqrt\gamma} \left[ \begin{matrix} - \beta^k & - \alpha \gamma^{jk} \\ - \alpha \delta^k_i & - \beta^k \delta^j_i \end{matrix} \right] \left[ \begin{matrix} \sqrt\gamma \Pi \\ \sqrt\gamma \Psi_j \end{matrix} \right]_{,k} = \left[ \begin{matrix} K \Pi \alpha - \beta^k {\Gamma^j}_{kj} \Pi + \Psi^j \alpha_{,j} - \alpha \Psi_i \Gamma^i - \alpha {\Gamma^j}_{kj} \Psi^k - f \alpha \\ \alpha_{,i} \Pi - \alpha {\Gamma^j}_{ij} \Pi + \Psi_j {\beta^j}_{,i} - \beta^k {\Gamma^j}_{kj} \Psi_i \end{matrix} \right]$

$\left[ \begin{matrix} (\sqrt\gamma \Pi)_{,0} \\ (\sqrt\gamma \Psi_i)_{,0} \end{matrix} \right] + \left[ \begin{matrix} - \beta^k & - \alpha \gamma^{jk} \\ - \alpha \delta^k_i & - \beta^k \delta^j_i \end{matrix} \right] \left[ \begin{matrix} \sqrt\gamma \Pi \\ \sqrt\gamma \Psi_j \end{matrix} \right]_{,k} = \sqrt\gamma \left[ \begin{matrix} K \Pi \alpha - \beta^k {\Gamma^j}_{kj} \Pi + \frac{(\sqrt\gamma)_{,0}}{\sqrt\gamma} \Pi + \Psi^j \alpha_{,j} - \alpha \Psi_i \Gamma^i - \alpha {\Gamma^j}_{kj} \Psi^k - f \alpha \\ \alpha_{,i} \Pi - \alpha {\Gamma^j}_{ij} \Pi + \Psi_j {\beta^j}_{,i} - \beta^k {\Gamma^j}_{kj} \Psi_i + \frac{(\sqrt\gamma)_{,0}}{\sqrt\gamma} \Psi_i \end{matrix} \right]$

Evaluating our volume time derivative:
$(\sqrt\gamma)_{,0} = \frac{1}{2 \sqrt\gamma} \gamma_{,0}$
$= \frac{1}{2 \sqrt\gamma} \gamma \gamma^{ij} \gamma_{ij,0}$
$= -\sqrt\gamma \gamma^{ij} \alpha K_{ij}$
$= -\sqrt\gamma \alpha K$

$\left[ \begin{matrix} (\sqrt\gamma \Pi)_{,0} \\ (\sqrt\gamma \Psi_i)_{,0} \end{matrix} \right] + \left[ \begin{matrix} - \beta^k & - \alpha \gamma^{jk} \\ - \alpha \delta^k_i & - \beta^k \delta^j_i \end{matrix} \right] \left[ \begin{matrix} \sqrt\gamma \Pi \\ \sqrt\gamma \Psi_j \end{matrix} \right]_{,k} = \sqrt\gamma \left[ \begin{matrix} - \beta^k {\Gamma^j}_{kj} \Pi + \Psi^j \alpha_{,j} - \alpha \Psi_i \Gamma^i - \alpha {\Gamma^j}_{kj} \Psi^k - f \alpha \\ \alpha_{,i} \Pi - \alpha {\Gamma^j}_{ij} \Pi - \alpha K \Psi_i + \Psi_j {\beta^j}_{,i} - \beta^k {\Gamma^j}_{kj} \Psi_i \end{matrix} \right]$

$\left[ \begin{matrix} \tilde{\Pi}_{,0} \\ \tilde{\Psi}_{i,0} \end{matrix} \right] + \left[ \begin{matrix} - \beta^k & - \alpha \gamma^{jk} \\ - \alpha \delta^k_i & - \beta^k \delta^j_i \end{matrix} \right] \left[ \begin{matrix} \tilde{\Pi} \\ \tilde{\Psi}_j \end{matrix} \right]_{,k} = \left[ \begin{matrix} - \beta^k {\Gamma^j}_{kj} \tilde{\Pi} + \tilde{\Psi}^j \alpha_{,j} - \alpha \tilde{\Psi}_i \Gamma^i - \alpha {\Gamma^j}_{kj} \tilde{\Psi}^k - f \alpha \sqrt\gamma \\ \alpha_{,i} \tilde{\Pi} - \alpha {\Gamma^j}_{ij} \tilde{\Pi} - \alpha K \tilde{\Psi}_i + \tilde{\Psi}_j {\beta^j}_{,i} - \beta^k {\Gamma^j}_{kj} \tilde{\Psi}_i \end{matrix} \right]$

This still doesn't look much different. We did manage to shift the K source term from $\Pi$ to $\Psi$.

---

Back to our original representation, what if we want to use the contravariant vector $\Psi^i$ as a state variable.

$\Psi^i = \gamma^{ij} \Psi_j$
${\Psi^i}_{,\mu} = (\gamma^{ij} \Psi_j)_{,\mu}$
${\Psi^i}_{,\mu} = {\gamma^{ij}}_{,\mu} \Psi_j + \gamma^{ij} \Psi_{j,\mu}$
${\Psi^i}_{,\mu} = -\gamma^{ij} \gamma_{jk,\mu} \Psi^k + \gamma^{ij} \Psi_{j,\mu}$
spacelike:
${\Psi^i}_{,m} = -\gamma^{ij} \gamma_{jk,m} \Psi^k + \gamma^{ij} \Psi_{j,m}$
${\Psi^i}_{,m} = -\gamma^{ij} (\Gamma_{jmk} + \Gamma_{kmj}) \Psi^k + \gamma^{ij} \Psi_{j,m}$
${\Psi^i}_{,m} = -{\Gamma^i}_{mk} \Psi^k - \Psi^k {\Gamma_{km}}^i + \gamma^{ij} \Psi_{j,m}$
so $\Psi_{j,k} = \gamma_{jm} {\Psi^m}_{,k} + \Gamma_{jkm} \Psi^m + \Psi^m \Gamma_{mkj}$
timelike:
${\Psi^i}_{,0} = -\gamma^{ij} \gamma_{jk,0} \Psi^k + \gamma^{ij} \Psi_{j,0}$
${\Psi^i}_{,0} = 2 \alpha {K^i}_j \Psi^j + \gamma^{ij} \Psi_{j,0}$
so $\Psi_{i,0} = \gamma_{ij} {\Psi^j}_{,0} - 2 \alpha K_{ij} \Psi^j$

$\left[ \begin{matrix} \Pi_{,0} \\ \Psi_{i,0} \end{matrix} \right] = \left[ \begin{matrix} \Pi_{,k} \beta^k + \alpha \Psi_{j,k} \gamma^{jk} + \Psi^j \alpha_{,j} + K \Pi \alpha - \alpha \Psi_i \Gamma^i - f \alpha \\ \alpha \delta^k_i \Pi_{,k} + \beta^k \delta^j_i \Psi_{j,k} + \Psi_j {\beta^j}_{,i} + \alpha_{,i} \Pi \end{matrix} \right]$

...substitute our new derivatives:

$\left[ \begin{matrix} \Pi_{,0} \\ \gamma_{ij} {\Psi^j}_{,0} - 2 \alpha K_{ij} \Psi^j \end{matrix} \right] = \left[ \begin{matrix} \Pi_{,k} \beta^k + \alpha (\gamma_{jm} {\Psi^m}_{,k} + \Gamma_{jkm} \Psi^m + \Psi^m \Gamma_{mkj}) \gamma^{jk} + \Psi^j \alpha_{,j} + K \Pi \alpha - \alpha \Psi_i \Gamma^i - f \alpha \\ \alpha \delta^k_i \Pi_{,k} + \beta^k \delta^j_i (\gamma_{jm} {\Psi^m}_{,k} + \Gamma_{jkm} \Psi^m + \Psi^m \Gamma_{mkj}) + \Psi_j {\beta^j}_{,i} + \alpha_{,i} \Pi \end{matrix} \right]$

$\left[ \begin{matrix} \Pi_{,0} \\ {\Psi^i}_{,0} \end{matrix} \right] = \left[ \begin{matrix} \Pi_{,k} \beta^k + \alpha {\Psi^k}_{,k} + K \Pi \alpha + \Psi^j \alpha_{,j} + \alpha {\Gamma^j}_{kj} \Psi^k - f \alpha \\ \alpha \gamma^{ij} \Pi_{,j} + \beta^k {\Psi^i}_{,k} + \gamma^{ij} \alpha_{,i} \Pi + \gamma^{im} \gamma_{mj,k} \Psi^j \beta^k + \Psi_j {\beta^j}_{,m} \gamma^{im} + 2 \alpha {K^i}_j \Psi^j \end{matrix} \right]$

$\left[ \begin{matrix} \Pi \\ \Psi^i \end{matrix} \right]_{,0} + \left[ \begin{matrix} - \beta^k & - \alpha \delta^k_j \\ - \alpha \gamma^{ik} & - \beta^k \delta^i_j \end{matrix} \right] \left[ \begin{matrix} \Pi \\ \Psi^j \end{matrix} \right]_{,k} = \left[ \begin{matrix} K \Pi \alpha + \Psi^j \alpha_{,j} + \alpha {\Gamma^j}_{kj} \Psi^k - f \alpha \\ \gamma^{ij} \alpha_{,i} \Pi + \gamma^{im} ( \Psi^n \gamma_{nm,k} \beta^k + \Psi^n \gamma_{nj} {\beta^j}_{,m} ) + 2 \alpha {K^i}_j \Psi^j \end{matrix} \right]$

As you can see, changing the state variables from covariant to contravariant only exchanges the delta and the metric terms in the flux. Our new flux Jacobian is just a transpose of our old flux Jacobian.

Now, what if we choose to use a locally orthonormal basis:
$e_I = {e^i}_I e_i = {e^i}_I \partial_i$
$\gamma_{ij} = {e_i}^I {e_j}^J \delta_{IJ}$
$\gamma^{ij} = {e^i}_I {e_j}^I = \delta^i_j$
Where the comma derivative is interchangeable with the application of the basis: $x_i = e_i(x)$ and $x_I = e_I(x) = {e^i}_I \partial_i (x)$
(Notice that I am not normalizing the time component, as the folks who work with ADM tetrads do ... but maybe I should ... ).

$\Psi^I = {e_i}^I \Psi_i$
spatial derivative:
${\Psi^I}_{,J} = {\Psi^I}_{,j} {e^j}_J$
${\Psi^I}_{,J} = (\Psi^i {e_i}^I)_{,j} {e^j}_J$
${\Psi^I}_{,J} = {e_i}^I {\Psi^i}_{,j} {e^j}_J + \Psi^i {{e_i}^I}_{,j} {e^j}_J$
so ${\Psi^i}_{,j} = {e^i}_I {\Psi^I}_{,J} {e_j}^J - \Psi^k {e^i}_I {{e_k}^I}_{,j}$
time derivative:
${\Psi^I}_{,0} = {{e_i}^I}_{,0} \Psi_i + {e_i}^I \Psi_{i,0}$
so $\Psi_{i,0} = {e^i}_I {\Psi^I}_{,0} - {e^i}_I {{e_k}^I}_{,0} \Psi_k$

$\left[ \begin{matrix} \Pi_{,0} \\ {\Psi^I}_{,0} \end{matrix} \right] = \left[ \begin{matrix} \Pi_{,k} \beta^k + \alpha {e^k}_I {\Psi^I}_{,k} + K \Pi \alpha + \alpha_{,j} \Psi^j - \alpha {e^k}_I {{e_m}^I}_{,k} \Psi^m + \alpha {\Gamma^j}_{kj} \Psi^k - f \alpha \\ \alpha {e_i}^I \gamma^{ik} \Pi_{,k} + \beta^k {\Psi^I}_{,k} + {e_j}^I \gamma^{ij} \alpha_{,i} \Pi - \beta^k {{e_m}^I}_{,k} \Psi^m + {e_i}^I \gamma^{im} \gamma_{mj,k} \beta^k \Psi^j + {e_i}^I \gamma^{im} {\beta^j}_{,m} \Psi_j + {{e_k}^I}_{,0} \Psi_k + 2 \alpha {K^I}_J \Psi^J \end{matrix} \right]$

A lot of those source terms can probably cancel if you convert the derivatives of ${e_i}^I$ into 3-connections.
Notice that I'm still keeping the derivative $e_k$ a coordinate derivative, and not converting it to the non-coordinate basis $e_K$.

$\left[ \begin{matrix} \Pi \\ \Psi^I \end{matrix} \right]_{,0} + \left[ \begin{matrix} - \beta^k & - \alpha {e^k}_J \\ - \alpha e^{kI} & - \beta^k \delta^I_J \end{matrix} \right] \left[ \begin{matrix} \Pi \\ \Psi^J \end{matrix} \right]_{,k} = \left[ \begin{matrix} K \Pi \alpha + \alpha_{,j} \Psi^j - \alpha {e^k}_I {{e_m}^I}_{,k} \Psi^m + \alpha {\Gamma^j}_{kj} \Psi^k - f \alpha \\ e^{iI} \alpha_{,i} \Pi - {e^{jI}}_{,k} \beta^k \Psi_j + e^{iI} {\beta^j}_{,i} \Psi_j + {{e_k}^I}_{,0} \Psi_k + 2 \alpha {K^I}_J \Psi^J \end{matrix} \right]$