As said in the comments above, it's easier to think of this as a word length problem on $\mathfrak{S}_S$, the symmetric group of $S$, in terms of certain transpositions (which we will also call "changes").

For $k = 2$, the formula is $2 N_1 N_2 - N_1 - N_2$. More generally:

Lemma: There is an element $\sigma \in \mathfrak{S}_S$ with minimal length exactly $|S| \sum_i (1 - \frac{1}{N_i})$ in terms of changes.

Let $\sigma \in \mathfrak{S}_S$ be the element that takes each $(a_1, a_2, ..., a_k)$ to $(a_1 + 1 (mod N_1), a_2 + 1 (mod N_2), ..., a_k + 1 (mod N_k))$. The minimum length expression of $\sigma$ in terms of changes has length $|S| \sum_i (1 - \frac{1}{N_i})$. The proof is similar to the proof that the long cycle in $\mathfrak{S}_n$ has minimum length $n - 1$ in terms of transpositions: each change can only increase one coordinate of the output for one input; the expression above is the number of increased coordinates in $\sigma$, and as such the minimal length is at least the above expression.

In more detail: For any $\sigma \in \mathfrak{S}_S$, let $m(\sigma) = \sum_{s \in S} |\{1 \leq i \leq k|\sigma(s)_i > s_i\}|$.

Claim: for any $\sigma \in \mathfrak{S}_S$ and change $\tau$, $m(\tau \sigma) \leq m(\sigma) + 1$.

Proof: The change only affects the sum for two elements $s, s'$ in $S$. Furthermore, it only affects one coordinate. Finally, for one of $s, s'$, that coordinate is decreased. Therefore, $m$ can increase by at most $1$.

Induction then shows that the minimum length of any $\sigma \in \mathfrak{S}_S$ is at least $m(\sigma)$. It's not hard to see that for our chosen $\sigma$, $m(\sigma) = |S| \sum_i (1 - \frac{1}{N_i})$, and we are done. $\square$

As a side note, this bound is sharp, that is, there is a way to express this $\sigma$ with that length: we can focus on one coordinate at a time. For each "axis" in the direction of some coordinate $i$, correcting that coordinate is equivalent to a long cycle on the $N_i$ elements, which therefore takes $N_i - 1$ changes. We therefore need $\frac{|S|}{N_i} (N_i - 1) = |S| (1 - \frac{1}{N_i})$ changes to correct the $i$th coordinate, which we can then leave alone. That totals to the above expression.

Lemma: For $k = 2$, any element $\sigma \in \mathfrak{S}_S$ can be written using at most $2 N_1 N_2 - N_1 - N_2$ changes.

Proof: When $N_1 = 1$, this says that any element of $\mathfrak{S}_{N_2}$ can be written using at most $N_2 - 1$ transpositions, which is already known to be true. As such, we have our base case.

The essential idea is that we consider the elements of $S$ as boxes in a rectangle; we start with the boxes lexicographically ordered, and we want to end up in some arbitrary reordering; we're allowed "moves" that switch two boxes either in the same row or column. The proof is done by correctly "filling" the last row, and bounding the number of moves we take to do so. We therefore only care about the boxes that will eventually go into that last row. There are $N_1$ of these.

We start by making sure there is one of these boxes in each column. We do this by repeating the following step: if there is an empty column, move the box that is supposed to be in that column to that column. This can only be done finitely many times (and we will later discuss how many times).

Once that is done, each column has exactly one box, by pigeonhole principle. Move every box to the last row. This takes at most $N_1$ changes.

Finally, as all boxes are in the last row, they can be moved to the correct position through changes that switch last-row boxes.

Claim: This takes at most $2 N_1 - 1$ changes.

Proof: First, we note that the first step doesn't affect at least one element. That is because at least one element is moved into the last empty column; as it doesn't leave an empty column, there must be some element there that wasn't moved.

The set of columns that were never empty is then also the set of columns affected by the third step. This set is nonempty, so (by the usual argument for word length with respect to arbitrary transpositions), the maximum number of transpositions in the third step is at most the size of the set minus 1.

As such, the total number changes in the first and third step is at most the number of columns that were ever empty plus the number of columns that were never empty minus 1, for a total of $N_1 - 1$. The second step takes at most $N_1$ steps, so the total number of steps is at most $2 N_1 - 1$.

Combining this with the base case, the induction is clear, and we are done. $\square$

I think that for $k > 3$, the above argument can be generalized, with enough care, to show that the lower bound in the first lemma is sharp, but don't immediately see it. The proof here essentially relied on the fact that we can choose any order in which to fix the elements of the last row, which is no longer true in the higher-dimensional case. I think a general proof has a good chance of working, using the facts that:

1) $M(N_1, N_2, ..., N_{k - 1}, 1) = M(N_1, N_2, ..., N_{k-1})$, and

2) $M(N_1, N_2, ..., N_k + 1) - M(N_1, N_2, ..., N_k) = N_1 N_2 ... N_{k - 1} + M(N_1, N_2, ..., N_{k-1})$

The first term should represent the number of moves in the $k$ direction that may be necessary to get everything to the last hyperplane (as in the second step in the above proof); the second term should represent the number of moves to rearrange everything in the other directions - though, as here, some of those rearrangements may have to happen outside of that plane.