I'm not sure if this point of view is taken up in the many references which are named here, but I'll say something about an "elementary" way to discover the exterior derivative which sounds like ordinary calculus. Let's take on the point of view that a $k$-form is something you integrate over a $k$-dimensional submanifold. If you imagine $k$-dimensional submanifolds as being composed of a $k$-dimensional blanket of little $k$-parallelograms, then this is a geometrically natural point of view since the $k$-form will assign a (small) number to each of these parallelograms. To actually realize a submanifold as such a "blanket" is to give a parameterization. (These parallelograms are oriented; this picture is different from surface integration of scalar functions in Riemannian geometry where one simply imagines some distribution of mass on the manifold and the integral is completely measure-theoretic.)
In one-variable calculus, when $f$ is a function, $df$ tells you the change in $f$ per small change in its input, and if you integrate it over a curve from $a$ to $b$, it expresses the total change in $f$ from $a$ to $b$. Now, a one form $\eta$ is integrated not over points but rather over curves. Still, you can ask, how does $\int_\gamma\eta$ change when you perturb $\gamma$? Well, if you deform a closed curve $\gamma_a$ into another curve $\gamma_b$, the difference between the integrals over $\gamma_b$ and $\gamma_a$ is some derivative we can call "$d\eta$" integrated over the surface swept between the two.
Picturing the case where $\gamma_a$ and $\gamma_b$ bound an annulus is a good thing to consider here; this interpretation tells you how to orient the boundary of the annulus if you want to think of $\int_\Sigma d\eta = \int_{\gamma_b} \eta - \int_{\gamma_a} \eta$ as being $\int_{\partial \Sigma} \eta$. On the other hand, you can take the point of view that the orientation for $\Sigma$ is determined by the requirement that we start at $\gamma_a$ and go to $\gamma_b$ (much like the case for $df$ of a function). You can then contract the inner circle to a point to recover Stokes' theorem for a disk -- the integral over the inner circle will vanish in the limit by the linearity and continuity of the form.
It's not completely necessary that the curve (or $k$-dimensional submanifold) you deform is closed, but by rule the boundary should remain fixed during the deformation or you will miss out on part of the boundary.
Using a specific example like a square/cube, we can get a coordinate representation for $d\eta$ through the fundamental theorem of calculus. (For $0$ forms, every point is closed, so we did not need to worry about the word "closed" before.)
It is easy to see many properties. For example, let's take $\eta$ to be a $1$-form in $3$-space; then $d^2 \eta$ is clearly $0$. Let $\gamma$ be a circle, and let $\Sigma_a$ and $\Sigma_b$ be the upper and lower hemispheres of a ball $B$ whose equator is $\gamma$. Then $\int_{\Sigma_a} d \eta = \int_\gamma \eta = \int_{\Sigma_b} d \eta$ by Stokes' theorem for a disk. On the other hand, the integral of $d^2\eta$ over the ball $B$ is just $\int_{\Sigma_b} d \eta - \int_{\Sigma_a} d\eta = 0$ because you can sweep out $B$ by deforming $\Sigma_a$ to $\Sigma_b$ with the boundary fixed. Since $\int_B d^2 \eta = 0$ for every ball, $d^2 \eta$ is identically $0$. When you execute this proof for a square, you see that mixed partials commute.
I would like to know if the product rule can easily be seen through this interpretation, but I have not thought enough about it to see it clearly yet.