Laplace operator and tensor calculus:

Not from the first principles, but based on physical intuition proof looks as follows. Consider the action for scalar field: $$ S = int d^D x sqrt{g} g^{mu nu} partial_mu phi partial_nu phi $$ This is the only viable GR covariant expression for the action of scalar field without free indices, and $d^D x sqrt{g}$ is an invariant volume element. Integrating this expression by parts, one gets: $$ S = -int d^D x phi partial_mu (sqrt{g} g^{mu nu} partial_nu phi) = - int d^D x sqrt{g} phi frac{1}{sqrt{g}}partial_mu (sqrt{g} g^{mu nu} partial_nu phi) = -int d^D x phi Delta phi $$ Where we have assumed that boundary terms vanis, and recovered in the last equality the invariant volume element.

Well, the Laplace operator is a composite operator:

$$ Delta F = div grad F = nablacdotnabla F $$

and as you wrote

$$ (grad F)^r = (nabla F)^r = frac{partial F}{partial x^k},g^{rk} = V^r $$

You obtain the divergence by contraction of the derivation operator $nabla$ and we emphasize that the contraction has to be performed on the covariant derivative:

$$ div boldsymbol{V} = nabla_iV^i = V^i_{phantom{i};,i}= frac{partial V^i}{partial x^i} + V^r; Gamma^i_{ir} $$

By use of a property of the levi-Civita connection coefficients

$$ Gamma^i_{ki} = frac{1}{2} g^{ij} frac{partial g_{ij}}{partial x^k} = frac{1}{2g} frac{partial g}{partial x^k} = frac{partial ,log sqrt{|g|}}{partial x^k} $$

you can write further

$$ div boldsymbol{V} = nabla_iV^i = V^i_{phantom{i};,i}= frac{partial V^i}{partial x^i} + V^r; Gamma^i_{ir} = frac{partial V^r}{partial x^r} + V^r; frac{partial ,log sqrt{|g|}}{partial x^r} = frac{1}{sqrt{|g|}}; frac{partial}{partial x^r} (sqrt{|g|}; V^r) $$

Finally, substituting $V^r$ gives the desired result:

$$ Delta F = div grad F = frac{1}{sqrt{|g|}}; frac{partial}{partial x^r} (sqrt{|g|}; frac{partial F}{partial x^k},g^{rk} ) $$

Here's a quick derivation of this problem:

As you said: $$Delta F= nabla .nabla F$$ And always using your steps: $$nabla . F=nabla_iv^i$$ And for those who don't know why he involved the "($g^{ik}partial_k F$)" well that's the contravariant components of the gradient operator. begin{align} Delta F&= nabla_iv^i &=nabla_ileft(g^{ik}partial_kFright) &=g^{ik}nabla_ileft(partial_k Fright) end{align}

Recall: $$nabla_i(partial_k F)=partial_{ik}F-Gamma_{ik}^lpartial_lF$$

Hence: $$Delta F=g^{ik}(partial_{ik}F-Gamma_{ik}^lpartial_lF)$$

Another recall: :) $$nabla . mathbf{V}=frac{1}{sqrt{vert gvert}}partial_ileft(v^i sqrt{vert gvert}right)quad{(1)}$$

Involving the contravariant components of $mathbf{grad}F$ in $(1)$ we got the following: $$bbox[silver,5px,border:2px solid teal] {Delta F=frac{1}{sqrt{vert gvert}}partial_ileft(sqrt{vert gvert} g^{ik}partial_kFright)}$$ and that's true, because when $g^{ik}=delta^{ik}$ we get the classic expression of the Laplacien operator: $$Delta F=partial_{kk}F.$$