  Multivariate spatial modeling is key to understanding the behavior of
materials downstream in a mining operation. The ore recovery depends on the
mineralogical composition, which needs to be properly captured by the model to
allow for good predictions. Multivariate modeling must also capture the
behavior of tailings and waste materials to understand the environmental risks
involved in their disposal. However, multivariate spatial modeling is
challenging when the variables show complex relationships, such as non-linear
correlation, heteroscedastic behavior, or spatial trends. This contribution
proposes a novel methodology for general multivariate contexts, with the idea
of disaggregating the global non-linear behavior among variables into the
spatial domain in a piece-wise linear fashion. We demonstrate that the complex
multivariate behavior can be reproduced by looking at local correlations
between variables at sample locations, inferred from a local neighborhood, and
interpolating these local linear dependencies by using a non-stationary version
of the Linear Model of Coregionalization. This mixture of locally varying
linear correlations is combined to reproduce the global complex behavior seen
in the multivariate distribution. The main challenge is to solve appropriately
the interpolation of the known correlation matrices over the domain, as these
local correlations defined at sample locations can be endowed with a manifold
structure, on which the Euclidean distance is not a suitable metric for
interpolation of such correlations. This is addressed by using tools from
Riemannian geometry: correlation matrices are interpolated using a weighted
Fr\'echet mean of the correlations inferred at sample locations. An application
of the procedure is shown in a real case study with good results in terms of
accuracy and reproduction of the reference multivariate distributions and
semi-variograms.
