Assign coordinates to the squares of the grid so that initially, $1$ occupies position $(1, 1)$, $n$ occupies position $(1, n)$, and $n^2$ occupies position $(n, n)$. Clearly any grid configuration that can be reached by the application of some number of operations is given uniquely by two permutations $\pi_r$ and $\pi_c$, where the entry originally at $(i, j)$ is now at position $(\pi_r(i), \pi_c(j))$. If row $i$ initially had average element $(n^2 + 1)/2$ the average element of row $\pi_r(i)$ is $$\sum_{k = 1}^{n}{\#(\pi_r(i), k)} = \sum_{k = 1}^{n}{\#(\pi_r(i), \pi_c(k))}$$ which is the sum of the entries initially in row $i$ so row $\pi_r(i)$ has average element $(n^2 + 1)/2$ as well. Analogous reasoning gives the same conclusion for columns. It remains only to consider diagonals. It's easy to see that initially, both main diagonals have average element $(n^2 + 1)/2$. The sum of the entries on the main diagonals after the operations have been performed is $$\sum_{k = 1}^{n}{\#(\pi_r(k), \pi_c(k))} = \sum_{k = 1}^{n}{n[\pi_r(k) - 1]+ \pi_c(k)} = \frac{n^3 + n}{2}$$ and $$\sum_{k = 1}^{n}{\#(\pi_r(k), \pi_c(n + 1 - k))} = \sum_{k = 1}^{n}{n[\pi_r(k) - 1]+ \pi_c(n + 1 - k)} = \frac{n^3 + n}{2}$$ so that both main diagonals have average element $(n^2 + 1)/2$ as before.

Assign coordinates same as hyperbolictangent. Then we see that we have a special property:$# (i_{1},j_{1}) + # (i_{2}i,j_{2}) = # (i_{2},j_{1}) + #(i_{1},j_{2})$.This property is preserved at each step and it maintains a constant sum of the two main diagonals.Similar reasoning goes for columns.