I'm going to make a quick argument that won't fill in the details, but sketches the important detail:

Let x = $(x_1, \ldots, x_n)$ and $dx = dx_1 dx_2 \ldots dx_n$.

Suppose you wish to integrate $$ \int_X df \, dx $$ where $X$ is an $(n+1)$-dimensional region. If we let $X_x$ be the one-dimensional region defined by a constant value of $x$, then generalizing Fubini's theorem, we can write this as an iterated integral $$ \int_Y \left( \int_{X_x} df \right) dx $$ where $Y$ is some suitable $n$-dimensional space.

The integral $\int_{X_x} df$ is just the alternating sum of the values $f(P)$ where $P$ iterates over the endpoints of the curves comprising $X_x$, where the upper endpoints are added and the lower endpoints are subtracted. It's convenient to write this as an integral over a zero-dimensional surface: $\int_{\partial X_x} f$.

Consequently, the original integral can be written as $$ \int_Y \left( \int_{\partial X_x} f \right) dx $$ and again essentially by Fubini's theorem, we can identify this with $$ \int_{\partial X} f \, dx $$

Consequently, defining $d(f \, dx)$ as $df dx$ is exactly the right thing to do to generalize the fundamental theorem of calculus to get Stoke's theorem.

To start with 0-forms, $df$ codes how $f$ varies. In fact, it does this in a way that is, IMO, more natural than partial derivatives.

For example, if I want to know how $z = x^2 y$ varies with $x$ — actually that's a lie, I want to know how $z$ varies with $x$ as $y$ is held constant — then I compute an exterior derivative, setting $dy=0$:

$$ dz = 2xy \, dx + x^2 \, dy \equiv 2xy \, dx \pmod{dy} $$

Similarly, as (tangent) vectors are dual to one-forms, we can see that the exterior derivative is the thing you combine with a vector to get a directional derivative.

This is further supported by path integrals; if $\gamma$ is a path from $P$ to $Q$, then $\int_\gamma \, df = f(Q) - f(P)$; so again we see that $df$ is an encoding of how $f$ varies, and the path integral is how we accumulate the variation into a finite difference.

We can argue that $d(df)$ should be zero, ans the variation in the variation of $f$ is second derivative information, and differential forms are only intended to capture first derivative information. Similarly for $df \, df$.

Stokes' theorem expresses the analog of the fundamental theorem of calculus in higher dimensions, giving a way to see the exterior derivative of a differential form as encoding the higher degree variation.

Alternatively, we can appeal to Fubini's theorem to reduce to the one-dimensional case: here's a sketch.

Let x = $(x_1, \ldots, x_n)$ and $dx = dx_1 dx_2 \ldots dx_n$.

Suppose you wish to integrate $$ \int_X df \, dx $$ where $X$ is an $(n+1)$-dimensional region. If we let $X_x$ be the one-dimensional region defined by a constant value of $x$, then generalizing Fubini's theorem, we can write this as an iterated integral $$ \int_Y \left( \int_{X_x} df \right) dx $$ where $Y$ is some suitable $n$-dimensional space.

The integral $\int_{X_x} df$ is just the alternating sum of the values $f(P)$ where $P$ iterates over the endpoints of the curves comprising $X_x$, where the upper endpoints are added and the lower endpoints are subtracted. It's convenient to write this as an integral over a zero-dimensional surface: $\int_{\partial X_x} f$.

Consequently, the original integral can be written as $$ \int_Y \left( \int_{\partial X_x} f \right) dx $$ and again essentially by Fubini's theorem, we can identify this with $$ \int_{\partial X} f \, dx $$

Consequently, defining $d(f \, dx)$ as $df dx$ is exactly the right thing to do to generalize the fundamental theorem of calculus to get Stoke's theorem.