A continuum mechanical porous media model for vertebroplasty: Numerical simulations and experimental validation

The outcome of vertebroplasty is hard to predict due to its dependence on complex factors like bone cement and marrow rheologies. Cement leakage could occur if the procedure is done incorrectly, potentially causing adverse complications. A reliable simulation could predict the patient-specific outcome preoperatively and avoid the risk of cement leakage. Therefore, the aim of this work was to introduce a computationally feasible and experimentally validated model for simulating vertebroplasty. The developed model is a multiphase continuum-mechanical macro-scale model based on the Theory of Porous Media. The related governing equations were discretized using a combined finite element–finite volume approach by the so-called Box discretization. Three different rheological upscaling methods were used to compare and determine the most suitable approach for this application. For validation, a benchmark experiment was set up and simulated using the model. The influence of bone marrow and parameters like permeability, porosity, etc., was investigated to study the effect of varying conditions on vertebroplasty. The presented model could realistically simulate the injection of bone cement in porous materials when used with the correct rheological upscaling models, of which the semi-analytical averaging of the viscosity gave the best results. The marrow viscosity is identified as the crucial reference to categorize bone cements as ‘high- ’or ‘low-’ viscosity in the context of vertebroplasty. It is confirmed that a cement with higher viscosity than the marrow ensures stable development of the injection and a proper cement interdigitation inside the vertebra.


Introduction
Annually, more than 1.4 million clinical vertebral fractures occur worldwide (Johnell and Kanis 2006). Vertebral fractures lead to back pain, loss of height, immobility, reduced pulmonary function, and increased mortality. Vertebroplasty is a common medical procedure for treating and preventing vertebral fractures. In this procedure, a so-called 'bone cement' is injected into the porous interior of the vertebra, where it undergoes curing, providing additional structural strength to the vertebra (Jensen et al. 1997). Vertebroplasty provides quick pain relief in 80-90 % of the cases (McGraw et al. 2002). Generally, the practitioners use haptic feedback and X-ray imaging to guide the cement injection. However, this is difficult due to several factors, including the non-constant viscosity of the bone cement owing to its non-Newtonian nature and the curing process. The bone marrow in the vertebra is also a non-Newtonian fluid. Severe complications could occur if the cement leaks into the blood vessels or the spinal canal due to improper injection, e.g. pulmonary embolism or paralysis (Bernhard 2003;Ratliff et al. 2001). In this regard, a reliable simulation of the procedure could help practitioners predict the risks and determine the safe ranges for operating parameters like injection pressure, cement viscosity, curing time, flow rate, etc., specific to each case. The decrease in reliance on X-ray imaging would also reduce the patient's radiation exposure.
In the past, various approaches have been used for simulating vertebroplasty. Theoretical approaches, (e.g. Bohner et al. (2003)), usually have limited use-case because many idealizing assumptions are necessary. The earliest works using computational approaches were limited to small two-dimensional segments and assumed the porous structure to be a bundle of capillary tubes, e.g. in Beaudoin et al. (1991). With improvement in computational capabilities, more complicated models appeared, such as branching pipe networks (Lian et al. 2008), computational fluid dynamics (CFD) models Landgraf et al. 2015) or Lattice-Boltzmann models (Zeiser et al. 2008). The main drawback in these models is that a porescale flow resolution inside a complicated porous structure quickly becomes computationally expensive. Therefore, most pore-scale models simulate flow only over a small part or a representative segment of the vertebra. For a simulation regarding the entire vertebra, a macro-scale model is a more feasible option. In general, there is substantial research on the mathematical modelling of such multiphase porous media flow problems at the macro-scale; however, most of it is in the context of groundwater flow. It is possible to extend and adapt such models to vertebroplasty, but a major challenge is to upscale the viscosities of the non-Newtonian fluids for use in the macro-scale flow equations. Widmer Soyka et al. (2013) developed such a macro-scale model using a continuum length scale Reynolds number for upscaling the viscosity. However, they employed a front-tracking volume of fluid method that assumes a sharp cement-marrow interface (Widmer 2011), which may not always be the case. In this regard, the Theory of Porous Media (TPM) (Ehlers 2002(Ehlers , 2009) provides a solid framework for implementing a multiphase continuum-mechanical macro-scale model to simulate the fluid flow in a porous medium as well as the deformations in the porous medium, without needing to idealize the interface between the fluids. The application of TPM application extends to, for example, the biomechanics of soft tissues such as liver (Ricken and Lambers 2019), brain (Ehlers et al. 2022), and cartilage (Wang et al. 2018). TPM has previously also been used to develop a preliminary model for vertebroplasty (Bleiler et al. 2015).
The aim of this work was to develop a computationally feasible model based on TPM, and to extend it beyond previous studies (Bleiler et al. 2015) for better suitability of simulating vertebroplasty by including a more appropriate choice of model variables and constitutive equations for non-Newtonian rheology, investigating suitable rheology upscaling models, and numerically treating the governing equations by a combined finite element-finite volume approach using the Box discretization to address the different natures of the fluid flow and the solid deformation problem in a stable framework. Furthermore, the work aims to validate the model using a benchmark problem, which is investigated experimentally as well as simulated computationally, and carry out further simulations to determine the influence of various parameters, with specific regard to avoiding cement leakage.

Multiphase modelling based on the Theory of Porous Media
The Theory of Porous Media (TPM) provides a framework for modelling the mechanics of a multiphase continuum at the macro-scale (Ehlers 2002(Ehlers , 2009. It originates from the theory of mixtures (Bowen 1984), where the porescale properties of the phases are homogenized over a so-called representative elementary volume (REV), such that the macroscopic properties of this volume are representative of the entire domain. In our case, the overall aggregate contains three phases: the solid trabecular bone S , the bone cement C , and the bone marrow M , where the fluid phases = {C, M} are immiscible since the bone cement is hydrophobic. In the TPM, this is further extended by including the individual amount of each phase using the concept of volume fractions. This is shown in Fig. 1a. The volume fraction n = dv ∕dv , where dv is the volume of the phase dv and dv is the total volume of the REV, sum up to one, i.e. ∑ n = 1 . Thus, the medium does not have any void spaces. Similarly, the saturations s = n ∕n F of the fluid phases are introduced.
Under loading conditions, such a system undergoes a change in its state with time. The state at the current time t is called the current configuration, whereas that at the reference time t = 0 is the reference configuration. Each spatial point with a position vector x in the current configuration is occupied by all three phases (superimposed continua). Each of these three phases may originate from a different point with position vector X in the reference configuration, depending on their respective placement function , such that x = (X , t) . The material time derivative of each phase is defined as � x = d( (X , t))∕dt . The deformation of the bone is then defined in a Lagrangian setting as u = x − X , while the movement of the fluid phases is described in a 'modified' Eulerian setting by their seepage velocities, given as w = � x − � x S . The kinematics are shown in Fig. 1a.
These variables are governed by a set of balance equations, viz., the mass balance and the momentum balance equations for each of the phases. These are obtained after applying the assumptions that all phases are materially incompressible, no mass exchange occurs between the phases, the porous medium is fully saturated, and the constituents are non-polar (stress tensors are symmetric). Additionally, quasi-static processes are assumed, body forces are neglected, and all processes are considered isothermal. This yields the volume balance law for the phases as For the solid, this becomes simply where n S 0S is the initial volume fraction and F S = dx∕dX S is the deformation gradient. For fluids, the volume balance law in the modified Eulerian setting becomes Additionally, we have the overall momentum balance given as where = S E − pI . The solid extra stress S E is related to the deformation by constitutive relations, and p is the hydrostatic pore pressure. For the considered application, the capillary number is large ( Ca ≫ 10 −3 ), implying that the capillary forces are negligible compared to the viscous forces. Therefore, we neglect the pressure difference between the two fluids arising due to capillary effects and assume the same pore pressure p for both fluid phases. The resulting stress of the solid skeleton is obtained by S = S E − n S pI , and the fluid stress is given by = −n pI . Here, we neglect the fluid extra stresses E arising from viscous effects compared to the local interaction forces, as can be concluded from dimensional analysis (Markert 2007).
Furthermore, the Darcy flow equations are used to relate the seepage velocities of the fluid phases to the pressure gradient. These equations can be derived in the framework of the TPM by applying constitutive relations (4) div = 0, where K S is the intrinsic permeability of the porous medium, r is the relative permeability of the fluid phases arising due to the resistance from the other fluid, and R is the dynamic viscosity. For a detailed overview of the TPM, the reader is referred to the works of Ehlers (Ehlers 2002(Ehlers , 2009).

Numerical treatment
The governing equations are solved over the entire domain Ω of the porous medium using the primary variables solid displacement u S , pore pressure p, and bone marrow saturation s M . Since we neglected the capillary forces for our intended application, we used saturation as one of the primary variables, making the formulation independent of the capillary pressure-saturation relation, in contrast to one with two pore pressures as primary variables (like in Bleiler et al. (2015)). The volume balance for the solid (Eq. 2) is implemented directly such that it is satisfied at all points, while those for the fluids (Eq. 3) and the momentum balance (Eq. 4) are satisfied weakly, i.e. only their weighted integrals over the domain must be zero. This is written as where u S , p , and s M are test functions. Applying the Gaussian divergence theorem, the surface integrals for the Neumann boundary conditions appear via where n is the surface normal vector. Equation 10 is spatially discretized and solved using the Bubnov-Galerkin method by the finite element discretization. In Eqs. 11 and 12, the volume integral terms consist of a temporal evolution term, a flow term, and a contribution from the solid deformation, therefore requiring both temporal and spatial discretization. For the temporal discretization, we used the implicit Euler method. For flow problems, especially those hyperbolic by nature, it is known that discretization using the Bubnov-Galerkin finite element method leads to numerical instabilities (Zienkiewicz and Heinrich 1978). Therefore, we used the Box spatial discretization scheme and applied mass-lumping for the temporal evolution term. This yields a setting which behaves similar to the finite volume discretization. The flow equations were then solved using the fully upwind Galerkin scheme. More details about the Box discretization can be found in Huber and Helmig (2000). For the last remaining term, i.e. the contribution from the solid deformation, we again used the finite element discretization.
The system was discretized using eight-noded hexahedral elements with linear shape functions for all primary variables. The numerical implementation of the model was done on the monolithic solver framework PANDAS 1

Constitutive models
The governing balance relations need to be supplemented by constitutive relations to include the specific material behaviour of the constituents. For the solid part, i.e. the trabecular bone, the relationship between the stress and the strain is described using the linear Hookean law where S E is the solid extra stress, S is the strain tensor, Lamé and Lamé are the Lamé material parameters.
For the two fluids, constitutive relations are required for the change in viscosity with shear rate. We used the Carreau model for this purpose, given for a fluid as where 0 is the viscosity limit at zero shear rate, ∞ is the viscosity limit at very high shear rates, rh is the reciprocal of the shear rate where the behaviour transitions to non-Newtonian, n rh is the flow behaviour index, and ̇ is the shear rate. Note that the superscript ' ' stands for the fluid phase (C or M for cement or marrow respectively) and the subscript 'rh' stands for 'rheological parameter' to avoid confusion with similar symbols used in other places.
The viscosity discussed up to this point is a pore-scale quantity. Unlike for a Newtonian fluid, the pore-scale viscosity cannot be directly used in the macro-scale equations (Eqs. 5, 6) for a non-Newtonian fluid because the shear rates, and hence the viscosities vary with the geometries of the pores through which the fluid flows. The viscosity at the macro-scale must be a representative average of the pore-scale viscosities obtained through upscaling. Here, we compared three different upscaling models to evaluate their applicability to model vertebroplasty at the macro-scale. The first two, namely the Cannella model (Cannella and Seright 1988) and the Hirasaki and Pope model (Hirasaki and Pope 1974), are based on the capillary-tube-bundle model and semi-empirical by nature. In both approaches, the effective shear rate is obtained from the equation given as where ̇e ff is the effective upscaled shear rate, and ℂ is a semi-empirical constant, the value of which varies depending on the internal geometry of the porous medium. The Cannella model uses the value ℂ = 6.0 , while the Hirasaki and Pope model uses ℂ = 0.69 . The effective shear rate is then plugged in Eq. 14 to obtain the effective viscosity at the macro-scale. Alternatively, the third upscaling method under consideration was the average viscosity model by Eberhard et al. (2019), which is semi-analytical by nature in contrast to the previous two. This approach analytically computes the average viscosity of the fluid flowing through a representative pore channel of characteristic radius R char . This characteristic radius R char is determined from the pore size distribution of the porous medium. The details about the derivations of the equations and their analytical solution can be found in Eberhard et al. (2019).
Finally, constitutive relations are required for the relative permeabilities. Relative permeability is a factor applied to account for the reduction in permeability due to mutual hindrance caused by the fluid phases. Hence, the constitutive relation is usually a function of saturation. In this work, the Brooks-Corey model (Brooks and Corey 1964) was used, given as where bc is the uniformity parameter. A large value means that pore sizes in the porous medium do not strongly vary, while a small value represents non-uniformity in the pore sizes.

Validation with benchmark problem
A simple problem, based on the actual vertebroplasty procedure, was formulated as a benchmark to validate the numerical model. This is described with a simplified schematic in Fig. 2a. In the problem, 2 ml of bone cement is injected into a piece of aluminium foam using a surgical syringe and cannula at two different flow rates, viz., 0.1 ml/s and 0.4 ml/s. Aluminium foam was chosen since its internal structure is very similar to a vertebra (Loeffel et al. 2008).

Experimental setup
To carry out this benchmark experiment, a custom-made bone cement injector was designed and manufactured, as shown in Fig. 2b. The injector consisted of a carriage which was driven by a stepper motor using a ball screw with 5-mm feed per revolution. The stepper motor had a resolution of 1.8 • corresponding to 200 steps per revolution, amounting to a calculated carriage resolution of 0.025 mm. A 200 N load cell was mounted on the carriage to measure the forces applied to the plunger of the syringe. The bone cement of brand name Vertecem V+, purchased as a non-sterile batch (OSARTIS GmbH, Germany), was used for the injection. Preparing the bone cement requires mixing the polymer powder and the monomer liquid provided in the Vertecem V+ Cement Kit (DePuy Synthes). The components were directly mixed in a 10-ml syringe to easily transfer it to the 2-ml syringes of the kit for injection, using a method described later in Sect. 2.2.3. The time from the start of mixing to the start of injection was measured using a stopwatch, which was about 200-210 s. The bone cement was injected by the injector at a prescribed rate into the aluminium foam, as shown in Fig. 2c. The aluminium foam was cut into a cylinder with a height of 30 mm and a diameter of 40 mm, for a size similar to that of a vertebra. The aluminium foam was placed in a polymethyl methacrylate (PMMA) housing with a hole in the centre for inserting the cannula and

Numerical implementation
The flow into the aluminium foam was simulated using our numerical model, the geometry and boundary conditions for which are shown in Fig. 2d. The geometry and mesh were   (2008) Fig. 2a. The geometry was discretized by a mesh of 2540 hexahedral elements.

Material parameters
The parameters and their values used for the benchmark problem are summarized in Table 1. The aluminium foam was made of standard material Aluminium 6101, thus the material parameters were easily available from literature. Since the aluminium foam consisted of nearly uniformly sized pores, bc = 3 was taken as a reasonable estimate. The rheological parameters of the bone cement include the upper and lower viscosity limits C 0 and C ∞ , the relaxation time C rh for the transition to shear-thinning behaviour, and the flow behaviour index n C rh for the shear-thinning. The upper viscosity limit and its change with time were determined by rheological measurement. We used multiple methods of preparation of bone cement for rheological measurements to simulate real kit application. After multiple trials, the successful and reproducible method consisted of the following steps: 1. Pre-weighing 2.6 gs of the polymer powder in a 10-ml beaker. 2. Separately, 10 ml of the monomer liquid was prepared in a batch glass vial. 3. At time t = 0 s, 1.0 ml monomer liquid was dropped into the beaker with the polymer powder using a positive displacement pipette.
4. Immediately at that time (t = 0 s), a timer was started to count 20 s of gentle stirring using a polyetheretherketone stirring rod (counting 20 rotations). 5. After 20 s, part of the sample was gently dropped onto the rheometer bottom plate and the top plate was gently lowered down. The humidity control hood was used to avoid sample evaporation, and silicone oil was gently spread around the sample directly before each measurement.
Understanding the timing of polymerization is crucial; therefore, we evaluated the average time taken from the start of the mixing of cement components to the start of the rheological measurement with the above-established procedure, which was 118 ± 10 s (based on a total of 26 trials). The rheological measurements were performed on MCR302 rheometer from Anton Paar. Parallel plate geometry (PP) with a 25-mm diameter top plate and a 1.5 mm gap was used. Time sweep measurement using plate-plate with a total measurement time of 20 min was run in a rotational mode at a shear rate of 0.2 s −1 . Each point was recorded every 5 s. The results are shown in Fig. 3a. The viscosity at 3.5 min, i.e. the cement preparation time in the injection experiments, was used to obtain C 0 = 1930 Pa s. The actual value of the lower viscosity limit C ∞ is difficult to obtain due to its occurrence at shear rates too high to be able to measure reliably, although the viscosity has been observed to reduce by at least two orders of magnitude without plateauing (Krause et al. 1982;Lepoutre et al. 2019). Therefore, the value C ∞ = 1.93 Pa s was assumed such that a 1000times reduction in viscosity is allowed. The cement showed Fig. 3 (a) Measurements used for obtaining rheological parameters of bone cement: change in viscosity with time due to curing (left) and injection pressure required for cement flow through the injection equipment without the aluminium foam (right) (b) Pore-network extraction from micro-CT image of the aluminium foam a sharp increase in viscosity for about first 2 − 2.5 min after the start of mixing, after which the curing slowed down for the next 9-10 min, followed by rapid curing at the end. In clinical applications and our benchmark experiment, the injection occurs in the period from 3 to 12 min, where the viscosity increase was found to be linear at 0.32 Pa s/s. This increase would have a negligible effect on the viscosity during the period of injection for our experiments. Hence, the dependence on time was neglected.
To obtain the remaining two parameters, C rh and n C rh , the experimental setup as explained in Sect. 2.2.1 was used without the aluminium foam to obtain the injection force required for flow through the assembly made of the syringe, the nozzle, and the cannula, at flow rates 0.1 and 0.4 ml/s. The injection pressure was then obtained by dividing the force by the area of the syringe. The injection pressures are shown in Fig. 3b. The value immediately after the initial ramp is taken. These injection pressures and their respective flow rates were inserted in the analytical solution for the pressure required by a Carrreau fluid flowing through a tube at a given flow rate, as given in Sochi (2015). This gave a system of two equations, which could be solved to obtain the two remaining unknowns C rh = 1.38 and n C rh = 0.30 . All parameters are listed in Table 1. The remaining parameters, viz., porosity, permeability, and the characteristic radius R char (used in the average viscosity rheological model), are dependent on the pore geometry of the porous medium. A micro-CT scan of the aluminium foam was carried out to resolve the pore structure of the aluminium foam, from which the porosity could be computed. The micro-CT scanner used here was a VivaCT40 from SCANCO Medical AG at an isotropic resolution of 19μm using 70 kV, 114μA and 200 ms integration time. To determine the permeability and the characteristic radius, we used a pore-network model (Joekar-Niasar et al. 2010). The micro-CT image was used to extract a pore-network, as shown in Fig. 3c, using an open-source Python toolkit PoreSpy (Gostick et al. 2019). The permeability was obtained from this pore-network by simulating a Stokes flow through this pore-network using the open-source porous media flow solver DuMux (Koch et al. 2021). Note that the permeability obtained here was isotropic, i.e. nearly the same in all directions, hence we could use a scalar value instead of the tensor in Eqs. 5 and 6. Similarly, the characteristic radius R char was obtained from the volume-weighted average of the radii of the pores and the throats.

Clinically relevant simulations and further investigations
To investigate a more clinically relevant setting, we used the same geometry and boundary conditions as in the benchmark problem with Cannella model and 0.4 ml/s flow rate, while suitably replacing the values of material parameters.
Firstly, we replaced the mechanical properties of the aluminium foam with those of trabecular bone (Wu et al. 2018).
Since the mechanical properties of the trabeculae vary depending on factors like specimen, location, and condition, mean values found in the literature were used. Furthermore, instead of a previously empty porous medium, we did trials assuming the presence of bone marrow. In reality, bone marrow could consist of red bone marrow, which due to its similarity to blood, is non-Newtonian by nature, and yellow bone marrow, which has Newtonian rheology (Gurkan and Akkus 2008). The rheological behaviour could vary depending on the red-to-yellow bone marrow ratio, which is then dependent upon the age and health conditions of the patient. Given this uncertainty, we did the simulations over a range of viscosities, considering both non-Newtonian and Newtonian rheologies. For non-Newtonian rheologies, we assumed realistic values for the parameters based on the literature (Gurkan and Akkus 2008;Davis and Praveen 2006;Metzger et al. 2014;Jansen et al. 2015), while the viscosity limits were set such that M 0 ∕ M ∞ = 10 . The rest of the parameters used were the same as in the benchmark experiment. The parameters and their values are summarized in Table 1.
Apart from the above, further trials were carried out to investigate the effect of permeability, porosity, and the Brooks-Corey uniformity parameter bc . In reality, these parameters are interdependent, e.g. a higher porosity also leads to a higher permeability. However, here we wanted to investigate them in isolation to study their standalone effects. The chosen values for permeability and porosity were those obtained from literature for human vertebrae (Baroud et al. 2004;Daish et al. 2017;Ochia and Ching 2002). As far as the author is aware, no information is available in the literature for the Brooks-Corey uniformity parameter of human vertebrae. Therefore, a general but wide range of values from 0-10 was chosen.

Validation with benchmark problem
Looking at the cement distribution first in Fig. 4a, similar cement distribution patterns, i.e. fully saturated ( s C ≈ 1 ) with a sharp front, were obtained from the simulation and in the experiment of the benchmark problem. However, the injection forces, shown in Fig. 4b, differed from the experiments depending on which rheology upscaling model was used. The difference also did not remain constant with time. Table 2 shows the percentage error averaged over the time, barring the 10% at the start and the end to ignore any effect of the difference in climb/drop times. The average viscosity model gave results closest to the experiments. The Cannella model gave results with more error but in a similar range, whereas the Hirasaki model performed the worst with much higher magnitudes of force.

Clinically relevant simulations
From Fig. 5d, it was observed that except at the peripheries, the viscosities were the same everywhere within the fluids despite the non-Newtonian rheology. The cement viscosity after shear-thinning was 5.9 Pa s. Figure 5a shows a stark difference in the injection force development in cases with marrow viscosity higher than this value of 5.9 Pa s, compared to the other way round. For cases where marrow viscosity was higher, the force peaked near the start and declined afterwards, whereas when marrow viscosity was lower, the force showed a small gradual increase with time.
The difference in the injection force between the non-Newtonian and the Newtonian marrow cases was only in magnitude, with higher values observed in the case of Newtonian marrow since it does not undergo shear-thinning. Another difference caused due to the marrow viscosity is shown in Fig. 5c. Higher marrow viscosity caused the cement distribution to be diffused, as observed from low cement saturation values, and to spread farther, nearly touching the end of the porous medium. For lower marrow viscosity, the cement cloud was fully saturated and compact, similar to the benchmark experiment. Note that with both the aluminium foam and the trabecular bone as the solid skeleton, the negligible deformations in the order of 10 −10 m and 10 −9 m, respectively, were obtained. The results can be referred to in the dataset (Trivedi et al. (2022)). Figure 6 shows how the injection force at the end of injection varies with permeability, porosity, and the uniformity parameter bc . Only permeability showed a significant influence out of the three parameters. Porosity had nearly no influence. For bc , there was some increase in injection force at bc = 0.1 , but that is already an extreme and unrealistic case. However, since bc represents the pore size  uniformity in the porous medium, we also investigated its effect on the cement distribution. Therefore, we considered two cases: with lower marrow viscosity (10 −5 Pa s) and with higher marrow viscosity (10 3 Pa s) relative to the cement viscosity since we observed a difference in cement distribution for these cases earlier in Fig. 5c. Accordingly, we compared the injection force graphs and cement distributions for various values of bc for both these cases, as  Fig. 6b. We observed no significant difference in the injection forces due to bc as long as bc > 0.1 . However, the cement distribution did show dependence on bc for the case of higher marrow viscosity, although this was not the case when the marrow viscosity was lower than the cement. The exception here is again the extreme case of bc = 0.1.

Discussion
The matching cement patterns in the experiment and the simulations for the benchmark problem show that the presented model can correctly describe and replicate the flow boundary conditions from the experiments. In terms of injection force, the average viscosity model for upscaling the rheologies gave the best results, followed by the Cannella model. The error, although not negligible, is expected due to uncertainty caused by factors like the sensitivity of the cement curing to the mixing conditions and exposure to air, which are difficult to control given the constraints of the rheological measurement and the benchmark experiments. The Hirasaki and Pope model performed the worst among the three. In contrast, the study by Eberhard et al. (2019) found the results from Hirasaki and Pope and the average viscosity models similar and more accurate than the Cannella model. The reason for better accuracy of the Hirasaki and Pope model in their study could be that they used a 'granular' structure made of monodisperse spheres, which is also what Hirasaki and Pope (1974) used to arrive at the value of ℂ = 0.69 (in Eq. 15) in their model. Comparatively, the 'fibrous'-type structure of the aluminium foam is considerably different. Therefore, the Cannella and the average viscosity models are better suited than the Hirasaki and Pope model to upscale the viscosity for the case of aluminium foam, and therefore, also for vertebrae. The average viscosity model gave the most accurate results not only in the work of Eberhard et al. (2019), but also in ours, despite the difference in the porous materials used, hinting at its potential applicability to a relatively wide range of porous media. The reason for its wider applicability compared to the semi-empirical models could be because it requires an additional parameter (the characteristic radius R char ) derived from the pore geometry, whereas the other two models use fixed values for the constant ℂ . Interestingly, the injection force increased with time in the simulations even though the time-dependent increase in viscosity due to curing was neglected in the model. Therefore, the increase in the force must have to do with the nature of the cement distribution. Also, note that in both experiment and simulation, a four-times increase in flow rate requires less than double the injection force. The nonlinear dependence is due to the decrease in the cement viscosity by shear-thinning at higher flow rates.
The presence of bone marrow had a considerable impact on the results. Whether the marrow is non-Newtonian or Newtonian affects only the magnitude of the injection force, the lesser force being in the non-Newtonian case due to the shear-thinning viscosity. The injection force development and cement distribution are significantly affected by whether the cement viscosity is higher or lower than the marrow viscosity. There are three clear risks if the cement viscosity falls below that of the marrow: (i) unintuitive injection force development, (ii) improper cement filling, and (iii) additional dependence on the pore size uniformity. Firstly, the required injection force is highest at the start, followed by a substantial dip. The peak force needed at the start requires the practitioners to apply more effort early on, which could lead to sudden excess cement injection as the required force suddenly dips. Hence, such injection force behaviour is highly unintuitive for practitioners relying on haptic feedback. Secondly, the injected cement only partially fills the porous bone before spreading further. The marrow is not fully displaced because its higher viscosity makes it less movable than the cement. Therefore, the cement spreads in an unstable fingering pattern in a phenomenon known as 'viscous fingering'. Such a filling pattern causes poor interdigitation of the cement with the trabeculae. It also causes the cement to spread farther, potentially leaking outside the cortical shell before the required volume is injected. Finally, since the fingers develop depending on the pore geometries, an additional dependence of the cement distribution on the pore size uniformity parameter bc is introduced. This additional dependence adds another factor of uncertainty to the procedure outcome. On the other hand, these risks could be avoided if the cement viscosity stays higher than the marrow's, yielding advantages like (i) a gradually rising injection force, which is more stable and intuitive for practitioners; (ii) the cement completely filling the porous bone, providing better interdigitation and leading to better mechanical strengthening; (iii) no strong dependence on the pore size uniformity.
The advantages of high-viscosity cement have been emphasized in earlier works (Baroud et al. 2006;Loeffel et al. 2008), but none have provided a clear definition for 'high-' or 'low-viscosity' cement. From the results of this work, we can state that the marrow viscosity is the reference relative to which the cement viscosity must be higher, especially after the shear-thinning from the injection. The shear-thinning causes the cement and marrow viscosities to come closer to their respective ∞ compared to 0 . In this regard, the power law exponent n and the lower viscosity limit ∞ of the cement and the marrow are crucial. For practitioners, this information could be useful for choosing the bone cement and its curing time depending on the patient's marrow composition, or for making decisions like whether techniques like marrow aspiration need to be employed.
There are also some limitations to our current study. As mentioned earlier, the bone cement viscosity and curing behaviour are sensitive to mixing conditions and exposure to air. We also did not consider the temperature dependence, which would cause the viscosities to change from the heat released by curing, or the difference in the room and body temperature in a clinical setting. More limitations could arise in clinical settings, e.g. we applied boundary conditions such that the fluids can freely flow out of the walls, whereas in vivo, they might enter into the blood vessels inside the vertebra. Modelling these in vivo conditions would require dependence of the boundary conditions on the distribution of the blood vessels in the vertebra. Moreover, influential parameters like permeability and marrow viscosity are hard to determine in in vivo conditions. These parameters could be estimated from in vivo measurements like clinical CT images  or patient characteristics based on clinical data. However, these factors contribute to the uncertainty in the outcome of the simulation results. In this regard, a study for quantifying the uncertainties in the parameters could be quite beneficial. Simulations using such estimations could help practitioners estimate at least a pilot range for the operating parameters, e.g. injection pressure, cement viscosity, curing time, flow rate, etc., before performing the actual procedure instead of completely relying on haptic feedback and frequent X-ray images at the time of the procedure. The reduced reliance on X-ray imaging would decrease the patient's exposure to harmful radiation. The model could also serve as a training tool for relatively young, inexperienced practitioners. In all the simulations carried out in this study, the deformations were negligible: ≈ 10 −10 m for aluminium foam, ≈ 10 −9 m for trabecular bone. Despite this, the capability to compute deformations is advantageous in cases like a weakened vertebra where the deformations may not be negligible. Extending the model to include the effect of existing fractures on vertebroplasty is possible and especially clinically relevant for future work since vertebroplasty mostly happens on fractured vertebrae. Similarly, the inclusion of temperature effects in the model has relevance since tissue necrosis is another risk associated with vertebroplasty. Lastly, the physics of this model can be suitably modified for other injection and infusion processes in the biomedical field.

Conclusion
We presented a continuum-mechanical model based on the Theory of Porous Media for simulating vertebroplasty. The model can simulate the injection of bone cement in porous materials like vertebrae and yield realistic results, as was evident from the benchmark experiment used for validation. The average viscosity model is found to be the most suitable approach for upscaling the rheology in the macroscale framework. The Cannella model could potentially be used, but the Hirasaki and Pope model is not suitable. Our simulations show that the cement must have a higher viscosity than the marrow to ensure stable development of injection force and proper filling of cement. In this regard, the marrow's and the cement's relative rheologies play a crucial role in the outcome of vertebroplasty and in avoiding cement leakage. We expect the current model to support future developments of vertebroplasty simulations that are closer to clinical reality and expect its possibilities to be extended towards modelling other fluid injection mechanisms in biomedical fields.