Your hypothesis about the application of an algorithmic rule to derive the equation of gradient flow is correct. To contextualize this question, the formulas reported in the OP describe an interesting nonparametric method to reconstruct surfaces from an unorganized point set $S$. This issue has relevant practical applications in a number of scenarios where interpolation of scattered point in three dimensions is needed, including advanced medical imaging (e.g. 3D organ reconstruction from 2D data in computed tomography, nuclear magnetic resonance, and more recently 3D echography), engineer modeling, computer graphics and vision, 3D scanning, and so on. The approaches used in these contexts can be categorized in two main groups. The first group includes parametric methods, where points are usually connected via Voronoi diagrams or Delaunay triangulations to obtain a triangulated/tetrahedral mesh representing a topological embedding of a 2D parameter domain in $\mathbb{R}^3$. These methods, however, are difficult to be applied to complex topologies. Other parametric methods are based on partial differential equations, but again require a strict global parametrization that is hard to be formulated for complex shapes. The second group include nonparametric methods, which substantially attempt to identify a smooth function that applied to the set of interest is close to the zero set. In other words, a scalar function $f$ is defined on a fixed grid to generate an implicit surface that fits the point set usually in progressive steps, with the final shape representing a particular level set of the function. These methods have many advantages over the parametric ones in terms of flexibility, information needed, and ability to represent complex shapes. They may differ with respect to the function $f$ and to the measure used to quantify closeness during the stepwise generation of the final surface.
In the nonparametric method proposed by Osher and Fedkiw, the function $f$ is the signed distance function given by $d(x)= dist(x,S)$, defined as the shortest distance between the element $x$ of our evolving surface and any of the points in the set $S$. On the other hand, as a measure of closeness, the Authors define the "surface energy" using the following functional reported in the OP
$$\displaystyle
E(\Gamma) = \left[ \int_\Gamma d^p(x)\ ds \right]^\frac{1}{p} \quad, \qquad 1 \le p \le \infty
$$
where $\Gamma$ is an arbitrary surface, and the energy is obtained by integrating the distance function over all elements $ds$ of the surface. Note that the parameter $p$ refer to different types of distances (e.g., the Euclidean distance, the so-called Manhattan distance, and so on) that can be employed in this method. In most cases, values of $p=1$ or $p=2$ are typically used. Also note that, from a strictly mathematical point of view, this formulation is substantially a variant of the concept of Dirichlet's energy, a well known measure of function variability that for a generic function $u(x)$ and a constant $1\leq q \leq \infty$ is given by
$$\displaystyle \int_{\Omega}|\nabla u (x)|^q dx$$
where $\nabla u (x)$ is the gradient vector field of the function $u(x)$.
Just as the Dirichlet's energy integral is commonly used in optimization analyses to find a minimizer for given boundary values, the Authors use their definition of energy to find a local minimizer and progressively obtain the final surface model by a process of continuous deformation. In particular, they interpret the evolving surface as a stretched elastic membrane around the points of the set $S$ and use a gradient descent method to identify the best balance between potential energy (distance from set points) and elastic energy (surface area). To achieve this, they propose the following formula to compute the variation in surface energy associated with a small perturbation:
$$\displaystyle\frac{dE(\Gamma)}{d\Gamma} =\frac{1}{p} \left[ \int_\Gamma d^p(x)\ ds \right]^{\frac{1}{p} - 1}
\left[p d^{p-1}(x) \nabla d(x) \cdot N + d^p(x)\kappa \right]$$
The last equation highlights the concept that the variations in surface energy depend on changes in potential and elastic energies, represented by the two terms in the square brackets, respectively. To estimate the change in potential energy, the Authors calculate the product of the derivative of the distance function ($pd^{p-1}(x)$) to its gradient vector field $\nabla d(x)$ and the normal unit $N$ (for the potential energy, we are clearly interested in normal distances). On the other hand, to estimate the change in the elastic energy, they borrow from physical chemistry the concept of surface tension. This is a classical property of liquids arising from unbalanced molecular cohesion forces, as a result of which the liquid surface behaves as an elastic membrane that tends to resist external forces and develops an inward force aimed at minimizing the surface to volume ratio. The mathematics underlying the surface tension phenomenon for real liquids are rather complex, but Osher and Fedkiw take the concept of surface minimization and apply it to their "elastic" surface, assuming that perturbations can lead to a change in the surface area, and then in the elastic energy, that is proportional to both the distance $d^p(x)$ and the curvature $k$. Although this assumption may be convenient from a computational standpoint and rather "intuitive" (the area and the elastic energy of a "stretched membrane" enclosing a 3D point set increase if we deform it with higher distances and more evident curvatures, and conversely decrease at lower distances and lower curvatures that approach the "resting state" of the membrane), it should be pointed out that it represents a considerable simplification of the true dynamics of elastic energy. The definition of surface tension given by the Authors to the quantity $d(x)\kappa$ may also be somewhat misleading, since for a given pressure gradient the "true" surface tension of liquids $T$ is inversely (and not directly) related to the mean curvature $H$ according to the Young-Laplace law $\displaystyle T=\Delta P/(2H)$.
From the formula for the calculation of surface energy variations reported above, the Euler-Lagrange equation for this system is obtained:
$$ \displaystyle d^{p-1}(x) \left[p \nabla d(x) \cdot N + d(x)\kappa \right]=0$$
The LHS corresponds to the quantity in the square brackets of the previous equation, with $ d^{p-1}(x)$ moved outside the brackets. The solution of this differential equation requires a rather complex PDE analysis, but ensures a balance between the potential and elastic energies, which is the fundamental condition for a stationary state. The first term $p \nabla d(x) \cdot N$ minimizes potential energy and moves the surface closer to the set $S$, whereas the second term $d(x)\kappa$ (the "surface tension" mentioned above) minimizes elastic energy and surface area, making the surface stiff at higher distances/curvatures and more flexible at lower distances/curvatures.
On the basis of all considerations above, Osher and Fedkiw finally utilize the well known gradient descent method to generate a progressive, continuous deformation of the surface towards the final shape. The gradient descent method is a technique commonly used to find a local minimum of a function, in which we begin with an initial estimate of the solution, calculate the gradient of the function in that point, take a step in the negative direction of the gradient, and repeat the process until the minimum is reached. This is based on the evidence that a differentiable multivariable function $f(x)$ decreases in a faster way if we take a step in the direction of its negative gradient, i.e. moving from a point $ x=a$ to the point $x=b=a-\alpha\nabla f(x)$, where $\alpha$ is a predefined factor affecting the magnitudes of steps, and $-\nabla f(x)$ is the negative gradient of $f(x)$. Iterating the sequence $x_{n+1}=x_n-\alpha \nabla f(x)$, we get $f(x_{n+1}) \leq f(x_n)$, so that the sequence eventually converges towards the minimum. To make a simple example for a univariable function: considering the function $y=x^2-2x+5$ and a step factor $\alpha=0.1$, we can start with an initial guess of $x_1=3$ (corresponding to $y_1=8$), calculate the negative value of the derivative in this point ($-y'=-(2x-2)$, so $-y'(3)=-4$), and get the next value $x_2=3-4\cdot0.1=2.6$, which corresponds to $y_2=6.56$ and is closer to the true minimum given by $x=1$ and $y=4$. To apply this method to the surface reconstruction problem, the Authors start with an initial surface enclosing all points of the set $S$, and then obtain progressive deformations of the surface following the gradient descent of the energy functional. In this procedure, where the successive surface deformations are captured using the so-called level set method (a technique frequently adopted in image processing, particularly useful for moving interfaces) they define as "gradient flow" the resulting progressive variation of surface energy. This is given by
$$\displaystyle \frac{dE(\Gamma)}{dt} =-\frac{1}{p} \left[\int_\Gamma d^p(x) ds\right]^{\frac{1}{p} - 1}
[ \nabla d^p(x)\cdot N] N$$
which, taking into account the previous expressions given for $\frac{dE(\Gamma)}{d\Gamma}$ and the Euler-Lagrange equation, can be written as
$$\displaystyle \frac{dE(\Gamma)}{dt} =-\frac{1}{p} \left[ \int_\Gamma d^p(x)\ ds \right]^{\frac{1}{p} - 1}\left[p d^{p-1}(x) \nabla d(x) \cdot N + d^p(x)\kappa \right] N$$
and then
$$\displaystyle \frac{d\Gamma}{dt} = -\left[ \int_\Gamma d^p(x)\ ds \right]^{\frac{1}{p} - 1} d^{p-1}(x) \left[ \nabla d(x) \cdot N + \frac{1}{p}d(x)\kappa \right] N $$
which is the equation reported in the OP, and where the minimizer or "stationary" solution is given by the Euler-Lagrange equation. Note that here the concept of "flow" and the use of $dt$ simply refer to the energy variations over the successive temporal steps of surface generation. Also note the minus sign, reflecting the application of the gradient descent method and the progressive decrease in surface energy.
As a final comment, it should be highlighted again that the method proposed by Osher and Fedkiw is only one of the several methods available for surface reconstruction from unorganized data sets. Although it may shows several advantages in terms of flexibility, suitability for complex topologies, smoothness of results, and numerical algorithm efficiency, it also has some limitations (e.g., critical dependency on the choice of the initial surface, risk of collapse, slower computational efficiency of the gradient descent method in "valley" regions).
Best Answer
Since $H^{-1} := (H_0^1)'$ is the dual of a Hilbert space, we have a canonical isomorphism $H_0^1 \simeq H^{-1}$ given by $f^\flat = \langle f, \cdot \rangle_{H^1_0};$ i.e. a function $f \in H^1_0$ corresponds to the linear functional $$f^\flat(\phi)= \int_\Omega \nabla f \cdot \nabla \phi$$ on $H_0^1.$ From integration by parts we see that $f^\flat(\phi) = -\langle \Delta f, \phi\rangle_{L^2};$ so the isomorphism $\flat : H_0^1\to H^{-1}$ can be viewed as the distributional Laplacian $-\Delta.$ This isomorphism is important to us because it is how we define the inner product on the dual space $H^{-1}:$ by definition we have $$\langle f^\flat, g^\flat\rangle_{H^{-1}} = \langle f,g \rangle_{H^1_0},$$ or more explicitly $$\langle u,v \rangle_{H^{-1}} = \langle \Delta^{-1}u,\Delta^{-1}v\rangle_{H_0^1}.$$
The directional derivative of $E$ is $$DE|_u(v) = \frac{d}{ds}\Big|_{s=0} E(u+sv)= \int_\Omega \nabla u\cdot \nabla v=\langle u,v \rangle_{H^1_0}. $$
The $H^{-1}$ gradient of $E$ is defined by the relation $\langle\nabla E|_u, \phi\rangle_{H^{-1}}=DE|_u(\phi);$ i.e. $$\langle \Delta^{-1} \nabla E|_u, \Delta^{-1}\phi \rangle_{H^1_0}=\langle u, \phi \rangle_{H^1_0}.$$ Integrating by parts on the LHS gives $$\langle\Delta^{-2}\nabla E|_u,\phi\rangle=\langle u,\phi\rangle$$ for all $\phi,$ and thus we conclude $\nabla E|_u = \Delta^2 u.$