<PERSON>'s equations in curved spacetime are written in the form $$\begin{split}\nabla_a F^{ab} &= - 4\pi J^b,\ \nabla_{[a} F_{bc]} &= 0,\end{split}$$ with $F$ the <PERSON> two-form, $J^a$ the current four-vector, $\nabla$ the covariant derivative and $[]$ denotes antisymmetrization of the indices. In terms of exterior calculus they become: $$ \begin{split} d\star F &= 4\pi \star J\ dF&=0,\end{split}$$ with $\star$ the <PERSON> dual, which sends p-forms to $4-p$-forms in dimension 4. If we integrate the left side of the first equation over a space-like hypersurface of dimension 3, $\Sigma$, with normal time-like vector $t^a$, then <PERSON>' theorem yields $$\int_\Sigma d\star F = \int_S \star F,$$ with $S$ the boundary of $\Sigma$ with normal $n_a$. Since $\Sigma$ is space-like and $(\star F){cd} = \frac{1}{2} F^{ab} \epsilon{abcd}$, $S$ is also space-like and one component of $F$ must be time-like, therefore $\star F = F^{ab} t_a n_b dS$. 578
This is easy to see also if one takes the restriction of the dual on $S$ in local coordinates, all the 2-forms are space-like. But we know that $E_a = F_{ab} t^b$, hence $$\int_S \star F = -\int_S F^{ba} t_a n_b d S = -\int_S E_b n^b d S.$$ Now we integrate the right side, $$\int_\Sigma \star J = \int_\Sigma J^a \epsilon_{abcd}=\int_\Sigma J^a t_a d \Sigma = -q,$$ and after combining this and the previous one, we obtain: $$\int_S E_a n^a d S = 4\pi q,$$ <PERSON>'s law applies also in curved spacetime. Note that $\epsilon_{abcd}$ is the Levi-Civita (volume) tensor, not the symbol. In local coordinates its components are the product of the symbol with $\sqrt{|\det{g_{\mu\nu}}|}$.
For the case of the magnetic field, $B_a = - \frac{1}{2} \epsilon_{abcd} F^{cd} t^b$, only the space-like components of $F_{ab}$ are used, and the magnetic field part of the <PERSON> tensor is $$F = F_{12} dx^1\wedge dx^2 + F_{23} dx^2\wedge dx^3 + F_{31} dx^3\wedge dx^1,$$ and the components of the field are $B^1 = - |\det g_{\mu\nu}|^{-1/2} F_{23}$ etc, therefore the integrals become $$0=\int_S F = -\int_S \sqrt{|\det g_{\mu\nu}|} (B^1 dx^2\wedge dx^3 + B^2 dx^3\wedge dx^1 + B^3 dx^1\wedge dx^2) = - \int_S B^a n_a dS.$$ 578