Normal variation of a positive harmonic function

Let $u$ be a positive harmonic function on a Lipschitz domain $\Omega\subset \mathbb{R}^{d+1}$. Following the works of Bourgain ('93 papers in Duke and Matem. Zametki, where $\Omega=\mathbb{R}^2_+$) we show that its variation along normals to the boundary, i.e. the integral
\begin{equation}\notag \displaystyle \int_0^1\left|\frac{\partial}{\partial t} u(\xi+t\vec{N}(\xi))\right|\,dt \end{equation}
where $\vec{N}(\xi)$ is the inward normal vector to $\partial \Omega$ at $\xi$, is finite for many $\xi\in\partial\Omega$.
To be more precise the set of such points turns out to be of full Hausdorff dimension in any boundary ball.
We show this by constructing a family of probability measures $\nu_{\varepsilon}$ that encode the points of bounded normal variation and prove that the local dimension of their supports tends to $d$ as $\varepsilon$ goes to zero. This, in turn, is provided by a careful analysis of a certain differential equation, whose properties are the main theme of the talk.
As an example we give a description of such points for the harmonic measure of a Cantor-type set of positive length on $[0,1]$ – the domain $\Omega$ here is the upper half-plane.