The case where $N=2$ is fairly easy. For $M\ge N$ we'll have an array
$$\begin{array}{ccccc}
a_{1,1} & a_{1,2} & a_{1,3} &\dotsm & a_{1,M}\\
a_{2,1} & a_{2,2} & a_{2,3} &\dotsm & a_{2,M}\\
a_{3,1} & a_{3,2} & a_{3,3} &\dotsm & a_{3,M}\\
&\dotsm & & & \\
a_{N,1} & a_{N,2} & a_{N,3} &\dotsm & a_{N,M}\\
\end{array}$$
Define $t(k,M)$ to be the number of integer solutions of $x_1+x_2+\dotsm+x_m=k$ with all $x_i\ge 0$. The stars and bars theorem tells us that
$$
t(k,M)=\binom{M-1+k}{M-1}
$$
Then if we let $C(N,M)$ denote the number of solutions to the original problem, it's not hard to see that for $N=2$,
$$\begin{align}
C(2,M) &= \sum_{i=0}^M\left[t(i,M)\sum_{j=0}^it(j,M)\right]\\
&= \sum_{i=0}^M\binom{M-1+i}{M-1}\sum_{j=0}^i\binom{M-1+j}{M-1}\\
&= \sum_{i=0}^M\binom{M-1+i}{M-1}\left[\binom{M-1}{M-1}+\binom{M}{M-1}+\dotsm+\binom{M-1+i}{M-1}\right]\\
&= \sum_{i=0}^M\binom{M-1+i}{M-1}\binom{M+i}{M}
\end{align}$$
The last step above uses a fairly well-known result on the sum of a column of Pascal's triangle. So finally we'll have
$$
C(2,M) =\binom{M-1}{M-1}\binom{M}{M}+\binom{M}{M-1}\binom{M+1}{M}+\dotsm+\binom{2M-1}{M-1}\binom{2M}{M}
$$
For $N=M=2$ we have $C(2,2)=1+6+18=25$, which agrees with your example.
Added yet again. For a (not particularly useful) general form, it's not hard to see that
$$
C(N,M)=\sum_{0\le i_1\le i_2\le\dotsm i_N\le M} t(i_1,M)t(i_2,M)\dotsm t(i_N,M)
$$
In this expression, each of the $i_k$ represents the possible sum of the $k$-th row.
I asked on Math.SE whether the expression for $C(2,M)$ had a closed form and the consensus was that it didn't (although it can be expressed as a hypergeometric series, for what that's worth).
