Assuming a metric-compatible connection (which is what we use in GR), then we can write
$$\nabla^2 \phi = g^{ab}\nabla_a \nabla_b \phi = g^{ab}\nabla_a (\partial_b\phi)$$
$$= g^{ab}\bigg(\partial_a\partial_b \phi - \Gamma^c_{ab}\partial_c \phi\bigg)$$
Where we have used the fact that $\phi$ is a scalar field and $\nabla\phi$ is a covector field. Given the Levi-Civita connection coefficients, and metric components, this can be expanded straightforwardly.
To see the conditions under which this formula is valid, note that the coordinate-independent action of the Laplace operator on a scalar field is to take the divergence of its gradient, i.e.
$$\nabla^2 \phi \equiv \mathrm{div}\big(\mathrm{grad}(\phi)\big)$$
The (vector) gradient of a scalar field is a vector field with components
$$\mathrm{grad}(\phi)^a \equiv g^{ab}\partial_b \phi$$
On the other hand, the (covariant) divergence of a vector field is given by $$\mathrm{div}(V) = \nabla_a V^a = \partial_a V^a + \Gamma^a_{a b}V^b$$
Putting those pieces together, we have that
$$\nabla^2\phi = \mathrm{div}\big(\mathrm{grad}(\phi)\big) = \partial_a\big[\mathrm{grad}(\phi)^a\big] + \Gamma^a_{a b} \big[\mathrm{grad}(\phi)^b\big]$$
$$= \partial_a \big[g^{ab}\partial_b \phi\big] + \Gamma^a_{ab} \big[g^{bc} \partial_c \phi\big] = g^{ab}\partial_a\partial_b \phi + \big(\partial_a g^{ab}+ \Gamma^a_{ac} g^{cb}\big)\partial_b\phi$$
where we have sneakily relabeled the indices in the last line. If you know what the metric components and connection coefficients are, then you can just plug them in and go from there. In Cartesian coordinates, the entire second term vanishes (the $g$'s are constant and the $\Gamma$'s are zero), but in general this is of course not true.
If we're working with a metric-compatible connection, then
$$\nabla_a g^{bc} = \partial_a g^{bc} + \Gamma^b_{ad} g^{dc} + \Gamma^c_{ad} g^{bd} = 0$$
for all $a,b,c$. In particular,
$$\nabla_a g^{ab} = \partial_a g^{ab} + \Gamma^a_{ac} g^{cb} + \Gamma^b_{ac} g^{ac} = 0 \implies \partial_a g^{ab} + \Gamma^a_{ac} g^{cb} = - \Gamma^b_{ac} g^{ac}$$
and so our expression can be rewritten
$$\nabla^2 \phi = g^{ab}\partial_a\partial_b \phi - \Gamma^b_{ac} g^{ac}\partial_b \phi = g^{ab} \nabla_a \partial_b \phi = g^{ab}\nabla_a \nabla_b \phi$$
where we've used the fact that $\nabla_b\phi = \partial_b \phi$ for a scalar field.
