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.