If a domain $\Omega$ has boundary of class $C^k$, $k\geq 2$, then in fact the distance function $d$ to the boundary of $\Omega$ is of class $C^k$ in a neighborhood of the boundary. This is exactly what is proved in Lemma 14.16 in [1] mentioned by OP.
In the case of a convex set we have the following result. Note that we have a better regularity than the one claimed by OP since we do not need to assume that the boundary is of class $C^{k+1}$.
Theorem. Assume that $K\subset\mathbb{R}^n$ is convex with non-empty interior and that $\partial K\in C^k$, $k\geq 2$. Then for any $p>k$ we have $d_K^p\in C^k(\mathbb{R}^n)$.
Proof. Let $\nu:\partial K\to\mathbb{R}^n$ be the unit outer notrmal. If $\partial K$ is represented (locally) as a graph of $\varphi\in C^k$, then $\nu$ can be represented in terms of $D\varphi$ and hence $\nu\in C^{k-1}$.
Consider the function
$$
\Phi:\partial K\times\mathbb{R}\to\mathbb{R}^n,
\quad
\Phi(x,t)=x+\nu(x)t.
$$
Clearly $\Phi\in C^{k-1}$. It follows from the computation in Lemma 14.16 in [1] that the Jacobian $J_\Phi$ of $\Phi$ can be expressed in terms of the principal curvatures and in fact
$$
J_\Phi(x,t)>0
\quad
\text{for all}
\quad
(x,t)\in \partial K\times [0,\infty).
$$
Hence $\Phi$ is a diffeomorphism in a neighborhood of $\partial K\times\{0\}$. However, the normal lines never intersect outside $K$, so it follows that $\Phi$ is actually a diffeomorphism in an open set $U$ that contains $\partial K\times [0,\infty)$.
Clearly $V=\Phi(U)$ is an open subset of $\mathbb{R}^n$ that contains $\mathbb{R}^n\setminus\overset{\circ}{K}$.
Let $\hat{d}_K$ be the signed distance to $\partial K$ which is positive and equals $d_K$ in $\mathbb{R}^n\setminus K$ and is negative in $\overset{\circ}{K}$.
For $y\in V$, let $\pi(y)\in\partial K$ be the unique closest point so
$$
y=\pi(y)+\nu(\pi(y))\hat{d}_K(y),
\quad
\text{i.e.,}
\quad
y=\Phi(\pi(y),\hat{d}_K(y)),
$$
and hence $(\pi(y),\hat{d}_K(y))=\Phi^{-1}(y)$ which proves that $\hat{d}_K\in C^{k-1}(V)$.
This is smaller regularity than we wanted, but a nice trick allows us to show that actually $\hat{d}_K\in C^{k}(V)$. Indeed, $\nabla \hat{d}_K$ points in directional in which the function growths fastest which is $\nu(\pi(y))$. Since the distance growths linearly in that direction, we have that $\nabla \hat{d}_K(y)=\nu(\pi(y))\in C^{k-1}$, and hence $\hat{d}_K\in C^k(V)$.
Now
$$
d_K=\begin{cases}
\hat{d}_K & \text{in } \mathbb{R}^n\setminus\overset{\circ}{K}\\
0 & \text{in } K.
\end{cases}
$$
For that reason reason we lose regularity of the function at the boundary. However, the function $d_K^p$, where $p>k$ has all partial derivatives of order up to $k$ equal zero on the boundary of $K$ and it follows that $d_K^p\in C^k$.
[1] Gilbarg, D. Trudinger, N. S. Elliptic partial differential equations of second order. Reprint of the 1998 edition. Classics in Mathematics. Springer-Verlag, Berlin, 2001.
