The following result is due to Asplund [A, p.235].
Theorem. If $\varnothing\neq E\subset\mathbb{R}^n$ is closed, and $d(x)=\operatorname{dist}(x,E)$, then $f:\mathbb{R}^n\to\mathbb{R}$ defined by $f(x)=|x|^2-d(x)^2$ is convex.
Proof. The proof is very easy, but clever. We have $$ f(x)=|x|^2-\inf_{y\in E}|x-y|^2=|x|^2+\sup_{y\in E}(-|x-y|^2)=\sup_{y\in E} \big(2\langle x,y\rangle -|y|^2\big). $$ Therefore, $f$ is a supremum of a family of affine functions, and hence it is convex. $\Box$
That means $d^2$ has in some sense the same regularity properties as a convex function and that translates into properties of the distance function. The Hessian of a convex function is a positive definite Radon measure and that says something about the Hessian of $d$ when we are at a positive distance to $E$.
If $d>0$, then $|x|^2-f(x)>0$ and hence $d(x)=\sqrt{|x|^2-f(x)}$ is a composition of a function which is the difference of two convex functions with a smooth function $0<y\mapsto\sqrt{y}$ so smoothness properties of convex functions translate into similar properties of $d$. In particular it is known that a convex function conincides with a $C^2$ function outiside a set of small measure (for references see the proof of Theorem 1(b) in [H]) and hence as a corollary we obtain:
Corollary. If $\varnothing\neq E\subset\mathbb{R}^n$ is closed, then for any $\varepsilon>0$ there is a function $g\in C^2(\mathbb{R}^n\setminus E)$ such that $$ |\{x\in\mathbb{R}^n\setminus E:\, f(x)\neq g(x)\}|<\varepsilon. $$
[A] E. Asplund, Cebyvsev sets in Hilbert space. Trans. Amer. Math. Soc. 144 (1969), 235-240.
[H] P. Hajłasz, On an old theorem of Erdős about ambiguous locus. Colloq. Math. 168 (2022), 249–256.