The symmetry of the partition function of some square ice models

We consider the partition function $Z(N;x_1,\dots,x_N,y_1,\dots,y_N)$ of the square ice model with domain wall boundary conditions. We give a simple proof that $Z$ is symmetric with respect to all its variables when the global parameter $a$ of the model is set to the special value $a=e^{i\pi/3}$. Our proof does not use any determinant interpretation of $Z$ and can be adapted to other situations (e.g., to some symmetric ice models).