This is version two of my proof. The OP discovered a sign error in my first attempt that revealed my argument to be circular. The correct proof is below.
Not surprisingly, this has to do with the signature of the spacetime metric not being positive definite. Furthermore, this issue is very subtle. It should be noted that this result is quoted partially in Hawking-Ellis, The large-scale structure of spacetime (1973), correctly in Wald, General Relativity (1984), incorrectly in Carroll, Spacetime and Geometry (2003) and hinted at in Straumann, General Relativity (2013). None of these contain a proof. The linked lecture notes have the correct result. Carroll has the "outward" and "inward" directions switched. The book by Gourgoulhon, Special Relativity in General Frames, contains a partial proof in $4$ dimensions, which we adopt here for $n$ dimensions.
We first state the preliminaries. Proofs can be found in the excellent book by Lee, Introduction to Smooth Manifolds (2013). We have provided the theorem numbers for convenience. (At certain points, results in Riemannian geometry must be generalized with some care.)
Let $M$ be a smooth $n$-manifold with boundary $\partial M$. Let $\mathrm{d}$ be the exterior derivative on $M$, $i_X$ the interior derivative wrt. the vector field $X$ and $L_X$ the Lie derivative wrt. $X$. Let $g$ be a pseudo-Riemannian metric on $M$ with Levi-Civita connection $\nabla$ and $(x^\mu)$ be oriented coordinates on $M$. Let $\{E_\mu\}_{\mu=1}^n$ denote an orthonormal frame. The canonical orientation is a top-degree differential form $\mu$ for which $\mu(E_1,\dotsc,E_n)=1$.
We next define "outward" and "inward". A vector $v\in T_pM$ ($v$ should not lie entirely in $T_p\partial M$) is said to point outward if there exists a curve $\gamma:(-\epsilon,0]\to M$ such that $\gamma(0)=p$ and $\dot\gamma(0)=v$. It is said to point inward if there exists a curve $\gamma:[0,\epsilon)\to M$ such that $\gamma(0)=p$ and $\dot\gamma(0)=v$. (Here $\epsilon$ is some positive number.)
Some Results. $\partial M$ is a smooth embedded codimension $1$ submanifold and has a uniquely defined (unit) normal vector. (Thm. 5.11, Prop. 15.33.) The notion of "outward" and "inward" is an equivalence relation on the set of normal vectors. A normal vector $N$ is inward iff $-N$ is outward. (Prop. 5.41, Prop. 15.33.) $\mu$ is a canonical orientation iff $\mu=\sqrt{\lvert\det g_{\mu\nu}\rvert}\,\mathrm{d}x^1\wedge\cdots\wedge\mathrm{d}x^n$. (Prop. 15.31.) If $\omega$ is an orientation, $\iota:S\hookrightarrow M$ is the boundary of a subset of $M$ and $Y$ is a vector that points out of $S$, then $\iota^*(i_Y\omega)$ is a consistently oriented orientation on $S$. (Prop. 15.24) The divergence $\operatorname{div}X:=\nabla_\mu X^\mu$ satisfies $\operatorname{div}X\,\mu=L_X\mu$. (Page 425.) Cartan's formula: $L_X=i_X\circ\mathrm{d}+\mathrm{d}\circ i_X$. (Thm. 14.35.)
Now let $(\mathcal{M},g)$ be the spacetime manifold and metric. Now $\partial \mathcal{M}$ is the boundary and $\iota:\partial \mathcal{M}\hookrightarrow \mathcal{M}$ is the inclusion. See Sect. 2.7 of Hawking-Ellis for a discussion of
Hypersurfaces. If $\iota^*g=:h$ is Riemannian, $\partial\mathcal{M}$ is said to be spacelike and $N$ is timelike. If $h$ is Lorentzian, $\partial\mathcal{M}$ is said to be timelike and $N$ is spacelike.
In the following, let $\tilde\mu$ be the canonical orientation on $(\partial\mathcal{M},h)$ and let $N$ be, for now, the outward pointing unit normal. (Since $g$ is Lorentzian, this means that $\langle N,N\rangle=\pm1$.)
By the Cartan formula, for any vector field $X$,
$$\operatorname{div}X\,\mu=L_X\mu=\{\mathrm{d},i_X\}\mu=\mathrm{d}(i_X\mu),$$
and by Stokes' theorem
$$\int_\mathcal{M}\operatorname{div}X\,\mu=\int_{\partial \mathcal{M}}\iota^*(i_X\mu).$$
It turns out that $\iota^*(i_X\mu)=i_X\mu\restriction_{\partial\mathcal{M}}$ is not so easy to determine.
It is clear if $\partial \mathcal{M}$ is not null, then we have the smooth orthogonal splitting $$T_p\mathcal{M}=T_p\partial \mathcal{M}\oplus \operatorname{span}N,$$
for all $p\in\partial\mathcal{M}$. So let $$ X=X^\top+X^\bot, X^\top\in T_p\partial\mathcal{M},X^\bot\in\operatorname{span}N.$$
Then $X^\bot=\alpha N$ and $X=X^\top+\alpha N$. Taking the inner product with $N$ gives $$\langle N,X\rangle=\alpha\langle N,N\rangle=\sigma\alpha\implies X=X^\top+\sigma\underbrace{\langle N,X\rangle}_\beta N,$$
where $\sigma=+1$ if $N$ is spacelike and $\sigma=-1$ is $N$ is timelike. Thus $$i_X\mu=i_{X^\top}\mu+\sigma\beta i_N\mu.$$
Note now that technically everything we are doing is being pulled back onto $\partial\mathcal{M}$, which is equivalent to the restriction to $\partial\mathcal{M}$. Let $\{E_\mu\}_{\mu=2}^n$ be an orthonormal frame on $\partial \mathcal{M}$. Then $X^\top$ is a linear combination of these vectors, $X^\top=\sum_{\mu=2}^n c_\mu E_\mu$ and
$$i_{X^\top}\mu(E_2,\dotsc,E_n)=\sum_{\mu=2}^nc_\mu\mu(E_\mu,E_2,\dotsc,E_n)=0,$$
due to the total antisymmetry of $\mu$. (If any two slots of a totally antisymmetric tensor are filled by the same vector, it vanishes.) By linearity we can extend this to all vectors in $T_p\partial\mathcal{M}$. Thus $i_{X^\top}\mu\restriction_{\partial \mathcal{M}}\equiv 0$ and we have
$$i_X\mu=\sigma\beta i_N\mu\quad(\text{on $\partial\mathcal{M}$}).$$
Let $\{E_\mu\}_{\mu=2}^n$ be as above, then $\{N,\{E_\mu\}_{\mu=2}^n\}$ is an orthonormal frame on $\mathcal{M}$. (One might have to adjust the sign of an $E_\mu$ so that $\{N,\{E_\mu\}_{\mu=2}^n\}$ is consistently oriented.) Then
$$i_N\mu(E_2,\dotsc,E_n)=\mu(N,E_2\dotsc,E_n)=1,$$
which implies that $i_N\mu$ is the canonical orientation $\tilde\mu$ by Prop. 15.24 and Prop. 15.31 of Lee. Note once again that we take $N$ to be pointing outwards here. Putting everything together, we obtain
$$\int_{\mathcal{M}}\operatorname{div}X\,\mu=\sigma\int_{\partial\mathcal{M}}\langle N,X\rangle\,\tilde\mu.$$
We now have two cases.
Case 1. $N$ is spacelike. Then $\sigma=+1$ and we are done.
Case 2. $N$ is timelike. Now $\sigma=-1$. But then
$$-\int_{\partial\mathcal{M}}\langle N,X\rangle\,\tilde\mu=\int_{\partial\mathcal{M}}\langle -N,X\rangle\,\tilde\mu=\int_{\partial\mathcal{M}}\langle N',X\rangle\,\tilde\mu,$$
where $N'$ is the inward pointing normal vector, by Prop. 5.41 in Lee.
Note that changing the sign of $N$ does not alter $\tilde\mu$, since for it to be consistently oriented we took the normal to be pointing outwards. You could write $\tilde\mu=i_{-N'}\mu$ if you wanted to.
Prop 15.31 in Lee gives $$\mu=\sqrt{\lvert g\rvert }\,\mathrm{d}x^1\wedge\cdots\wedge\mathrm{d}x^n,\quad\tilde\mu=\sqrt{\lvert h\rvert}\,\mathrm{d}y^1\wedge\cdots\wedge\mathrm{d}y^{n-1},$$
where $(y^1,\dotsc,y^{n-1})$ are consistently oriented coordinates on $\partial\mathcal{M}$.
We have thus shown:
$$\int_\mathcal{M}\nabla_\mu X^\mu\sqrt{\lvert g\rvert }\,\mathrm{d}x^1\wedge\cdots\wedge\mathrm{d}x^n=\int_{\partial\mathcal{M}}N_\mu X^\mu\sqrt{\lvert h\rvert}\,\mathrm{d}y^1\wedge\cdots\wedge\mathrm{d}y^{n-1},$$
where $N$ points inward if $\partial\mathcal{M}$ is spacelike and outward if $\partial\mathcal{M}$ is timelike. Note that more generally, one can take the integrals over any $n$-dimensional $\mathcal{U}\subset\mathcal{M}$ with boundary $\partial\mathcal{U}$. This is because $(\mathcal{U},g\restriction_\mathcal{U})$ is a spacetime in its own right, and the proof given above works for $\mathcal{M}$ replaced by $\mathcal{U}$ and $g$ by $g\restriction_\mathcal{U}$.
As for the example on page 456 of Carroll: The region $R$ is bounded by spacelike hypersurfaces, so we must take the normal vector pointing inwards. Then the normal on $\Sigma_2$ is future directed and the normal on $\Sigma_2$ is past directed. But he is defining his integrals $\int_{\Sigma_1}$ and $\int_{\Sigma_2}$ with future directed normal vectors. So the integral $\int_{\Sigma_2}\star J$ picks up a minus sign from the past-directed normal.
