Let's have a surface $\Sigma(t)$ moving in the three-dimensional space $\mathbb{R}^3$, with boundary $\partial \Sigma$ as in the below figure:
Let's have a space-dependent scalar $c(\mathbf{x})$, with $\mathbf{x}$ the three-dimensional position vector.
I use the Leibniz rule to compute:
$$\dfrac{\mathrm{d}}{\mathrm{d} t} \int_{\Sigma(t)} c(\mathbf{x}) \mathrm{d} \Sigma = \int_{\Sigma(t)} \dfrac{\partial c(\mathbf{x})}{\partial t} \mathrm{d} \Sigma + \int_{\partial\Sigma} (c \mathbf{v}_b) \cdot \mathbf{n} \ \mathrm{d} s$$
The second term is null because the scalar field $c$ is fixed in time. If also the boundary $\partial \Sigma$ is fixed $\mathbf{v}_b=0$ and:
$$\dfrac{\mathrm{d}}{\mathrm{d} t} \int_{\Sigma(t)} c(\mathbf{x}) \mathrm{d} \Sigma = 0 \ ,$$
which is clearly not the case for my expanding moving surface. Where am I wrong?
Comparing your expression with the "Three-dimensional, time-dependent case" on wiki (resp.): $$\color{blue}{\frac{d}{dt}\int_{\Sigma(t)}c(\boldsymbol{x})d\Sigma = \int_{\Sigma(t)}\frac{\partial c(\boldsymbol{x})}{\partial t}d\Sigma + \int_{\partial\Sigma}(c\boldsymbol{v}_b)\cdot \boldsymbol{n}ds}$$ $$\text{vs.}$$ $$\frac{d}{dt}\int\int_{\Sigma(t)}\boldsymbol{F}(\boldsymbol{r},t)\cdot d\boldsymbol{A}=\int\int_{\Sigma(t)}(\boldsymbol{F}_t(\boldsymbol{r},t)+[\nabla\cdot\boldsymbol{F}(\boldsymbol{r},t)]\boldsymbol{v})\cdot d\boldsymbol{A}-\oint_{\partial\Sigma(t)}[\boldsymbol{v}\times\boldsymbol{F}(\boldsymbol{r},t)]\cdot d\boldsymbol{s}$$
It seems we should have: $$\color{red}{\boldsymbol{F}(\boldsymbol{r},t)}:= c(\boldsymbol{x}),\text{ }\text{ }\color{red}{d\boldsymbol{{A}}}:=\boldsymbol{n}(\boldsymbol{x})d\Sigma, \text{ }\text{ and }\text{ }\text{ }\color{red}{d\boldsymbol{s}}:=\boldsymbol{n}(\boldsymbol{x})ds$$ with $$(\boldsymbol{F}_t(\boldsymbol{r},t)+[\nabla\cdot \boldsymbol{F}(\boldsymbol{r},t)]\boldsymbol{v})\cdot d\boldsymbol{A} = \color{blue}{\frac{\partial c(\boldsymbol{x})}{\partial t}\boldsymbol{n}(\boldsymbol{x})d\Sigma}\text{ }\text{ }\text{ }\text{ }(\star1)$$ and $$-[\boldsymbol{v}\times \boldsymbol{F}(\boldsymbol{r},t)]\cdot d\boldsymbol{s}= \color{blue}{(c(\boldsymbol{x})\boldsymbol{v}_b)\cdot \boldsymbol{n}(\boldsymbol{x})ds}.\text{ }\text{ }\text{ }\text{ }(\star2)$$ The vector velocities: $\boldsymbol{v}$ and $\boldsymbol{v_b}$ are of the surface and the boundary (I think these are spatially dependent). Annoyingly, I think also $\boldsymbol{r}$ is in place of $\boldsymbol{x}$.
Note: I added the $n(\boldsymbol{x})$ into your $d\boldsymbol{A}$ terms, since the integral represents flux of a field through the surface, we are measuring the component (dot product) of the field aligning with the normal. This (unit normal) is scaled by the volume element in the integrand.
But, with the given definitions $\color{red}{(red)}$, calculating the left hand sides gives: $$(\star1) = \bigg(\frac{\partial}{\partial_t}c(\boldsymbol{x})+[\nabla\cdot c(\boldsymbol{x})]\boldsymbol{v}(\boldsymbol{x})\bigg)\cdot \boldsymbol{n}(\boldsymbol{x})d\Sigma = \color{purple}{[\nabla\cdot c(\boldsymbol{x})]\boldsymbol{v}(\boldsymbol{x})\cdot \boldsymbol{n}(\boldsymbol{x})d\Sigma}$$ and $$(\star2) = \color{purple}{-[\boldsymbol{v_b}(\boldsymbol{x})\times c(\boldsymbol{x})]\cdot \boldsymbol{n}(\boldsymbol{x})ds}.$$
The RHS's differ by one term (one zeroes out) and the dot term is left. This raises a question of how to take dot and cross product with a vector field and a scalar field (do we reduce to multiplication?). If this is the case, then by your observation with $v_b=0$ then the only surviving term would be: $$\frac{d}{dt}\int_{\Sigma(t)}c(\boldsymbol{x})d\Sigma = \int_{\Sigma(t)}[\nabla\cdot c(\boldsymbol{x})]\boldsymbol{v}(\boldsymbol{x})\cdot \boldsymbol{n}(\boldsymbol{x})d\Sigma$$