In Clifford Algebra to Geometric Calculus[3] section 7-2, Hestenes and Sobczyk note that the vector derivative can expressed as
$$
\nabla f(x) = \lim_{\mathscr R\to 0}\frac1{|\mathscr R|I}\int_{\partial\mathscr R}\mathrm d^{n-1}x'\,f(x'). \tag{$*$}
$$
(This is not the exact notation they use, but still the same thing.) This is for an $n$-dimensional space, where $I$ is the unit pseudoscalar, $\mathrm d^{n-1}x$ is the $(n-1)$-vector-valued surface measure, and $\mathscr R$ is meant to be a set containing $x$ with volume $|\mathscr R|$ and boundary $\partial \mathscr R$; the limit is meant to convey the idea that we're shrinking $\mathscr R$ around $x$.
(No, this is not rigorous; I address this somewhat further down.)
As for interpretation, I think "represent[s] the way the field varies in a generalized sense" is exactly what this says. It really is a direct generalization of the one dimensional derivative; when $n = 1$:
$$
\nabla f(x)
= \lim_{l\to 0}\frac1{l\cdot1}\int_{\partial[x-l/2, x+l/2]}\mathrm d^0x'\,f(x')
= \lim_{l\to 0}\frac{f(x+l/2) - f(x-l/2)}l
$$
Here, $\mathscr R$ is an interval $[x-l/2,x+l/2]$, $\partial\mathscr R$ is its end points, and $\mathrm d^0x'$ would be the discrete measure for the end points. We have to remember that $\mathrm d^0x'$ is an oriented measure, and the "normal vector" for $\mathscr R$ at $x-l/2$ is $-1$, and at $x+l/2$ is $1$.
A crucial difference in the general case though is the fact that $\mathrm d^{n-1}x'$ is an $(n-1)$-vector. Using the convention that $\mathrm d^{n-1}x'\,\hat N = |\mathrm d^{n-1}x'|I$ where $\hat N$ is the unit vector normal to $\partial\mathscr R$ at $x'$ (which is actually implicitly assumed by the way ($*$) is formulated), we see that
$$\begin{aligned}
\nabla f(x)
&= \lim_{\mathscr R\to 0}\frac1{|\mathscr R|}\int_{\partial\mathscr R}|\mathrm d^{n-1}x'|\hat Nf(x')
\\
&= \lim_{\mathscr R\to 0}\left[\frac1{|\mathscr R|}\int_{\partial\mathscr R}|\mathrm d^{n-1}x'|\,\hat N{\rfloor} f(x') + \frac1{|\mathscr R|}\int_{\partial\mathscr R}|\mathrm d^{n-1}x'|\,\hat N\wedge f(x')\right].
\end{aligned}$$
Recall for $A$ a blade that $\hat N{\rfloor}A\,A^{-1}$ (where ${\rfloor}$ is the left contraction[2]) is the component of $\hat N$ parallel to $A$, which we can interpret as answering the question "to what degree does $A$ contain $\hat N$?" Also recall that $\hat N\wedge A\,A^{-1}$ is the component of $\hat N$ orthogonal to $A$ in $\hat N\wedge A$, which we can interpret as answering the question "to what degree is $A$ orthogonal to $\hat N$?"
Think of the case that $\mathscr R$ is a sphere centered at $x$, and consider the subspace $\mathscr A$ parallel to $A$ that contains $x'$ (and of course $\hat N$ is normal to $\mathscr R$ at $x'$). Then $\mathscr A$ contains $x$ ("slices through $x$") iff $\hat N\wedge A = 0$ (i.e. $\hat N$ is parallel to $\mathscr A$ and $|\hat N{\rfloor}A|$ is maximal); and $\mathscr A$ is tangent at $x'$ (as $x'$ varies, $\mathscr A$ "wraps around" $x$) iff $\hat N{\rfloor} A = 0$ (i.e. $\hat N$ is orthogonal to $\mathscr A$ and $|\hat N\wedge A|$ is maximal).
So, the $\hat N{\rfloor}f(x')$ integral is telling us "how much is $f$ slicing through $x$ locally on average", and the $\hat N\wedge f(x')$ integral is telling us "how much is $f$ wrapping around $x$ locally on average." When $f$ is vector-valued, this is exactly the divergence and curl (times a pseudoscalar).
Clearly ($*$) is not rigorous; quoting Hestenes and Sobczyk:
The limit is taken by shrinking $\mathscr R$ and hence its volume
$|\mathscr R|$ to zero at point $x$. We allow the limit to be proportional to a delta function or its derivatives, so the limit is well defined in the sense of distribution theory, ref. [GS]. A precise discussion of the limiting process is too involved to be given here. It requires, however, only standard arguments of analysis.
The reference [GS] is
Gelfand, I. and Shilov, G., Generalized Functions, Vol. I, Academic Press, New York (1964).
I have yet to see this definition worked out rigorously anywhere. The issue is really making the limiting process rigorous. Integration seems to be straightforward, at least in flat Euclidean space; see [1,3,5], but do take note of [4].
However, I will demonstrate an interpretation that showcases how it could recapture $\nabla = \sum_{i=1}^n e_i\frac\partial{\partial x_i}$ for an orthonormal basis $\{e_i\}_{i=0}^n$ with coordinates $\{x_i\}_{i=0}^n$.
Let $\mathcal G_n \supset \mathbb R^n$ be an $n$-dimensional Euclidean geometric algebra over $\mathbb R$, i.e. $e_i^2 = 1$ and $e_i\cdot e_j = 0$ for each $i \not= j$. Let
$$\begin{gathered}
\mathscr R(x; s) = \{x' \in \mathbb R^n \;:\; |(x'-x)\cdot e_i| < s/2,\; i = 1,\cdots, n\},
\\
\partial\mathscr R(x; s) = \bigoplus_{i=1}^n S^+_i(x; s)\oplus S^-_i(x; s),\quad S^\pm_i(x; s) = \{x' \in \mathbb R^n \;:\; (x'-x)\cdot e_i = \pm s/2\}
\end{gathered}$$
which is a hypercube centered on $x$ with edge length $s \in (0, 1]$ and each face orthogonal to some $e_i$; the $\oplus$ in the expression for the boundary is to indicate that each component can be integrated separately, since their pair-wise intersections have no area.
We will take $\lim_{\mathscr R(x; s)\to 0} F(\mathscr R(x; s)) = L$ for multivectors $F(\mathscr R), L$ to mean that
$$
\lim_{s\to0^+}||F(\mathscr R(x; s)) - L|| = 0,\quad ||X|| = \sqrt{\langle X\widetilde X\rangle}
$$
where $\langle\cdot\rangle$ is the scalar part projection and $\widetilde X$ is the reverse of $X$. It doesn't matter what norm you choose since $\mathcal G_n$ is a finite-dimensional vector space, and I will only use the facts
$$\begin{gathered}
\lim F = L \implies \lim aFb = aLb,
\qquad
\lim F = L \implies \lim a{\rfloor}F = a{\rfloor}L
\\
\lim(F + G) = \lim F + \lim G,
\end{gathered}$$
for constant vectors $a, b$ and where in the addition formula all limits are assumed to exist. For the contraction formula, note that
$$
a{\rfloor}F = \frac12(aF - \hat Fa),\quad
\langle(\hat F -\hat L)(\hat F - \hat L)\widetilde{}\rangle
= \langle(F - L)(F - L)\widetilde{}\rangle\widehat{}
= \langle(F - L)(F - L)\widetilde{}\rangle,
$$
where $\hat F$ is the grade involution of $F$.
Now let $f : \mathbb R^n \to \mathcal G_n$ be a function such that
$$
\nabla f(x) = \lim_{s\to0^+}\frac1{|\mathscr R(x;s)|I}\int_{\partial\mathscr R(x;s)}\mathrm d^{n-1}x'\,f(x')
$$
exists. Then
$$\begin{aligned}
\nabla f(x)
&= \lim_{s\to0^+}\frac1{s^n}\int_{\partial\mathscr R(x;s)}|\mathrm d^{n-1}x'|\,\hat Nf(x')
\\
&= \lim_{s\to0^+}\frac1{s^n}\sum_{i=1}^n\left[\int_{\partial S^+_i(x;s)}|\mathrm d^{n-1}x'|\,e_if(x') - \int_{\partial S^-_i(x;s)}|\mathrm d^{n-1}x'|\,e_if(x')\right]
\\
&= \lim_{s\to0^+}\frac1{s^n}\sum_{i=1}^ne_i\left[\int_{\partial S^+_i(x;s)}|\mathrm d^{n-1}x'|\,f(x') - \int_{\partial S^+_i(x;s)}|\mathrm d^{n-1}x'|\,f(x'-se_i)\right]
\\
&= \lim_{s\to0^+}\frac1{s^n}\sum_{i=1}^ne_i\int_{\partial S^+_i(x;s)}|\mathrm d^{n-1}x'|\,\bigl(f(x') - f(x'-se_i)\bigr)
\end{aligned}$$
Applying $e_k{\rfloor}$, we see that each limit within the sum exists, so
$$
\nabla f(x) = \sum_{i=1}^ne_i\lim_{s\to0^+}\frac1{s^n}\int_{\partial S^+_i(x;s)}|\mathrm d^{n-1}x'|\,\bigl(f(x') - f(x'-se_i)\bigr).
\tag{$**$}
$$
Let $J_i(x; s)$ be the expression inside the limit on the last line and $F_i(y; s) = f(y) - f(y-se_i)$. Then
$$\begin{aligned}
||J_i(x; s) - f_i(x)||
&\leq ||J_i(x; s) - \frac1s F_i(x; s)|| + ||\frac1s F_i(x; s) - f_i(x)||
\\
&= ||J_i(x; s) - \frac1s F_i(x; s)|| + o(1),
\end{aligned}$$
where $f_i(x) = \partial f(x)/\partial x_i$. The first term can be bounded as
$$\begin{aligned}
||J_i(x; s) - \frac1s F_i(x; s)||
&= ||\frac1{s^n}\int_{S^+_s(x;s)}|\mathrm d^{n-1} x'|\bigl(F_i(x'; s) - F_i(x; s)\bigr)||
\\
&\leq \frac1s\sup_{x'\in S^+_i(x;s)}||F_i(x'; s) - F_i(x; s)||
\end{aligned}$$
Inside the supremum,
$$\begin{aligned}
||F_i(x'; s) - F_i(x; s)||
&= ||f(x') - f(x) + f(x - se_i) - f(x' - se_i)||
\\
&= ||Df(x)[x'-x] + Df(x)[x - x'] + o(|x' - x|)||
\\
&= o(|x' - x|),
\end{aligned}$$
where $Df(x) : \mathbb R^n \to \mathbb R^n$ is the total derivative of $f$ at $x$ and $o(\cdot)$ is little $o$ notation. Since the furthest $x'$ can be from $x$ is the corner of the cube, we get
$$
|x' - x| \leq \frac s2\sum_{i=1}^ne_i = \frac{s\sqrt n}2,
$$
so $||F_i(x'; s) - F_i(x; s)|| = o(s)$, and altogether $||J_i(x; s) - f_i(x)|| = o(1) \to 0$ as $s \to 0^+$. Hence, returning to ($**$) we finally have
$$
\nabla f(x) = \sum_{i=1}^n e_if_i(x) = \sum_{i=1}^n e_i\frac{\partial f}{\partial x_i}(x).
$$
[1] C. Doran, A. Lasenby, Geometric Algebra for Physicists, Cambridge University Press, Cambridge, 2003.
[2] L. Dorst, The Inner Products of Geometric Algebra, in: L. Dorst, C. Doran, J. Lasenby (Eds.), Applications of Geometric Algebra in Computer Science and Engineering, Birkhäuser, Boston, 2002: pp. 35–46.
[3] D. Hestenes, G. Sobczyk, Clifford algebra to geometric calculus: a unified language for mathematics and physics, D. Reidel, Dordrecht, 1984.
[4] A. Macdonald, Sobczyk’s simplicial calculus does not have a proper foundation, arXiv:1710.08274 [Math]. (2017).
[5] G. Sobczyk, O.L. Sánchez, Fundamental Theorem of Calculus, Adv. Appl. Clifford Algebras. 21 (2011) 221–231.
