Data assimilation methods are commonly used to predict petrophysical properties of deep saline aquifers for carbon dioxide sequestration studies. However, data assimilation is usually computationally challenging for large-scale geological models with a large number of geophysical observations. A computationally efficient method is developed for large-scale data assimilation problems. In the proposed method, model parameters and observations are reduced by randomized tensor decomposition, and the inversion is performed in lower-dimensional spaces. Tensors are multiway arrays that generalize matrices to multiple dimensions and therefore provide a natural way to represent three-dimensional geological models and geophysical measurements. To exploit the inherent multidimensional structure of tensors, high-order singular value decomposition performs orthogonal decomposition of the tensor in the high-order space instead of using a flattened matrix representation. However, for complex geological models, the dimension of tensors becomes very large, making the decomposition computationally unfeasible. To overcome this limitation, the dimension of the tensors is first reduced by a randomized linear algebra algorithm which guarantees that most of the information of the original tensors is preserved with high probability in the reduced tensor, and then high-order singular value decomposition is applied on the smaller tensors with relatively small computational cost. By using randomized tensor decomposition for model and data re-parameterization, data assimilation is efficiently performed in the low-dimensional model and data space. The proposed method is demonstrated and validated on the three-dimensional data set synthesized based on the Johansen formation, a potential carbon dioxide storage site offshore in Norway.