There is the geometric definition of a divergence of a vertor field to be the following limit:

How does this definition turns out to be the del operator dot the vector field in cartesian coordinates?
Thanks, Eran.
I tried the calculations for the $ z $ component of the field and got the integrals:
$ \int_{y-\frac{L}{2}}^{y+\frac{L}{2}}\int_{x-\frac{L}{2}}^{x+\frac{L}{2}}F_{z}(x,y,z+L)dxdy-\int_{y-\frac{L}{2}}^{y+\frac{L}{2}}\int_{x-\frac{L}{2}}^{x+\frac{L}{2}}F_{z}(x,y,z)dxdy= $
$ =\int_{y-\frac{L}{2}}^{y+\frac{L}{2}}\int_{x-\frac{L}{2}}^{x+\frac{L}{2}}[F_{z}(x,y,z+L)-F_{z}(x,y,z)]dxdy $
Where of course, $L$ goes to $0$.
What can I do next?
Now that you're at this point, you need to approximate $F_z(x, y, z+L)$ in terms of $F_z(x, y, z)$ and $L$.
For instance, consider a Taylor series expansion:
$$F_z(x, y, z+L) = F_z(x, y, z) + L \left. \frac{\partial F_z}{\partial z} \right|_{x, y, z} + \frac{1}{2} L^2 \left. \frac{\partial^2 F_z}{\partial z^2} \right|_{x, y, z} + O(L^3)$$
This lets you argue that the integral on the $z$-faces takes the form
$$\int_{x-L/2}^{x+L/2} \int_{y-L/2}^{y+L/2} L \left. \frac{\partial F_z}{\partial z} \right|_{x, y, z} + O(L^2) \, dx \, dy$$
Use more Taylor expansions (on $x$ and $y$) to argue that, to lowest order in $L$, we get
$$\left. L^3 \frac{\partial F_z}{\partial z} \right|_{x, y, z}$$
Of course, the volume of the cube is $L^3$, so when we divide and take the limit, we get
$$\lim_{L \to 0} \frac{1}{L^3} L^3 \frac{\partial F_z}{\partial z} = \frac{\partial F_z}{\partial z}$$
I trust you can repeat the process for the $F_x, F_y$ components?