If you just calculate using the moving frame, you'll get the answer for the variation of the principal curvatures in a few lines:
$$
\delta\kappa_i = \mathrm{Hess}(u)(e_i,e_i) + \kappa_i^2\,u .
$$
Here, $\delta\kappa_i$ is the first $t$-derivative of $\kappa_i$ at $t=0$ (i.e., the 'first variation of $\kappa_i$'), $\mathrm{Hess}(u)$ is the quadratic form that is the Hessian of $u$ (using the induced metric on the surface), and $e_1$ and $e_2$ are the principal unit vector fields on the surface (i.e., $e_i$ is tangent to the $i$-th principal curve). In particular, since $H=\tfrac12(\kappa_1+\kappa_2)$ and $K = \kappa_1\kappa_2$, you'll get
$$
\delta H = -\tfrac12\Delta u + (2H^2{-}K)\,u
$$
and
$$
\delta K = -2H\,\Delta u - \mathrm{II}{\cdot}\mathrm{Hess}(u) + 2HK\,u,
$$
where $\Delta u$ is the Laplacian of $u$ (i.e., minus the trace of $\mathrm{Hess}(u)$ with respect to the first fundamental form), and $\mathrm{II}{\cdot}\mathrm{Hess}(u)$ is the inner product of the second fundamental form $\mathrm{II}$ and the Hessian of $u$ (interpreted as symmetric matrices using the orthonormal basis $e_i$, this inner product is just the trace of the product of the two corresponding matrices).
Added at request of the OP:
Sorry that this took so long. Here is how the calculation goes using the moving frame. (As usual, the work is explaining the notation. The calculation itself is easy.)
I'm only going to treat the case of Darboux surfaces, i.e., surfaces in $\mathbb{E}^3$ for which the principal curvatures are distinct. Let $X:\Sigma\times\mathbb{R}\to\mathbb{E}^3$ be a $1$-parameter family of Darboux immersions of an abstract, simply connected, oriented surface $\Sigma$. Let $N:\Sigma\times\mathbb{R}\to S^2$ be the function that gives the oriented unit normal vector field $N(\cdot,t):\Sigma\to S^2$ along the immersion $X(\cdot,t):\Sigma\to\mathbb{E}^3$ for each $t\in\mathbb{R}$. By reparametrizing the family $X$, I can assume that $\partial X/\partial t = u N$ for some (unique) function $u:\Sigma\times\mathbb{R}\to\mathbb{R}$. Since the individual immersions are Darboux, there will exist a smooth mapping $(e_1,e_2,e_3):\Sigma\times\mathbb{R}\to\mathrm{SO}(3)$ such that $e_3 = N$ and such that $e_i(\cdot,t)$ for $i=1,2$ is a Darboux (i.e., principal) frame field along $X(\cdot,t):\Sigma\to\mathbb{E}^3$. Now, setting
$$
\mathrm{d}X = e_1\,\omega_1 + e_2\,\omega_2 + e_3\,\omega_3\,,
$$
one has $\omega_3 = u\,\mathrm{d}t$ and the forms $\omega_1,\omega_2,\mathrm{d}t$ are a basis for the $1$-forms on $\Sigma\times\mathbb{R}$. The Darboux framing assumption is that
$$
\omega_{3i} = \kappa_i\,\omega_i + \mu_i\,\mathrm{d}t
$$
for $i=1,2$, where $\kappa_i$ are the (distinct) principal curvatures.
Now, one has the structure equations $\mathrm{d}e_a = e_b\,\omega_{ba}$, where $\omega_{ba} = -\omega_{ab}$ for $1\le a,b\le 3$ and (summation convention assumed)
$$
\mathrm{d}\omega_a = -\omega_{ab}\wedge\omega_b
\qquad\text{and}\qquad
\mathrm{d}\omega_{ab} = - \omega_{ac}\wedge\omega_{ca}\,.
$$
Now, $\mathrm{d}u = u_1\,\omega_1 + u_2\,\omega_2 + \dot u\,\mathrm{d}t $, so the equation $\mathrm{d}\omega_3 = -\omega_{31}\wedge\omega_1 -\omega_{32}\wedge\omega_2$ becomes
$$
(u_1\,\omega_1 + u_2\,\omega_2)\wedge\mathrm{d}t
= - \omega_{31}\wedge\omega_1 -\omega_{32}\wedge\omega_2
= -\mu_1\,\mathrm{d}t\wedge\omega_1 - \mu_2\,\mathrm{d}t\wedge\omega_2\,,
$$
so it follows that $\mu_i = u_i$. Moreover, since $\mathrm{d}(\mathrm{d}u)=0$, this gives
$$
\mathrm{d}u_1\wedge\omega_1 -u_1\omega_{12}\wedge\omega_2
+\mathrm{d}u_2\wedge\omega_2 -u_2\omega_{21}\wedge\omega_1 \equiv 0 \mod \mathrm{d}t,
$$
so this implies, by Cartan's Lemma, that, for some $u_{ij}=u_{ji}$, I have
$$
\left.
\begin{aligned}
\mathrm{d}u_1&\equiv u_2\omega_{21} + u_{11}\,\omega_1 + u_{12}\,\omega_2\\
\mathrm{d}u_2&\equiv u_1\omega_{12} + u_{12}\,\omega_1 + u_{22}\,\omega_2
\end{aligned}\right\}
\mod \mathrm{d}t
$$
Geometrically, this means that the Hessian of $u(\cdot,t)$ is the quadratic form
$$
\mathrm{Hess}\bigl(u(\cdot,t)\bigr) = u_{11}(\cdot,t)\,{\omega_1}^2
+2u_{12}(\cdot,t)\,\omega_1\omega_2 + u_{22}(\cdot,t)\,{\omega_2}^2.
$$
One also has
$$\mathrm{d}\omega_i
= - \omega_{ij}\wedge\omega_j -\omega_{i3}\wedge\omega_3
= -\omega_{ij}\wedge\omega_j + u\kappa_i\,\omega_i\wedge\mathrm{d}t.
$$
Finally, computing the exterior derivatives of the $\omega_{3i}$ for $i=1,2$ via the structure equation $\mathrm{d}\omega_{3i} = -\omega_{3j}\wedge\omega_{ji}$ yields
$$
\mathrm{d}(\kappa_i\,\omega_i + u_i\,\mathrm{d}t)
= -(\kappa_j\,\omega_j + u_j\,\mathrm{d}t)\wedge\omega_{ji}
$$
(where $j$ is summed over $1,2$). Expanding this, using the above equations, and comparing the coefficients of the $\omega_i\wedge\mathrm{d}t$ term on the RHS and LHS, one obtains the relation
$$
\dot\kappa_i = u_{ii} + {\kappa_i}^2\,u\,,\tag1
$$
where, for $i=1,2$, one has an expansion of the form
$$
\mathrm{d}\kappa_i = \kappa_{ij}\,\omega_j + \dot\kappa_i\,\mathrm{d}t
$$
(where $j$ is summed over $1,2$). Then equation (1) is the desired formula.