Abstract-In this study, the problems of parameter estimation for systems with scarce measurements are examined. In this case, an original geostatistical methodology is applied to generate bathymetric surface models of a faulted aquifer system within the Jeffara Basin in southeastern Tunisia. The modelling workflow is based on i) the spatial analysis of the data configuration and ii) the conceptual stacking pattern. This allows provision of key concepts and geostatistical approaches to be undertaken during geomodelling procedures. In fact, two constraints have been integrated: i) inequality data provided by the end of the boreholes and ii) the faults that compartmentalize the aquifer system. First, kriging with inequality was used for depth estimation of the Turonian reservoir Top. The results are compared with those from classical kriging and evaluated through the estimation quality, the adopted assumptions and the method limitations. Validation analyses show that the model developed with the inequality data leads to a significant gain in mapping accuracy and geological realism. In the second step, the geostatistical approach was used to model the successive bounding surfaces of the aquifer system units. Consequently, the kriging constrained by inequality data could be applied in a variety of hydrogeological parameters interpolations to provide significantly better informative maps that are useful for hydrodynamic modelling.