You have a matrix containing $T$ observations of each of $K$ random variables \begin{align} U = \begin{bmatrix} u_{11} & \dots & u_{1T} \\ \vdots & \ddots & \vdots \\ u_{K1} & \dots & u_{KT} \end{bmatrix}. \end{align} The covariance matrix of the $u_{kt}$'s is given by $\Sigma_u := E(u_t u_t')$ for $t = 1, \dots, T$ and where $ u_t = (u_{1t}, \dots, u_{Kt})'$ is a column vector by convention.
Now, we introduce the $vec$ operator which stacks columns of a matrix one over the next, i.e. \begin{align} vec(U) = \begin{bmatrix} u_{11} \\ \vdots \\ u_{1T} \\ \vdots \\ u_{K1} \\ \vdots \\ u_{KT} \end{bmatrix} = \begin{bmatrix} u_{1}' \\ \vdots \\ u_{K}' \end{bmatrix} \end{align} again using as convention that $u_k = ( u_{k1}, \dots, u_{kT} )'$ are column vectors for $k = 1, \dots, K$. It turns out that the covariance matrix of $vec(U)$ is \begin{align} E\left( vec(U) vec(U)' \right) &= I_T \otimes \Sigma_u \end{align} which is to say, a block diagonal $TK \times TK$ matrix with entries $E( u_t u_t' )$. While I can see the intuition, I don't know how to formally establish this result. Any help would be appreciated.