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_i, the collection $(a_i,b_i)$ covers $(0,1)$, and there exist $C^1$ functions $x_i,y_i$ which are strictly increasing in $s$ such that $t\in U,s\in (a_i,b_i)\implies \gamma_t(s)=\begin{cases} (x_i(t,s),f_i(t,x_i(t,s))) &\text{if } 1\leq i\leq n\\ (f_i(t,y_i(t,s)),y_i(t,s)) &\text{if } n For notational convenience, assume $0=a_1< a_2< b_1 i.e. that the intervals are in order, nicely trimmed and non-redundant, which saves us from having to through out/trim intervals and rearrange junk. Furthermore, we can assume that $x_i,y_i$ are affine functions in $s$, at least where their domain does not overlap with that of $x_j,y_j$ for $j>i$, 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$.