It is well-known that for a stochastic aperiodic matrix $M$, the sequence $(M^n)_n$ converges.
Here I would like to a have a more precise analysis. Consider now a sequence of stochastic matrices $(M_n)_n$, converging to $M$. We even assume that there exists $0 < \alpha < 1$ such that for all $n$, we have $||M_n - M|| \le \alpha^n$. Is it true that $||M_n^n - M^n||$ converges to $0$?
Since your question is about matrices (finite dimension) there is a simple proof. Suppose your matrices are $N\times N$. The key is to note that any column stochastic matrix has norm one as a map on $\mathbb{R}^N$ with the $\ell^1$ norm ($\|(v_1,\ldots,v_N)^T\|_{\ell^1} = \sum_{j=1}^N |v_j|$). Then, since $$M_n^n -M^n = \sum_{k=0}^{n-1} M_n^k (M_n - M ) M^k$$ and any power of a column stochastic matrix is column stochastic, we obtain $$\|M_n^n -M^n\|_{\ell^1 \rightarrow \ell^1} \ \le \ n \|M_n -M\|_{\ell^1 \rightarrow \ell^1}.$$ If you were asking about the $\ell^1$ to $\ell^1$ norm, this would answer your question. (Note that this works even for infinite stochastic matrices.)
Probably you were asking about the operator norm when we take the Euclidean ($\ell^2$) norm on $\mathbb{R}^N$. However, all norms are equivalent in finite dimensions. Explicitly, $$ \|\vec{v}\|_{\ell^2} \le \|\vec{v}\|_{\ell^1} \le \sqrt{N} \|\vec{v}\|_{\ell^2}$$ and so $$ \frac{1}{\sqrt{N}} \| M \|_{\ell^2 \rightarrow \ell^2} \le \|M\|_{\ell^1 \rightarrow \ell^1} \le \sqrt{N} \|M\|_{\ell^2 \rightarrow \ell^2} $$ for any matrix $M$. Thus $$\|M_n^n -M^n\|_{\ell^2 \rightarrow \ell^2} \ \le \ N n \alpha^n,$$ which goes to zero as you would like.
Chapter 3 of
"Non-negative Matrices and Markov Chains" by E. Seneta, Springer, reprint 2006,
for general non-negative matrices and chapter 4 (especially section 4.3 ff. for the general case) for stochastic matrices seem to be relevant here.
As you have a bound of the form $||M_n - M|| \leq \alpha^n$, where $\alpha \ < 1$, you can probably follow the reasoning beginning on p. 92 there to show, that $M$ has the form $pe^t$ (assuming column stochastic matrices), where $p$ is a probability vector, while $e = (1, \ldots,1)^t$. In that case you find $M^n \ = \ M$.
For all sufficiently large n you can then try to proceed as follows:
As the $M_n$ converge to $M$, their eigenvalue structure will allow to show that $((M_n)^{(n)})^t$ converges to some $\hat{M_n}$ for fixed $n$, and $||((M_n)^n)^t - \hat{M_n}||$ can be bounded in terms of $||((M_n))^t - \hat{M_n}||$.
$||\hat{M_n} - M||$ will be close to zero, maybe under certain extra conditions, because the $M_n$ converge to $M$. The necessary concepts (asymptotic homogenity etc.) seem to be in Seneta.
