Possibly a bit more light occurs, if we consider to express that iteration using the matrix-representation of the continued-fractions-evaluation.
Here we insert the coefficients of the cf in a sequence of matrices, whose form is
$ \qquad M_0(a)= \begin{bmatrix} 0 & 1\\1 & a \end{bmatrix}$
and for the product of matrices according to the coefficients in the cf allow the generalized notation
$ \qquad M_0(a,b,c,...,h)= M(a)*M(b)*...*M(h)$
The iteration refers to the partial products, used for a new continued fraction:
$ \qquad M_1(a) = M_0(a); M_1(b) = M_0(a)*M_0(b) = \begin{bmatrix} 1 & b\\a & a*b+1 \end{bmatrix} ; \ldots $
This is simple to program, for instance in Pari/GP. Unfortunately this does not hold exactly: the $M_1()$ created by this procedure are not compatible with the form $ \qquad \begin{bmatrix} 0 & 1\\1 & a \end{bmatrix}$
, they have general (integer) values in all four entries; so do not reflect the required "simple-continued-fraction" representation. So for to model the process as you described it in your question we need some normalization.
What I tried next was to insert the evaluated partial continued fractions as rational (or real) numbers instead of integer numbers in the positions of the a in $M_1(a)$ - so $M_1(a) = M_0(a); M_1(b) = M_0( {ab+1 \over b}), \ldots $ according to the evaluation of the partial convergents. But still this requires then another type of re-normalization (it does not lead to the same limit value), because the fractional parts of the coefficients must now be "shifted" to the remaining cont-frac-expression.
But perhaps we find here the effect, why the leading coefficients of the iterated continued fractions converge to some constant: because after evaluations of subsequent coefficients a certain threshold for the fractional part of a coefficient can no more be overcome (in the next iteration).But I do not see really clear here...