Unconventional geomaterials often exhibit multi-modal pore size distribution. We have developed a comprehensive framework for porous media exhibiting multiple porosity scales that are saturated with one or two types of fluids using mixture theory. Both the governing equations and constitutive laws have been clearly derived and identified, respectively. The effective stress 𝝈′ emerged from the energy balance equation is adoptable for both elastic and elastoplastic deformations, in which pore fractions and saturations play a central role. The proposed model is general in a sense that it works for both uncoupled simulation and coupled simulation. The field equations for uncoupled flow simulation are solved using the Laplace transform and numerical Laplace inversion methods. By visualizing the dimensionless results, we can gain a quantitative insight of the different stages in the depletion process of a naturally fractured reservoir. For coupled flow and geomechanics simulation, a strip load problem and a two-phase flow in a deformable 3D reservoir problem illustrate the impacts of plasticity, multiple porosity, inter-porosity exchange, and capillary pressure on the system response.