Mean value theorem for surface integrals: Shrinking surface of integration to $0$ Consider a disk of radius $R$, area $A$, somewhere in $\mathbb{R}^3$. Consider a vector function $\vec{F}: \mathbb{R}^3 \to \mathbb{R}^3$. I form the surface integral over the disk
$$ \iint_{\text{disk}} \vec{F} \cdot d\vec{S} \\
\iint_{\text{disk}} F_ndS
$$
where $F_n$ indicates the normal component to the oriented surface. I applied the Mean Value Theorem for integrals (which I'm assuming works for any integral of any dimension) to get
$$ F_n(P^*)A$$
where $P^*$ is some point on the disk. I form the equation
$$ \frac{1}{A}\iint F_ndS = F_n(P^*)$$
Now I take the limit
$$ \lim_{A\to 0}\frac{1}{A}\iint F_ndS = \lim_{A\to0}F_n(P^*) $$
What value will $F_n(P^*)$ converge to? Will it converge to the value over the center of the disk $F_n(P_0)$? I've done limits along paths in $\mathbb{R}^2$ or $\mathbb{R}^3$. However, I've never done limits of surface areas embedded in $\mathbb{R}^3$. Is such an operation well defined? If the disk shrinks uniformly, I can imagine $F_n(P^*) \to F_n(P_0)$ where $P_0$ is the center of the disk. However, if one side shrinks faster than the other, I can imagine $F_n(P^*)$ converging to any value we like depending on how we shrink the area. I could ask the same question for $1$-dimensional domains. Taking the limit after applying the Mean Value Theorem
$$ \lim_{h\to0} \frac{1}{h}\int_a^{a+h} f(x) dx = \lim_{h\to0} f(c) = f(a)$$ 
However, if I ask 
$$ \lim_{\ell \to 0} \frac{1}{\ell} \int_{a}^b f(x) dx$$
where $\ell = b- a$ is the interval length, is such an operation well defined?
 A: There is a mean value theorem for multiple integrals. For example, if $f:U \subset \mathbb{R}^2 \to \mathbb{R}$ is continuous and $U$ is compact and rectifiable, then there is a point $\xi \in U$, not necessarily unique, such that
$$\int_U f = f(\xi)\cdot \text{area}(U)$$
This is proved in the usual way, noting that if $f$ attains minimum (maximum) values $m (M)$, then $m \cdot \text{area}(U) \leqslant \int_Uf \leqslant M \cdot \text{area}(U).$
It seems you have in mind an integral over a "nice" surface that can be parameterized as $S = \{\mathbf{x} \in \mathbb{R}^3| \,\mathbf{x} = \mathbf{g}(s,t), \,(s,t)\in U \subset \mathbb{R}^2\}$, where $\mathbf{g} \in C^1(U)$ and $U$ is compact and rectifiable. 
The surface integral can be expressed as
$$\int_S F_n \, dS= \int\int_U F_n(\mathbf{g}(s,t)) \left\|\frac{\partial \mathbf{g}}{\partial s} \times \frac{\partial \mathbf{g}}{\partial t} \right\|\, ds\,dt,$$
and if $F_n$ is continuous, the entire integrand is continuous. Whence, there is a point $(\xi,\eta) \in U$ such that
$$\int_S F_n \, dS = F_n(\mathbf{g}(\xi,\eta)) \left\|\frac{\partial \mathbf{g}}{\partial s} \times \frac{\partial \mathbf{g}}{\partial t} \right\|_{(\xi,\eta)}\cdot\text{area}(U)$$
Is "shrinking" the surface is a well-defined operation? The answer is only if you specify the operation as such. For example, construct a decreasing sequence of subsets $S_k$ where $S_{k+1} \subset S_k \subset S$ and $\bigcap_{n=1}^\infty S_k = \{\mathbf{q} \}$.  There will be a corresponding sequence of subsets $U_k \subset U$ and a point $(s_q,t_q) = \mathbf{g}^{-1}(\mathbf{q})$.
Applying the mean value theorem we find a sequence $(\xi_k,\eta_k)$ converging to $(s_q,t_q)$ such that
$$\frac{1}{\text{area}(S_k)} \int_{S_k} F_n \, dS = F_n(\mathbf{g}(\xi_k,\eta_k)) \left\|\frac{\partial \mathbf{g}}{\partial s} \times \frac{\partial \mathbf{g}}{\partial t} \right\|_{(\xi_k,\eta_k)} \cdot \frac{\text{area}(U_k)}{\text{area}(S_k)},$$ 
which should converge as $k \to \infty$ to $F_n(\mathbf{q})$.
To reiterate, the particular value of the limit will depend on how the shrinking process is specified, as illustrated by taking your example further:
$$\lim_{h \to 0} \frac{1}{h}\int_a^{a+h} f(x) \, dx = f(a), \\ \lim_{h \to 0} \frac{1}{h}\int_{b-h}^{b} f(x) \, dx = f(b)$$
A: If the vector field ${\bf v}$ is continuous in a neighborhood of $P_0$ then for any disc $D$ with center $P_0$ and unit normal ${\bf  n}$ one has
$$\int_D{\bf v}\cdot{\bf n}\>{\rm d}\omega= {\bf v}(P_*)\cdot{\bf n}\ {\rm area}(D)$$
for some $P_*\in D$. This follows from the standard MVT for continuous real-valued functions on a path-connected set $D\subset{\mathbb R}^3$. It follows that
$$\lim_{r(D)\to 0+}{1\over{\rm area}(D)}\int_D{\bf v}\cdot{\bf n}\>{\rm d}\omega= {\bf v}(P_0)\cdot{\bf n}\ .$$
But you can also set up an $\epsilon/\delta$-proof that does not make use of the MVT.
