Characterization of the biomaterial flow through porous bone is crucial for the success of the bone augmentation process in vertebroplasty. The biofluid, biomaterial, and local morphological bone characteristics determine the final shape of the filling, which is important both for the post-treatment mechanical loading and the risk of intraoperative extraosseous leakage. We have developed a computational model that describes the flow of biomaterials in porous bone structures by considering the material porosity, the region-dependent intrinsic permeability of the porous structure, the rheological properties of the biomaterial, and the boundary conditions of the filling process. To simulate the process of the substitution of a biofluid (bone marrow) by a biomaterial (bone cement), we developed a hybrid formulation to describe the evolution of the fluid boundary and properties and coupled it to a modified version of Darcy’s law. The apparent rheological properties are derived from a fluid-fluid interface tracking algorithm and a mixed boundary representation. The region- specific intrinsic permeability of the bone is governed by an empirical relationship resulting from a fitting process of experimental data. In a first step, we verified the model by studying the displacement process in spherical domains, where the spreading pattern is known in advance. The mixed boundary model demonstrated, as expected, that the determinants of the spreading pattern are the local intrinsic permeability of the porous matrix and the ratio of the viscosity of the fluids that are contributing to the displacement process. The simulations also illustrate the sensitivity of the mixed boundary representation to anisotropic permeability, which is related to the directional dependent microstructural properties of the porous medium. Furthermore, we compared the nonlinear finite element model to different published experimental studies and found a moderate to good agreement (R2=0.9895 for a one-dimensional bone core infiltration test and a 10.94–16.92% relative error for a three-dimensional spreading pattern study, respectively) between computational and experimental results.