The purpose of seismic reservoir characterization is to predict the spatial distribution of the subsurface rock properties from a set of direct and indirect measurements. Commonly, rock property volumes are obtained in a two-steps approach. First, the elastic properties are inverted from seismic reflection data and then used to compute rock properties volumes by applying a calibrated rock physics models. Such an approach is not only time consuming but may also lead to biased results as the uncertainties related to seismic inversion may not be propagated through the entire geomodeling workflow. Here, we propose an iterative geostatistical shale rock physics seismic AVA inversion method to invert seismic reflection data directly for shale rock properties. The workflow consists of three main steps starting from shale properties model generation of brittleness, total organic carbon (TOC), and porosity using stochastic sequential simulation and cosimulation and calculation of a volume of shale. In the following step, elastic property volumes are calculated based on a calibrated shale rock physics model using the self-consistent approximation. The elastic models are then used for the calculation of the synthetic seismic data. In the final step, the misfit between synthetic and real data is calculated and used as part of the stochastic update of the model parameter space. The whole process is repeated until a minimum misfit between the observed and synthetic seismic is achieved. The proposed method is successfully tested on an onshore Lower Paleozoic shale gas reservoir in Northern Poland, and the predicted shale rock properties agree with those observed at a blind-well location.