Random Domain Decomposition (RDD) provides a powerful tool for quantifying uncertainty in flow simulations, when both the geologic makeup of a porous medium and its hydraulic parameters are under-specified by data. Its prior applications dealt with flows in porous media whose internal compositions are amenable to simple parameterizations. This study provides a means for probabilistic reconstruction of boundaries between geologic facies. We apply our general approach to multiple data sets to reconstruct highly permeable zones within an aquitard in the Bologna (Italy) aquifer system and demonstrate how it can be used in conjunction with RDD.