So I wanted to return to this and add another answer about some of what happens when we exclude $0$ from the natural numbers (and $n\mapsto 0$ from the starting functions).

In this case, this set of functions, ordered under $\preceq$, is in fact a well-quasi-order!  I don't know how to show it's a total order, which would require showing both antisymmetry and totality.  But we can show it's at least a well-quasi-order, and get an upper bound on its maximum linearization, aka type (which if it *is* well-ordered would then be a bound on its order type), without a lot of difficulty, by applying Kruskal's tree theorem and the habilitation thesis of Diana Schmidt.

Specifically, we can think of function specifications here as being encoded by ternary trees; each vertex either has degree 3 or degree 0, the degree 3 leaves are unlabeled, and the degree 0 leaves are labeled with one of $1$, $\mathrm{id}$, or $S$.  Actually, to make the argument work, we'll allow labeling leaves with any constant, not just $1$; I'm hopeful there's a way to make it work without that, but at the moment I don't see how to.  I'll order the labels by $\preceq$, so they form an $\omega+2$.

We then want to, using the order on trees from Kruskal's tree theorem, map this set of trees onto this set of functions in a monotonic way.  If we can do so, then the set of functions will be a well-quasi-order, and its type will be at most that of our set of trees.  The type of the latter, meanwhile, can be upper bounded using Schmidt's thesis.  Specifically -- just applying the results directly in the most naïve way -- we get an upper bound of $f^+ \left(\begin{array}{cc}\omega+2 & 1 \\0 & 3\end{array}\right)$ (in Klammersymbol notation), which I *think* offhand is equal to the ordinal $\phi_{1,0,0}(\omega+2)$ but I'm doing this quickly so I'm not certain of that.  In reality I suspect this bound is overly large and one can do better, but I'm going with the easy bound for now!

The problem, of course, is showing that the mapping is monotonic -- and this is a problem because in fact it *isn't*.  However, we can restrict to a subset on which it is, and which will still get us all of our functions.  I claim that we can do that in a way that makes the restricted mapping monotonic.

In order to verify this, we'll first need to prove some basic facts about the functions we can get.  Specifically:

1. Every function we in this set is monotonic.
2. Every function in this set is either constant, or satisfies $f(n)\ge n$ everywhere.  Moreover, in the second case, either $f$ is the identity, or it satisfies $f(n)>n$ everywhere.

Both these statements are of course true of our starting functions $0$ (or any constant), $\mathrm{id}$, and $S$.  Now, if we have three such functions in our set $f$, $g$, and $h$, we can check that they're true for $x\mapsto f^{g(x)}(h(x))$.

Specifically: If $f$ is constant, then (because $g(x)\ge 1$!) we just get $f$ back, so we get back a constant, which satisfies both constraints.  So assume $f(n)\ge n$.  Then, since $f$ is monotonic, $h$ is monotonic, and $g$ is monotonic and $f(n)\ge n$, we can conclude that our new function is monotonic.  Now we need to check the second constraint.

We get the following cases:
1. $f$ is constant -- already checked above, we get back $f$ (using the fact that $g(n)\ge 1$).
2. $g$ and $h$ both constant -- we get another constant.
3. $f=\mathrm{id}$, $h$ is constant -- we get back $h$, again a constant.
4. $f(n)>n$, $g(n)\ge n$ -- Since $f$ is monotonic, if we use the fact that $h(n)\ge 1$, we get $f^{g(n)}(h(n))\ge S^n(1)=n+1>n$.
5. $f(n)\ge n$, $h(n)\ge n$ -- Since $f$ is monotonic, we get $f^{g(n)}(h(n))\ge n$.  If $f=h=\mathrm{id}$ we get the identity function; if either $f(n)>n$ or $h(n)>n$ we get $f^{g(n)}(h(n))>n$

You can check that these are exhaustive, and that therefore the two constraints above always hold.

OK!  So now let's restrict our construction by disallowing the following cases:

1. $f$ is a constant.  In this case we just get back $f$, which we already had.
2. $f$ is the identity.  In this case we just get back $h$, which we already had.
3. $g$ and $h$ are constants.  In this case we will just get back another constant, though it may be different from $g$ and $h$.  This case is why I allowed arbitrary constants... if we only have $1$, and we bar this case, we can't make any other constants!

So these restrictions haven't cost us any surjectivity.  And, I claim, we make our map monotonic!

So, let's check the monotonicity of our map.  There are two things we have to check.  The first is that replacing $f$, $g$, or $h$ by a larger function (under $\preceq$) will make our function larger (under $\preceq$).  This holds without any restrictions.  For $h$ this holds because $f$ is monotonic; for $g$ this holds because $f$ is monotonic and $f(n)\ge n$ (and if we allow $f$ constant, it still holds because we just get $f$); and for $f$ this holds because $h$ is monotonic and all our allowed $f$ are monotonic as well.

The more problematic thing is that we need to check that $f^{g(x)}(h(x))\succeq f, g, h$.   But our restrictions will make this true.

First off, since $f(n)\ge n$, $f$ is monotonic, and $g(n)\ge 1$, we see that $f^{g(x)}(h(x))\ge h$.  Also, since $h(n)\ge 1$ and $f(n)>n$, we have $f^{g(x)}(h(x))\ge S^{g(n)}(1)=g(n)+1>g(n)$.

Finally we need to check that $f^{g(x)}(h(x))\succeq f$.  Since $g$ and $h$ are not both constant, either $h(n)\ge n$, or $g(n)\ge n$.  If $h(n)\ge n$, then the claim follows by the monotonicity of $h$ and $f$ (and the fact that $g(n)\ge 1$).  If $g(n)\ge n$, meanwhile, then (applying monotonicity of $f$, etc) $f^{g(n)}(h(n))\ge f^n(1)$, so it suffices to show that $f^n(1)\ge f(n)$.  But $f^n(1) = f(f^{n-1}(1)) \ge f(S^{n-1}(1)) = f(n)$, as required.

That completes the proof.

So:
1. This set of functions is at least well-quasi-ordered under $\preceq$; whether it is antisymmetric or total I do not know, but it is at least well-quasi-ordered.
2. The type (maximum linearization), and therefore the order type if it is well-ordered, is at most $f^+ \left(\begin{array}{cc}\omega+2 & 1 \\0 & 3\end{array}\right)$; this is probably a substantial overestimate, but I prioritized putting someting up here rather than optimizing the bound.

I realize this is many years later, but I hope this is helpful! :)