Let's enlarge the set of variables and consider more generally, for $a\in\mathbb{N}^n$ and $b\in\mathbb{N}^m$, the number of $n\times m$ matrices with non-negative integer coefficients whose $i$-th row sums to $a_i$ and whose $j$-th column sums to $b_j$, for $1\le i\le n$ and $1\le j\le m$: $$\mu(a,b):= \Big| \big\{v\in\mathbb{N}^{n\times m}\ :\ \sum_{1\le j\le m}v_{ij}=a_i\ , \sum_{1\le i\le n}v_{ij}=b_j\big\}\Big|\\ .$$ Then, for $x:=(x_1,\dots,x_n)$ and $y:=(y_1,\dots,y_m)$ $$\sum_{a\in\mathbb{N}^n\atop b\in\mathbb{N}^m}\mu(a,b)x^ay^b=\sum_{v\in\mathbb{N}^{n\times m} } \prod_{1\le i\le n \atop 1\le j\le m} x_i ^{\sum_jv_{ij}} y_j ^{\sum_iv_{ij}}=\prod_{1\le i\le n \atop 1\le j\le m}(1-x_iy_j)^{-1}\\ . $$ Thus $\mu$ is the discrete convolution of certain $nm$ log-concave functions $\delta _{ij}:\mathbb{N}^{n+m}\to[0,\infty)$ , namely the coefficients sequences of $(1-x_iy_j)^{-1}$. A discrete convolution of log-concave discrete functions $\mathbb{N}^p\to [0,\infty)$ is log-concave [edit: this has to be checked!].
In particular, for $n=m$ and $u=(1,\dots,1)$ the sequence $|M(n,k)|=\mu(ku,ku)$ is log-concave w.r.to $k\in\mathbb{N}\\ .$

Let's enlarge the set of variables and consider more generally, for $a\in\mathbb{N}^n$ and $b\in\mathbb{N}^m$, the number of $n\times m$ matrices with non-negative integer coefficients whose $i$-th row sums to $a_i$ and whose $j$-th column sums to $b_j$, for $1\le i\le n$ and $1\le j\le m$: $$\mu(a,b):= \Big| \big\{v\in\mathbb{N}^{n\times m}\ :\ \sum_{1\le j\le m}v_{ij}=a_i\ , \sum_{1\le i\le n}v_{ij}=b_j\big\}\Big|\, .$$ Then, for $x:=(x_1,\dots,x_n)$ and $y:=(y_1,\dots,y_m)$ $$\sum_{a\in\mathbb{N}^n\atop b\in\mathbb{N}^m}\mu(a,b)x^ay^b=\sum_{v\in\mathbb{N}^{n\times m} } \prod_{1\le i\le n \atop 1\le j\le m} x_i ^{\sum_jv_{ij}} y_j ^{\sum_iv_{ij}}=\prod_{1\le i\le n \atop 1\le j\le m}(1-x_iy_j)^{-1}\, . $$ Thus $\mu$ is the discrete convolution of certain $nm$ log-concave functions $\delta _{ij}:\mathbb{N}^{n+m}\to[0,\infty)$ , namely the coefficients sequences of $(1-x_iy_j)^{-1}$. A discrete convolution of log-concave discrete functions $\mathbb{N}^p\to [0,\infty)$ is log-concave [edit: this has to be checked!].
In particular, for $n=m$ and $u=(1,\dots,1)$ the sequence $|M(n,k)|=\mu(ku,ku)$ is log-concave w.r.to $k\in\mathbb{N}\, .$

Let's enlarge the set of variables and consider more generally, for $a\in\mathbb{N}^n$ and $b\in\mathbb{N}^m$, the number of $n\times m$ matrices with non-negative integer coefficients whose $i$-th row sums to $a_i$ and whose $j$-th column sums to $b_j$, for $1\le i\le n$ and $1\le j\le m$: $$\mu(a,b):= \Big| \big\{v\in\mathbb{N}^{n\times m}\ :\ \sum_{1\le j\le m}v_{ij}=a_i\ , \sum_{1\le i\le n}v_{ij}=b_j\big\}\Big|\\ .$$ Then, for $x:=(x_1,\dots,x_n)$ and $y:=(y_1,\dots,y_m)$ $$\sum_{a\in\mathbb{N}^n\atop b\in\mathbb{N}^m}\mu(a,b)x^ay^b=\sum_{v\in\mathbb{N}^{n\times m} } \prod_{1\le i\le n \atop 1\le j\le m} x_i ^{\sum_jv_{ij}} y_j ^{\sum_iv_{ij}}=\prod_{1\le i\le n \atop 1\le j\le m}(1-x_iy_j)^{-1}\\ . $$ Thus $\mu$ is the discrete convolution of certain $nm$ log-concave functions $\delta _{ij}:\mathbb{N}^{n+m}\to[0,\infty)$ , namely the coefficients sequences of $(1-x_iy_j)^{-1}$. A discrete convolution of log-concave discrete functions $\mathbb{N}^p\to [0,\infty)$ is log-concave. In particular, for $n=m$ and $u=(1,\dots,1)$ the sequence $|M(n,k)|=\mu(ku,ku)$ is log-concave w.r.to $k\in\mathbb{N}\\ .$

Let's enlarge the set of variables and consider more generally, for $a\in\mathbb{N}^n$ and $b\in\mathbb{N}^m$, the number of $n\times m$ matrices with non-negative integer coefficients whose $i$-th row sums to $a_i$ and whose $j$-th column sums to $b_j$, for $1\le i\le n$ and $1\le j\le m$: $$\mu(a,b):= \Big| \big\{v\in\mathbb{N}^{n\times m}\ :\ \sum_{1\le j\le m}v_{ij}=a_i\ , \sum_{1\le i\le n}v_{ij}=b_j\big\}\Big|\\ .$$ Then, for $x:=(x_1,\dots,x_n)$ and $y:=(y_1,\dots,y_m)$ $$\sum_{a\in\mathbb{N}^n\atop b\in\mathbb{N}^m}\mu(a,b)x^ay^b=\sum_{v\in\mathbb{N}^{n\times m} } \prod_{1\le i\le n \atop 1\le j\le m} x_i ^{\sum_jv_{ij}} y_j ^{\sum_iv_{ij}}=\prod_{1\le i\le n \atop 1\le j\le m}(1-x_iy_j)^{-1}\\ . $$ Thus $\mu$ is the discrete convolution of certain $nm$ log-concave functions $\delta _{ij}:\mathbb{N}^{n+m}\to[0,\infty)$ , namely the coefficients sequences of $(1-x_iy_j)^{-1}$. A discrete convolution of log-concave discrete functions $\mathbb{N}^p\to [0,\infty)$ is log-concave [edit: this has to be checked!]. 
In particular, for $n=m$ and $u=(1,\dots,1)$ the sequence $|M(n,k)|=\mu(ku,ku)$ is log-concave w.r.to $k\in\mathbb{N}\\ .$