Toroidal alternating sign matrices Consider $n\times k$ matrices with entries from $\{0,1,-1\}$ such that the sum in each row and each column is 0 and the non-zero numbers in each row/column alternate in sign (so, they alternate if we make a torus from an $n\times k$ table). Were such things studied? Enumerated?
 A: Recall Kuperberg's proof of the ASM conjecture [arXiv:math/9712207] via a relation to a quantum-integrable two-dimensional lattice model from classical statistical mechanics, namely a special case of the six-vertex model with domain-wall boundaries. See also the excellent book by Bressoud, Proofs and Confirmations (The Story of the Alternating-Sign Matrix Conjecture).
The same mapping will relate the problem you're interested in to a special case of the six-vertex model with toroidal boundary conditions. Let me denote the three weights of the (symmetric) six-vertex model by $a,b,c$ as usual for this model. For ASM it suffices to consider the special case with vertex weights $a=b$ and $c$ free, which is also known as the F-model in the physics literature. Its partition function for toroidal boundaries is a generating function for $x$-enumerations, including $x=1$ (counting the number of ASM) and $x=2$ (for domino tilings), of the ASM you're after.
Its partition function was found exactly long ago by Lieb [PRL 18 (1967) 1046]. The main difference with the case of usual ASM (or domain-wall boundaries in terms of the six-vertex model) is that for finite size the result is not given in a closed form, but only up to solving systems of algebraic equations known as the "Bethe-ansatz equations". For sufficently low size this could be done numerically. The large-size asymptotics of the partition function, however, is known in closed form due to Lieb; and since then surely also including the leading finite-size corrections. (See also the rigorous work of Duminil-Copin and collaborators on related models using the Bethe ansatz.)
Since the model is invariant under (discrete) translations in either direction there is no phase separation, and there are no limit shapes.

Edit. Let me outline how the problem of $x$-enumerating toroidal alternating-sign matrices, which I will abbreviate to tASM, is turned to that of solving coupled algebraic equations.
Step 1. Turn the problem into linear algebra. In physics this is called the transfer-matrix method. I also think that this is a computationally quite efficient way to find the number of tASM for sufficiently low $n,k$. Without loss of generality assume that $n\geq k$.
This goes as follows. Let me use the bijection that you mentioned between tASM and Eulerian (divergence-free) orientations of the square lattice with toroidal boundary conditions; the $\pm1$ correspond to the source/sink orientations around a vertex.
Consider the partition function (moment-generating function)
$$ Z = \sum_{C \,:\, \text{config}} \text{weight}(C) , \quad \text{weight}(C) = \prod_{v \, : \, \text{vertex}} \text{weight}(v) $$
for your favourite choice of vertex weights, e.g. $\text{weight}(v) = 1$ if orientation around $v$ is Eulerian, $=0$ else; or e.g. a refinement that counts the number of signs.
$$ Z = \text{tr}[t^k] $$
where $t$ is a $2^n \times 2^n$ matrix constructed as follows. Let $V = \mathbb{C} \, e_\uparrow \, \oplus \, \mathbb{C} \, e_\downarrow$ be the complex vector space with basis vectors $e_s$ labelled by the the two possible orientations of a single edge, $s \in \{\uparrow,\downarrow\}$. Its $n$-fold tensor product $W = V^{\otimes n}$ has basis $\{ e_{s_1} \otimes \dots \otimes e_{s_n} \}_{s_1,\dots,s_n}$ indexed by all possible configurations on $n$ edges, which we will think of as the $N$ vertical edges connecting two adjacent rows of vertices of the lattice. The transfer matrix is the linear operator $T \in \mathrm{End} \, W$ is defined as
$$\langle e_{s'_1} \otimes \dots \otimes e_{s'_n}, t \, e_{s_1} \otimes \dots \otimes e_{s_n}\rangle = \sum_{C} \prod_{v \, : \, \text{vertex} \,\in\, \text{row}} \text{weight}(v) $$
where the sum is over all configurations of orientations for the horizontal edges connecting the row of vertices, with periodic boundary conditions in the horizontal direction. As my description indicates we think of this as follows. Fix a row of $n$ vertices in the lattice. The matrix entries defined above give the weight (probability up to normalisation) where the configuration $e_{s_1} \otimes \dots \otimes e_{s_n}$ on the row of vertical edges directly below our row of vertices is transferred to the configuration $e_{s'_1} \otimes \dots \otimes e_{s'_n}$ on the row of vertical edges directly above it.
Define the $\mathfrak{sl}_2$-weight (in the sense of representation, not probability, theory) of a configuration on a row of vertical edges to be the number of $1\leq i\leq n$ such that $s_i = \,\downarrow$. Then, as the orientation is Eulerian and using the horizontal periodicity, $t$ preserves the $\mathfrak{sl}_2$-weight. The nonzero entries of $t$ precisely correspond to transfering a configuration with a given $\mathfrak{sl}_2$-weight to any other configuration with that same $\mathfrak{sl}_2$-weight. Moreover, the above sum defining the matrix entry of $t$ contains two terms if $s_i = s'_i$ for all $i$ (namely: all horizontal edges have the same orientation) and one term otherwise.
Next, $t^2$ transfers a configuration through two rows of the lattice (composition corresponds to summing over all possible configurations on the intermediate vertical edges). Likewise, the matrix entries of $t^k$ gives the (probabilistic) weight of all possible configurations on an $n\times k$ lattice with prescribed configurations on the vertical edges at the bottom and the top, and periodic boundary conditions in the horizontal direction. Finally, its trace (over $W$) amounts to imposing periodic boundary conditions in the vertical direction as well to get $Z$.
(In practice one defines $t$ in turn by a similar slicing of the row, thinking of it as consisting of a successing of operators --- acting e.g. from left to right --- encoding the weights of a single vertex. All of this is usually done via graphical notation, see e.g. the left column on p5 of the paper [arXiv:1702.05474] that I wrote with Rick Keesman, take `staggering' therein to be zero.)
Step 2. Compute the eigenvalues of $t$. Consider the case where the vertex weights are invariant under reversing the orientations on all edges. Then there are three vertex weights $a,b,c$ to specify. Overall rescalings only give an overall factor for $Z$ so are irrelevant. This leaves a two-parameter family of six-vertex models. It is known that at fixed value of the combination
$$ \Delta := \frac{a^2 + b^2 - c^2}{2ab}$$
the transfer matrices commute with each other at different values of the remaining parameter. This one-parameter family of commuting transfer matrices can be simultaneously diagonalised. By $\mathfrak{sl}_2$-weight conservation this may be done per fixed $\mathfrak{sl}_2$-weight. This is what can be achieved, as I mentioned above, via an ansatz for the form of the eigenvectors that is due to Bethe. If the $\mathfrak{sl}_2$-weight is $m$ then this involves $m$ free parameters that have to solve a set of coupled algebraic equations to guarantee that they are actually eigenvectors. Depending on the (probabilistic) weights the dominant eigenvalue, most relevant when $k$ grows, occurs either for $\mathfrak{sl}_2$-weight $m=0$ ('ferromagnetic') or instead at the 'equator' $m=n/2$ (if $n$ is even, 'antiferromagnetic'). In particular one can compute the dominant eigenvalue (via Perron--Frobenius) and its large-$n$ asymptotics. For the details see e.g. the recent preprint [arXiv:2012.11675] by Duminil-Copin et al.
