This paper proposes a methodology to construct linear coregionalization models with as many nested structures as desired, together with as few orthogonal factors as possible. The construction rests on the decomposition of the model coregionalization matrices into pairwise commuting matrices, followed by a factorization by principal component analysis. The proposed approach is illustrated through a case study in mineral resources evaluation and compared to the traditional fitting procedure, obtaining an equally good fit of the direct and cross variograms but with significantly less factors.