Real-world signals typically span across multiple dimensions, that is, they naturally reside on multi-way data structures referred to as tensors. In contrast to standard ``flat-view'' multivariate matrix models which are agnostic to data structure and only describe linear pairwise relationships, we introduce the tensor-valued Gaussian distribution which caters for multilinear interactions -- the linear relationship between fibers -- which is reflected by the Kronecker separable structure of the mean and covariance. By virtue of the statistical identifiability of the proposed distribution formulation, whereby different parameter values strictly generate different probability distributions, it is shown that the corresponding likelihood function can be maximised analytically to yield the maximum likelihood estimator. For rigour, the statistical consistency of the estimator is also demonstrated through numerical simulations. The probabilistic framework is then generalised to describe the joint distribution of multiple tensor-valued random variables, whereby the associated mean and covariance exhibit a Khatri-Rao separable structure. The proposed models are shown to serve as a natural basis for gridded atmospheric climate modelling.