What follows is the work I've done towards a solution, which is unfortunately far too large for a comment. I can make no guarantee that it goes anywhere.
I'll assume $f$ is $C^1$ and strictly convex. Assuming the level sets are compact curves, they must also be simple since a convex function has at most one local extreme. Consider values $t$ in some neighborhood $U$ of $t_0$ (which we will make as small as needed). Let $\gamma_t:[0,1]\to \mathbb R^2$ be a parametrization of $f^{-1}\{t\}$ and note that $\gamma_t$ is locally the graph of a $C^1$ function of $(t,x)$ or $(t,y)$ for fixed $t$. Thanks to compactness and various theorems in calculus, we have a finite set of $a_i,b_i\in [0,1]$ and $C^1$ functions
$$\begin{align}
f_{1}&:U\times (a_1,b_1)\to \mathbb R,\\
&\vdots\\
f_{n}&:U\times (a_n,b_n)\to \mathbb R\\
f_{n+1}&:U\times (a_{n+1},b_{n+1})\to \mathbb R\\
&\vdots\\
f_{n+m}&:U\times (a_{n+m},b_{n+m})\to \mathbb R
\end{align}$$
such that each $a_ii$, by choosing an appropriate parametrization for $\gamma_t$ (arc length is parametrization-independent). Let $L(t)=\ell(f^{-1}\{t\})$ and $a_{n+m+1}=b_{n+m}=1$. Note that
$$\begin{align}
L(t)=\ell(\mathrm{im}(\gamma_t)) &= \int_0^1\left\|\frac{d\gamma_t}{ds}\right\|ds\\
&=\sum_{i=1}^{n}\int_{a_i}^{a_{i+1}} \frac{dx_i}{ds}\sqrt{1+\left(\frac{df_i}{dx_i}\right)^2}ds+\sum_{i=n+1}^{n+m}\int_{a_i}^{a_{i+1}} \frac{dy_i}{ds}\sqrt{\left(\frac{df_i}{dy_i}\right)^2+1}ds
\end{align}$$
thus
$$\begin{align}
\frac{dL}{dt} &= \frac{d}{dt}\left(\sum_{i=1}^{n}\int_{a_i}^{a_{i+1}} \frac{dx_i}{ds}\sqrt{1+\left(\frac{df_i}{dx_i}\right)^2}ds+\sum_{i=n+1}^{n+m}\int_{a_i}^{a_{i+1}} \frac{dy_i}{ds}\sqrt{\left(\frac{df_i}{dy_i}\right)^2+1}ds\right)\\
&= \sum_{i=1}^{n}\int_{a_i}^{a_{i+1}} \frac{d}{dt}\left(\frac{dx_i}{ds}\sqrt{1+\left(\frac{df_i}{dx_i}\right)^2}\right)ds+\sum_{i=n+1}^{n+m}\int_{a_i}^{a_{i+1}} \frac{d}{dt}\left(\frac{dy_i}{ds}\sqrt{\left(\frac{df_i}{dy_i}\right)^2+1}\right)ds
\end{align}$$
which we can analyze piece-by-piece. Looking at the first integrand, note that $x_i$ is affine in $s$ in the region we care about so $\frac{dx_i}{ds}$ is a positive $C^1$ function $c(t)$. Thus we have
$$\begin{align}
\frac{d}{dt}\left(\frac{dx_i}{ds}\sqrt{1+\left(\frac{df_i}{dx_i}\right)^2}\right) &= \frac{d^2x_i}{ds dt}\sqrt{1+\left(\frac{df_i}{dx_i}\right)^2}+\frac{dx_i}{ds}\frac{df_i}{dx_i}\frac{d^2f_i}{dx_i dt}\left(1+\left(\frac{df_i}{dx_i}\right)^2\right)^{-1/2}\\
&= \frac{dc}{dt}\sqrt{1+\left(\frac{df_i}{dx_i}\right)^2}+c(t)\frac{df_i}{dx_i}\frac{d^2f_i}{dx_i dt}\left(1+\left(\frac{df_i}{dx_i}\right)^2\right)^{-1/2}\\
\end{align}$$
which maybe we can sort of understand. Additionally, since the region contained in $f^{-1}\{t\}$ is the intersection of a convex set with a plane, it is convex as well thus $f_i$ (which defines part of its boundary) is concave as a function of $x_i$ so $\frac{df_i}{dx_i}$ is decreasing. We can go back and change some of what we did earlier, moving around our choices of $f_i$ to avoid the points on the curve tangent to the axes, which means we miss a set of points of measure $0$ (which doesn't affect the integral) so $\frac{df_i}{dx_i}$ will have constant sign $\sigma_i$, but this could cause some problems. From here, the best approach is probably to try to understand each of the derivatives. The hardest are those with respect to $t$.