The physical property models obtained from geophysical inversions can be converted to a 3D quasi-geology model via a process termed geology differentiation. Recent works show that geology differentiation can help maximize the value of information contained in geophysical data. However, it remains largely unexplored as to how to quantify the uncertainties of a 3D quasi-geology model. We approach this problem by using a recently developed mixed Lp norm regularization and a priori physical property measurements. We use mixed $L_p$ norm joint inversion to construct a large sequence of physical property models based on the Gzz component of the airborne gravity gradient and magnetic measurements. The available physical property measurements are used to determine which physical property models to accept. We then construct a sequence of 3D quasi-geology models by performing geology differentiation for all the accepted models, which allows us to compute the probabilities of our geology differentiation results. We apply our approach to a set of field data collected over the Decorah area located in northeast Iowa. We successfully quantify the uncertainties of the spatial extents for the identified geological units and compute probabilities of geologic units at any location in our study area. The proposed workflow has broad implications for 3D geological model building based on multiple geophysical and/or rock sample measurements.