A new method is proposed for decomposing covariance matrices that appear in the parameter estimation phase of all multivariate stochastic models in hydrology. This method applies not only to positive definite covariance matrices (as do the typical methods of the literature) but to indefinite matrices, too, that often appear in stochastic hydrology. It is also appropriate for preserving the skewness coefficients of the model variables as it accounts for the resulting coefficients of skewness of the auxiliary (noise) variables used by the stochastic model, given that the latter coefficients are controlled by the decomposed matrix. The method is formulated in an optimization framework with the objective function being composed of three components aiming at (1) complete preservation of the variances of variables, (2) optimal approximation of the covariances of variables, in the case that complete preservation is not feasible due to inconsistent (i.e., not positive definite) structure of the covariance matrix, and (3) preservation of the skewness coefficients of the model variables by keeping the skewness of the auxiliary variables as low as possible. Analytical expressions of the derivatives of this objective function are derived, which allow the development of an effective nonlinear optimization algorithm using the steepest descent or the conjugate gradient methods. The method is illustrated and explored through a real world application, which indicates a very satisfactory performance of the method. ¿ 1999 American Geophysical Union