Estimating covariance parameters for multivariate spatial Gaussian random fields is computationally challenging, as the number of parameters grows rapidly with the number of variables, and likelihood evaluation requires operations of order O((np)^3). In many applications, however, not all cross-dependencies between variables are relevant, suggesting that sparse covariance structures may be both statistically advantageous and practically necessary. We propose a LASSO-penalized estimation framework that induces sparsity in the Cholesky factor of the multivariate Matérn correlation matrix, enabling automatic identification of uncorrelated variable pairs while preserving positive semidefiniteness. Estimation is carried out via a projected block coordinate descent algorithm that decomposes the optimization into tractable subproblems, with constraints enforced at each iteration through appropriate projections. Regularization parameter selection is discussed for both the likelihood and composite likelihood approaches. We conduct a simulation study demonstrating the ability of the method to recover sparse correlation structures and reduce estimation error relative to unpenalized approaches. We illustrate our procedure through an application to a geochemical dataset with p = 36 variables and n = 3998 spatial locations, showing the practical impact of the method and making spatial prediction feasible in a setting where standard approaches fail entirely.
