I have rewritten my answer to make it easier to follow and I also fixed some mistakes.
Let $S:[0,1]^2 \rightarrow \mathbb{R}^3$ be some parameterization of your surface (assume that the derivatives do not become zero anywhere, i.e. it's regular/an immersion). That is, the set of points of your surface is $\mathcal{M} = \{S(\beta, \gamma) \,:\, (\beta, \gamma) \in [0,1]^2\}$. Let $h:[0,1]^2\rightarrow \mathbb{R}$ be the heightmap. The unit normal $n:[0,1]^2\rightarrow \mathbb{R}^3$ is given by the normalized cross product of the partial derivatives: $n = \frac{\partial_\beta S \times \partial_\gamma S}{\|\partial_\beta S \times \partial_\gamma S\|}$. Then we displace the surface with the heightmap along the normal in order to construct a new surface $T(\beta, \gamma) = S(\beta, \gamma) + h(u(\beta, \gamma),v(\beta, \gamma))n(\beta, \gamma)$. We want to find the normal of $T$:
$$m = \partial_\beta T \times \partial_\gamma T$$
$$\partial_\beta T = \partial_\beta S+(\partial_\beta h)n + h(\partial_\beta n), \quad \partial_\gamma T = \partial_\gamma S+(\partial_\gamma h)n + h(\partial_\gamma n)$$
$$ \partial_\beta T \approx \partial_\beta S+(\partial_\beta h)n, \quad \partial_\gamma T \approx \partial_\gamma S+(\partial_\gamma h)n$$
$$m\approx (\partial_\beta S+(\partial_\beta h)n)\times(\partial_\gamma S+(\partial_\gamma h)n) = $$
$$\require{cancel}\partial_\beta S\times\partial_\gamma S + (\partial_\beta h)(n\times\partial_\gamma S) + (\partial_\gamma h)(\partial_\beta S\times n) + \cancel{(\partial_\beta h)(\partial_\gamma h)(n\times n)}$$
The above expression gives you an approximation of the (not necessarily unit length) normal of the displaced surface $T$. This is Blinn's approximation - it assumes that $h$ is small enough for $h(\partial_\beta n)$ and $h(\partial_\gamma n)$ to be ignored.
Let us consider how this can be applied in practice. There you usually have a triangle mesh. Then if we know how to modify the normal for every triangle, we know how to modify it for the whole mesh. Consider a specific triangle with vertices $p_0, p_1, p_2 \in \mathbb{R}^3$. We can parametrize the triangle through barycentric coordinates $\alpha, \beta, \gamma \geq 0, \alpha + \beta + \gamma = 1$:
$$S(\beta, \gamma) = \alpha p_0 + \beta p_1 + \gamma p_2 = \\
(1-\beta-\gamma) p_0 + \beta p_1 + \gamma p_2 = p_0 + \beta (p_1-p_0) + \gamma (p_2-p_0).$$
We can now compute the partial derivatives of $S$:
$$\partial_{\beta} S = p_1 - p_0, \, \partial_{\gamma} S = p_2 - p_0,$$
which turn out to be the edges of the triangle. Then the geometric normal is:
$$n_g = \partial_{\beta}S \times \partial_\gamma S, \, \hat{n}_g = \frac{n_g}{\|n_g\|}.$$
We also need $\partial_{\beta} h, \, \partial_\gamma h$ in order to be able to evaluate the perturbed normal. Let the texture coordinates at the 3 triangle vertices be $uv_0, uv_1, uv_2 \in \mathbb{R}^2$, then the reparametrization from barycentric coordinates to UV coordinates is given by:
$$uv(\beta, \gamma) = uv_0 + \beta (uv_1-uv_0) + \gamma (uv_2-uv_0).$$
We can now compute the partial derivatives:
$$\partial_\beta h(uv(\beta, \gamma)) = (uv_1-uv_0) \cdot (\partial_u h(uv(\beta, \gamma)),\partial_v h(uv(\beta, \gamma)))$$
$$\partial_\gamma h(uv(\beta, \gamma)) = (uv_2-uv_0) \cdot (\partial_u h(uv(\beta, \gamma)),\partial_v h(uv(\beta, \gamma)))$$
The finally the perturbed normal is:
$$m_g = \partial_{\beta}S\times\partial_{\gamma}S + (\partial_{\beta} h)\hat{n}_g\times\partial_{\gamma}S + (\partial_{\gamma}h)\partial_{\beta}S\times \hat{n}_g.$$
To compute $\partial_u h, \, \partial_v h$ you can use finite differences:
$$\partial_u h(u,v) \approx \frac{h(u+\epsilon,v) - h(u,v)}{\epsilon}, \, \partial_v h(u,v) \approx \frac{h(u, v+\epsilon) - h(u,v)}{\epsilon}$$
Here is some code to clarify how this would be implemented:
// Assume the following is given for a triangle
vec3 p0, p1, p2; // positions of the 3 vertices
vec2 uv0, uv1, uv2; // texture coordinates at the 3 vertices
// epsilon for the finite differences
float eps;
// we want to compute the normal at bary coords (beta, gamma)
float beta, gamma; // barycentric coordinates
// compute edges
vec3 e1 = p1-p0;
vec3 e2 = p2-p0;
// compute geometric normal
vec3 n_g = cross(e1,e2);
vec3 n_g_n = normalize(n_g);
// compute n x (d_gamma S)
vec3 t_g = -cross(n_g_n, e2);
// compute (d_beta S) x n
vec3 b_g = -cross(e1, n_g_n);
// compute partial derivatives of h
vec2 uv = uv0 + beta * (uv1-uv0) + gamma * (uv2-uv0);
vec2 h_uv = 1/eps * (vec2(h(uv+vec2(eps,0)), h(uv+vec2(0,eps))) - vec2(h(uv)));
float h_beta = dot(uv1-uv0, h_uv);
float h_gamma = dot(uv2-uv0, h_uv);
// compute coefficients
float coef_n = 1.0;
float coef_t = -h_beta;
float coef_b = -h_gamma;
// compute new geometric normal
vec3 m_g = coef_n * n_g + coef_t * t_g + coef_b * b_g;
Note that this matches how normal mapping is sometimes defined, with a sampled normal $(c_n, c_t, c_b) = (1, -\partial_\beta h, -\partial_\gamma h)$:
$$m_g = c_n (\partial_\beta S \times \partial_\gamma S) + c_t (\partial_\gamma S\times\hat{n}_g) + c_b (\hat{n}_g\times \partial_\beta S).$$
You may notice that the cross product arguments have been switched around in the last two term. This is because of how tangent space is defined when doing normal mapping. You will see I have done the same in the code (the minuses involved in $t_g, b_g$). The above doesn't take into account shading normals. When considering those, various approaches have been proposed. I can suggest looking into Simulation of Wrinkled Surfaces Revisited by Mikkelsen for more details on that.
