Statistically orthogonal - explanation? I did see the related question here but my question is more related to the actual explanation of the orthogonality itself.
So the following design is orthogonal (this is a latin square to be precise):
1234
2341
3412
4123

As noted in one of the afore mentioned answers:
"An experimental design is orthogonal if the effects of any factor balance out (sum to zero) across the effects of the other factors."
Can that be demonstrated on the design above? I still cannot see how that is orthogonal (and from a comment it should be immediately obvious when you look at it).
 A: I understand you have four row levels, four column levels and the third factor has four levels (1, 2, 3, 4). In total, there are 16 observations.
The model is:
$${\vec Y}_{ijk} = \alpha + \sum_{i=1}^4\alpha_i^{A}{\vec v}_i +  \sum_{j=1}^4\alpha_i^{B}{\vec w}_j +   \sum_{k=1}^4\alpha_k^{C}{\vec z}_k   $$
where $A$, $B$ and $C$ are respectively row, column and the third ("latin letter") factor. Vectors ${\vec v}_i$, ${\vec w}_j$ and ${\vec z}_k$ are such that contain 1 at positions corresponding to combinations of factors present in the design. For instance, you have an experiment in which $A=2$ (second row), $B=2$ (second column) and $C=3$, so in the same position occupied by 
$Y_{223}$ in vector ${\vec Y}_{ijk}$, vectors ${\vec v}_2$, ${\vec w}_2$ and
${\vec z}_3$ would have 1 and all the others 0.
Now if the usual constraints $\sum_{i=1}^4\alpha_i^{A} = \sum_{j=1}^4\alpha_j^{B} = \sum_{k=1}^4\alpha_k^{C} = 0$ are forced upon the coefficients, it is easy to see that any inner product such as
$$< {\vec v}_i, \sum_{k=1}^4\alpha_k^{C}{\vec z}_k>$$
will be necessarily zero. The "ones" in ${\vec  v}_i$ multiply the "ones" in ${\vec z}_k$ only once for each $k$, so:
$$< {\vec v}_i, \sum_{k=1}^4\alpha_k^{C}{\vec z}_k>  = \sum_{k=1}^4\alpha_k^{C}<{\vec v}_i,{\vec z}_k> = \sum_{k=1}^4\alpha_k^{C} = 0$$
the last equality because of the zero sum constraint on the coefficients.
EDIT (in response of comment):
It would require much writing to answer using your example (and would perhaps not fit across the page) so let us imagine a simpler case:
12
21

that is, two row, column, and "latin letter" levels. The response vector and the design matrix would be as follows:
$$
\pmatrix{Y_{111} \\ Y_{122} \\ Y_{212} \\ Y_{221} } 
=
\pmatrix{
    1 & 0 & 1 & 0 & 1 & 0 \cr
    1 & 0 & 0 & 1 & 0 & 1 \cr
    0 & 1 & 1 & 0 & 0 & 1 \cr
    0 & 1 & 0 & 1 & 1 & 0 }
$$
The columns of the matrix are respectively $v_1$, $v_2$. $w_1$, $w_2$, $z_1$ and $z_2$. So you see for instance that facing $Y_{212}$ (=in the same row)  you have "ones" in vectors $v_2$, $w_1$ and $z_2$. That's what I meant.
Let me add that this is something which doesn't lend itself to lengthy explanations online. I suggest you look at some of the (many) good books on this topic. My favourite is Seber's book (I used much the first edition, but there is now a second one). HTH.
Further EDIT:
Consider for instance $\alpha^B_1 = - \alpha^B_2 = 0.4$ so they add up to zero as required by the constraint. Then,
$$ <{\vec v}_1, \sum_{k=1}^4\alpha_k^{B}{\vec w}_k>  =  <\pmatrix{1 & 1 & 0 & 0}, 0.4\pmatrix{1 \\ 0 \\ 1 \\0} - 0.4\pmatrix{0 \\ 1 \\ 0\\ 1}> = 0.
$$
HTH.
