Geostatistical seismic inversion methods use stochastic sequential simulation as the model generation and perturbation technique. These stochastic simulation methods use a global variogram model to express the expected spatial continuity pattern of the subsurface elastic properties of interest. The conditioning to a single variogram model is not suitable for complex and non-stationary geological environments, resulting in poor inverted models unable to reproduce non-stationary features such as channels, folds, and faults. The proposed method uses a stochastic sequential simulation and co-simulation method able to cope with spatially varying information using local and independent variogram models. The information about the dip, azimuth, and ranges of the local variogram model is inferred directly from the observed data. First, local dip and azimuth structural volumes are computed from seismic attribute analysis. Then, local variogram models are fitted along the directions estimated from the previous step. This information is used as steering data during the inversion, acting as proxy of the true subsurface geological complexities. Application examples in synthetic and real datasets with complex geometries show the impact of using local anisotropy models in both the reproduction of the original seismic data and the reliability of the inverted models. The resulting inverted models show enhanced consistency, where small-to-large scale discontinuities and complex geometries are better reproduced, allowing reducing the spatial uncertainty associated with the subsurface properties. This work represents a step forward in integrating geological consistency into geostatistical seismic inversion, surpassing the limitation of using a single variogram model to reproduce complex geological patterns. © 2020, Springer Nature Switzerland AG.