<p>In the Earth system sciences, identifying physically meaningful patterns from dependent, high-dimensional, tensor-valued data is important for understanding variability in the Earth system. Orthonormal matrices are fundamental for statistical modeling of dependent data due to their role in matrix and tensor decomposition methods. While significant effort has been devoted to statistical inference for dependent matrix-valued data, these developments have not yet been extended to tensor-valued data. Here, we propose a general framework for dependent tensor-valued generalized linear mixed models (GLMMs) with a non-separable and non-stationary structured random effect. The random effect is parameterized using a tensor decomposition consisting of structured random orthonormal matrices, whose structures and associated coefficients are inferred during model estimation. Fully Bayesian inference enables joint estimation of all components of the dependent random effect and fixed effect within a unified statistical model. We validate our method with two synthetic examples and apply it to two real-world applications within atmospheric science, highlighting its capacity to uncover structured modes of variability across space, time, and vertical level. More broadly, our method provides a path forward for identifying structured patterns and their relative importance with Bayesian uncertainty quantification for arbitrary dependent, exponential family-distributed tensor-valued data.</p>