Background and assumptions

A complex process such as glacial ripping is influenced by a wide range of parameters, each with a range of possible values. To focus our modelling, this range of parameters needs to be restricted: some parameters are not explicitly modelled, but are discussed. Field evidence suggests that damage from glacial ripping was particularly effective below the ablation zone, close to the retreating margin of the Pleistocene ice sheets, and thus constituted an intensive phase of subglacial erosion just prior to deglaciation (Hall and others, Reference Hall2020, Reference Hall, Mathers and Krabbendam2021; Bukhari and others, Reference Bukhari2021; Krabbendam and others, Reference Krabbendam, Palamakumbura, Arnhardt and Hall2021, Reference Krabbendam, Hall, Palamakumbura and Finlayson2022). We thus focus on glaciological conditions beneath the last FIS during its late stage of final deglaciation. Direct glaciological controls exerted by this palaeo-ice sheet are evidently not available, but the Greenland ice sheet (GrIS) is also in a state of retreat, so we use observations from the ablation zone of the GrIS as an analogue for the retreating Pleistocene ice sheets.

Relevant assumptions are as follows:

1. (1) Basal thermal regime. The ice-sheet bed was ‘warm’ or thawed. Abundant abraded surfaces decorated with striae show that abrasion (and hence warm-based sliding) was active prior to deglaciation in eastern Sweden (e.g. Sohlenius and others, Reference Sohlenius, Hedenström and Rudmark2004); ice-sheet models also simulate warm-based conditions close to the retreating margin (e.g. Näslund and others, Reference Näslund, Rodhe, Fastook and Holmlund2003; Patton and others, Reference Patton2017). Similarly, the base of the GrIS ablation zone is thawed (Macgregor and others, Reference MacGregor2016), and ice motion shows a significant component (45–75%) of basal sliding (Ryser and others, Reference Ryser2014).

2. (2) Ice viscosity. We model ice as a viscous Newtonian medium, with the viscosity as established experimentally by Byers and others (Reference Byers, Cohen and Iverson2012). This makes the treatment simpler (see also Nye, Reference Nye1970; Hallet, Reference Hallet1979), but is also realistic for temperate ice (Colbeck and Evans, Reference Colbeck and Evans1973; Chandler and others, Reference Chandler, Hubbard, Hubbard, Murray and Rippin2008; Byers and others, Reference Byers, Cohen and Iverson2012; Krabbendam, Reference Krabbendam2016). In essence, we thus assume a layer of temperate ice, thicker than the highest rock obstacles. In Greenland, borehole temperature measurements show a 30–100 m thick basal layer of temperate ice at 20–50 km from the margin, whereas close to the margin the entire ice thickness is temperate (Lüthi and others, Reference Lüthi, Funk, Iken, Gogineni and Truffer2002; Ryser and others, Reference Ryser2014; Harrington and others, Reference Harrington, Humphrey and Harper2015; Harper and others, Reference Harper, Meierbachtol and Humphrey2019). The effects of cold ice, with its different rheological behaviour, is discussed but not modelled.

3. (3) Ice sliding velocity. Patton and others (Reference Patton2017) modelled maximum ice surface velocities of 200–400 m a⁻¹ during the deglaciation of eastern Sweden. In the GrIS ablation zone, basal sliding velocities are c. 50–75% of the surface velocities (Ryser and others, Reference Ryser2014), so we assume a value of ~300 m a⁻¹ as a maximum realistic sliding velocity.

4. (4) Ice thickness. Although the FIS reached a maximum thickness of c. 3000 m during late glacial maximum conditions (Quiquet and others, Reference Quiquet, Colleoni and Masina2016), field evidence suggests that glacial ripping was particularly active just prior to deglaciation, hence with low ice thickness. In the lower ground in eastern Sweden, ice retreat was dominated by calving in the Baltic Ice Lake, with water depth up to 180 m. We assume an ice thickness of 300 or 600 m.

5. (5) Basal water pressure fluctuations. Evidence for hydraulic jacking includes dilated and sediment-filled fractures (see section ‘Setting and summary of field evidence’). High water pressures, including overpressure when temporarily exceeding overburden pressure (P_i), have been measured or demonstrated at the base of the GrIS ablation zone in the following settings. Firstly, dramatic supraglacial lake drainage events can lift up the ice surface over several square kilometres (Das and others, Reference Das2008; Doyle and others, Reference Doyle2013) and these represent high-magnitude overpressure events, involving volumes >10⁶ m³ water, resulting in substantial ice–bed separation. Secondly, high-frequency (daily in the melt season), but short-lived (hours), water-pressure fluctuations between 80 and 110% of P_i have been measured in boreholes (Andrews and others, Reference Andrews2014; Claesson Liljedahl and others, Reference Claesson Liljedahl2016; Wright and others, Reference Wright, Harper, Humphrey and Meierbachtol2016; Harper and others, Reference Harper, Meierbachtol and Humphrey2019). The long-term average water pressure in these boreholes is c. 90–95% of P_i, although some boreholes show long-term close to or exceeding 100% of P_i. These water pressure fluctuations are highly localised events (out-of-phase with adjacent boreholes) and thus represent low-magnitude events, involving only small volumes of water, with only localised ice–bed separation. However, they occur repeatedly (tens of times) throughout the melting season. Thirdly, Andrews and others (Reference Andrews2014) noticed an additional setting, where moulins are well-connected to the glacier bed, preventing overpressure to build up. In this regime, frequent water pressure fluctuations occur between 60 and 98% of P_i, with a long-term average of c. 80% of P_i. All three processes are likely to have occurred during deglaciation of the FIS and may have led locally to hydraulic jacking in eastern Sweden. Thus, for modelling purposes we assume that high water pressure events, including overpressure events, occurred at the ice bed near and around the idealised obstacles under discussion. We do not model or discuss the specific dynamics or spatial extent of high water-pressure events in Sweden. We further assume that the normal stress at the ice bed equals the overburden pressure of the ice.

6. (6) We make the simplifying assumption that the fracture patterns of the bedrock are broadly orthogonal, and dominated by subhorizontal and subvertical fractures. Such patterns are common in basement rocks in Sweden and elsewhere (Talbot and Sirat, Reference Talbot and Sirat2001; DeWandel and others, Reference Dewandel, Lachassagne, Wyns, Maréchal and Krishnamurthy2006; Pless and others, Reference Pless2015; Krabbendam and others, Reference Krabbendam, Palamakumbura, Arnhardt and Hall2021), and near ubiquitous in flat-lying sedimentary rocks. However, we do model two different fracture patterns: (a) a scenario where all fractures are continuous in-plane, and (b) a scenario where subhorizontal fractures occur at different levels (‘stepped’) but joined to subvertical fractures, via T-junctions, thus being discontinuous in a single plane. This is similar to the ‘step-path’ failure mechanisms recognised in large rock-slope failures (e.g. Brideau and others, Reference Brideau, Yan and Stead2009). In-plane rock bridges along a single fracture (Kemeny, Reference Kemeny2003; Hooyer and others, Reference Hooyer, Cohen and Iverson2012; Elmo and others, Reference Elmo, Donati and Stead2018), and parts of the fracture sealed by fracture fills (Shang and others, Reference Shang, Hencher and West2016), also represent discontinuous fractures and may hamper hydraulic jacking. During overpressure events in which hydraulic jacking occurred, some rock bridges were likely broken by hydraulic fracturing, forming longer, more continuous fractures.

Modelling scenarios

The geometry of the ripped bedrock hills, their internal fracture networks and the resultant effects of glacial ripping in eastern Sweden are highly variable. We confine ourselves to model a number of simple scenarios, informed by natural examples, but including some end-member scenarios (Figs 3a–e):

1. (1) a blunt hemispherical obstacle of intact rock, without a basal fracture;

2. (2) a blunt hemispherical obstacle with a continuous, subhorizontal basal fracture;

3. (3) a blunt hemispherical obstacle with a discontinuous basal fracture, formed by steps in the basal fracture system;

4. (4) an elongate obstacle with subhorizontal fractures and a blunt stoss side, but a large flat top;

5. (5) a small rock step, caused by the differential uplift of a rock block, which then protrudes above the surrounding rock surface. This scenario is applicable to the smooth, abraded top surfaces of large hills and low relief rock surfaces.

Fig. 3. Modelling scenarios 1–5: conceptual geometries; modelled geometries with some parameters indicated – see also Table 1. Figure © Svensk Kärnbränslehantering. Reproduced with permission.

The basis of the modelling is the principle that if the drag force F _d exerted by ice flowing past a blunt obstacle exceeds the resisting forces F _r of that obstacle, then the obstacle moves or disintegrates:

(1)$$F_{\rm d} > F_{\rm r}$$

All variables and constants are also shown in Table 1.

Table 1. Variables, constants and parameters used

Drag forces for a hemispherical obstacle

Temperate ice is modelled as a viscous Newtonian medium creeping at low velocities (Hallet, Reference Hallet1979; Byers and others, Reference Byers, Cohen and Iverson2012). In this laminar flow regime with very low Reynolds number (Re ≪ 1), also known as Stokes regime, the drag force F _d on a spherical particle has been obtained by solving the Navier–Stokes equations (Stokes, Reference Stokes1951; Loth, Reference Loth2008):

(2)$$F_{\rm d} = 6\pi \eta Ur$$

where η is the viscosity, U the velocity of the medium and r the radius of a spherical particle. F _d is a combination of form drag, which depends on the shape of the particle because it arises from the pressure the fluid exerts on the cross-sectional area of the object perpendicular to the streamlines, and viscous friction (skin) drag, which is due to the tangential shear stress at the particle surface (Leith, Reference Leith1987). In a normal viscous medium, the creep velocity right at the contact is zero (non-slip boundary), and the tangential shear stress at the obstacle surface is substantial. This is patently not the case for a thawed ice–bed contact, where sliding is known to occur with a thin film of water between the obstacle and the ice. Thus, for temperate ice creeping past an obstacle, skin drag may be low. For this reason, some authors suggest that F _d = 4πηUr (e.g. Hallet, Reference Hallet1979; Cohen and others, Reference Cohen2005). However, Byers and others (Reference Byers, Cohen and Iverson2012) empirically confirmed that Eqn (2) is broadly valid (with values for the numerical prefactor of 5.5π; 6.9π and 7.5π) for temperate ice. Since we also use the experimental values of ice viscosity from the same experiments (see above), we herein use Eqn (2) as is. However, we cut the sphere in half to model a hemispherical obstacle (Fig. 4) so that:

(3)$$F_{\rm d} = 3\pi \eta Ur$$

Fig. 4. Principle of applying Stokes law for laminar flow around a sphere to a hemispherical obstacle: (a) laminar flow around a sphere and (b) hemispherical obstacle, in a half space, within a laminar flow field. Area of projected stoss side A _st and area of footprint of obstacle A_xy are indicated. Direction of ice flow = x. Figure © Svensk Kärnbränslehantering. Reproduced with permission.

We ignore the effects of regelation because the large size of the rock obstacles compared to smaller debris particles diminishes the heat flow effect through the obstacle and thus renders regelation inefficient over length scales greater than a few centimetres (see also Hallet, Reference Hallet1979; Cuffey and Patterson, Reference Cuffey and Patterson2010; Byers and others, Reference Byers, Cohen and Iverson2012).

Note that the ‘sliding laws’ as derived by for instance by Nye (Reference Nye1970; see discussion in Cuffey and Patterson, Reference Cuffey and Patterson2010), are broadly similar in concept, except that these sliding laws integrate obstacle size over a number of obstacles and thus employ a roughness parameter rather than the size of an individual obstacle, and use the overall basal shear stress, rather than the force exerted on an individual obstacle.

Drag forces for a non-hemispherical, elongate obstacle

Viscous drag forces for a non-hemispherical obstacle are more complicated since they depend on the shape of the object, which needs to be quantified with a shape descriptor, which treats the form and skin drag separately (e.g. Leith, Reference Leith1987; Ganser, Reference Ganser1993; Bagheri and Bonadonna, Reference Bagheri and Bonadonna2016; Dioguardi and others, Reference Dioguardi, Mele and Dellino2018). These approaches all assume a non-slip boundary, and thus are not appropriate in our case, because there is certainly slip along a thawed ice–bed contact. To approximate the drag forces on elongate obstacles, we model a geometry where we split the obstacle in (1) a blunt stoss side that faces ice flow, and (2) a flat horizontal top surface parallel to ice flow (Fig. 3d). We assume that (1) the drag force F _vd at the stoss side is controlled by the viscous drag acting upon that stoss face (with surface A _st) and controlled by a form of Stokes law and that (2) the drag forces F _cf on the flat top surface area are controlled by a simple friction law acting on the top surface (with surface A_xy). This ice–rock friction is complex in detail, and different from normal rock–rock friction, because of the presence of rock debris particles between the ice and the bed, variations in concentration and particle size of the debris at the ice bed, and ice deformation and melting around the particles (Hallet, Reference Hallet1979; Emerson and Rempel, Reference Emerson and Rempel2007; see also discussion in Schweizer and Iken, Reference Schweizer and Iken1992). Nevertheless, integrated over a large area, the concept of bulk Coulomb friction appears to be valid (Cohen and others, Reference Cohen2005; Emerson and Rempel, Reference Emerson and Rempel2007; McCarthy and others, Reference McCarthy, Savage and Nettles2017), and is used herein:

(4)$$F_{{\rm cf}} = \mu _{{\rm ir}}\;F_{{\rm iz}}$$

where μ _ir is the bulk friction coefficient on the ice–rock contact. We take μ _ir at 0.05, at the lower end of experimentally obtained values (from Cohen and others, Reference Cohen2005; Emerson and Rempel, Reference Emerson and Rempel2007; McCarthy and others, Reference McCarthy, Savage and Nettles2017); effects of higher friction coefficients as measured by Emerson and Rempel (Reference Emerson and Rempel2007) are not modelled but will be discussed. The normal force F _iz exerted by the overlying ice is a function of thickness and density of ice (cryostatic pressure P _i), the water pressure P _w and the horizontal surface area A_xy over which it operates as follows:

(5)$$F_{{\rm iz}} = ( P_{{\rm i\;}}-P_{{\rm w\;}}\;) \;A_{xy}$$

The horizontal surface A_xy is approximately equal to the surface area of the basal fracture of the rock obstacle. As we wish to explore the effects of relative overpressure P _w/P _i (in essence the potential effects of hydraulic jacking), this is rewritten as:

(6)$$F_{{\rm iz}} = gh_{\rm i}\rho _{\rm i}\left({1-\displaystyle{{P_{\rm w}} \over {P_{{\rm i\;}}}}} \right)A_{xy}$$

where g is the gravitational constant, h _i the ice thickness and ρ _i the density of ice. Total friction on the flat surfaces then becomes:

(7)$$F_{{\rm cf}} = \mu _{{\rm ir}}\;gh_{\rm i}\rho _{\rm i}\left({1-\displaystyle{{P_{\rm w}} \over {P_{{\rm i\;}}}}} \right)A_{xy}$$

Note that Coulomb friction F _cf will approach zero if P _w approaches P _i; since a frictional force cannot be negative, Eqns (5)–(7) are thus only valid for P _w ≤ P _i.

The viscous drag F _vd acting on the stoss side can be approximated by assuming that it is controlled by the square root of the vertical surface area that faces up-ice (the projected ‘stoss-side surface area’ A _st, which is half the surface area of a circle with radius r, or the width W multiplied by height H for a rectangular obstacle):

(8)$$A_{{\rm st}} = \displaystyle{1 \over 2}\pi \;r^{2\;}\quad {\rm so\;that}\,\,\;r = \sqrt {2\;A_{{\rm st}}/\pi \;} \;$$

The viscous drag component then becomes ((3) combined with (8)):

(9)$$F_{{\rm vd}} = 3\;\eta U\sqrt {2A_{{\rm st}}/\pi \;} $$

The total drag forces of a non-hemispherical obstacle is then ((7) and (9)):

(10a)$$F_{\rm d} = 3\eta U\sqrt {2\;A_{{\rm st}}/\pi \;} + \mu _{{\rm ir}}gh_{\rm i}\rho _{\rm i}\left({1-\displaystyle{{P_{\rm w}} \over {P_{\rm i}}}} \right)A_{xy}\;$$

For a rectangular cuboid-shaped obstacle, with width W, height H and length L this becomes:

(10b)$$F_{\rm d} = 3\;\eta U\sqrt {2HW/\pi \;} + \mu _{{\rm ir}}gh_{\rm i}\rho _{\rm i}\left({1-\displaystyle{{P_{\rm w}} \over {P_{{\rm i\;}}}}} \right)WL\;$$

Resisting forces for a hemispherical obstacle without fractures

If the rock obstacle is poorly fractured, for instance if the rock hill is smaller than the vertical fracture spacing in the local bedrock, or if no continuous subhorizontal fractures are present (Fig. 3a), the resisting force of intact rock F _rr is controlled by the intact rock strength as follows:

(11)$$F_{{\rm rr}} = \tau _{\rm r}A_{\rm r}$$

where τ _r is the shear strength of intact rock, and A _r the surface area of the potential shear plane of the intact rock. Shear strength values for intact rock are rarely measured, in contrast to uniaxial compressive strength (UCS) and tensile strength. Singh and others (Reference Singh2017) obtained shear strength of between c. 10 and 35 MPa for different gneisses. Shear strength is typically 20–30% of UCS, and about twice the tensile strength. Rock mechanic tests on rocks near Forsmark yielded a range of 160–370 MPa for UCS and 10–18 MPa for tensile strength (Glamheden and others, Reference Glamheden2007), suggesting shear strength range of 20–60 MPa. We take a conservative value of ~20 MPa for τ _r, as in nature rocks may contain micro-structures, which lower the shear strength.

Resisting forces for a hemispherical obstacle with continuous basal fracture

If a rock hill contains fractures – as most do – its rock mass strength is lower than that of intact rock. Rock mass strength is complex, so we consider the simple scenario of an obstacle with a single continuous subhorizontal fracture at its base (Fig. 3b). The resisting force F _rj along this fracture can be defined by:

(12)$$F_{{\rm rj}} = \mu _{{\rm rr}}F_{{\rm rz}\;}$$

where μ _rr is rock–rock friction coefficient along the basal fracture, and F _rz the normal force acting on the basal fracture. Rock–rock friction coefficient concerns the sliding of two rock blocks along a fracture (different from the ice–rock friction discussed previously) and varies with rock type and fracture roughness; for natural fractures in gneissic rocks Ramana and Gogte (Reference Ramana and Gogte1989) report a range of 0.64–0.77; Glamheden and others (Reference Glamheden2007) report a range of 0.48–0.77, with a mean of 0.67, for rocks at Forsmark. We take a value of 0.7. We ignore fracture cohesion, since hydraulic jacking will have broken any such cohesion.

The normal force F _rz has two components:

(13)$$F_{{\rm rz}} = F_{{\rm bz\;}} + F_{{\rm iz}}$$

where F _bz is the buoyant weight of the rock above the fracture and F _iz is any force (weight) exerted by the overlying ice (as per Eqn (6)). The drag force exerted by ice flowing vertically downwards due to basal melting is ignored: while it is potentially important for small (centimetre-sized) debris particles under conditions of fast basal melting (Cohen and others, Reference Cohen2005; Byers and others, Reference Byers, Cohen and Iverson2012), it becomes very small for metre-sized obstacles (Krabbendam and Hall, Reference Krabbendam and Hall2019).

The buoyant weight of the block is:

(14)$$F_{{\rm bz}} = V\;( {\rho_{\rm r}-\rho_{\rm w}\;} ) g\;$$

where V is the volume of the block, ρ _r the density of rock and ρ _w the density of water. The total resisting force is then (Eqns (13), (14) and (6)):

(15a)$$F_{{\rm rj}} = \mu _{{\rm rr}}\;\left[{V\;( {\rho_{\rm r}-\rho_{\rm w}} ) g + gh_{\rm i}\rho_{\rm i}\left({1-\displaystyle{{P_{\rm w}} \over {P_{{\rm i}\;}}}} \right)A_{xy}} \right]$$

For a hemispherical obstacle, as a direct function of r (with V = (2/3)πr ³ and A_xy = πr ²), this becomes:

(15b)$$F_{{\rm rj}} = \mu _{{\rm rr}}\;\left[{\displaystyle{2 \over 3}\pi r^3( {\rho_{\rm r}-\rho_{\rm w}} ) g + gh_{\rm i}\rho_{\rm i}\left({1-\displaystyle{{P_{\rm w}} \over {P_{{\rm i}}}}} \right)\pi r^2} \right]$$

For a rectangular cuboid-shaped obstacle, with width W, height H and length L, this becomes:

(15c)$$F_{{\rm rj}} = \mu _{{\rm rr}}\;\left[{H\;W\;L\;( {\rho_{\rm r}-\rho_{\rm w}} ) g + gh_{\rm i}\rho_{\rm i}\left({1-\displaystyle{{P_{\rm w}} \over {P_{{\rm i}}}}} \right)WL} \right]$$

Effect of variable transmissivity of basal fractures

Thus far, it is assumed that the water pressure within a basal fracture equalises instantaneously with pressure fluctuations at the ice–bed contact, and leads to hydraulic jacking if relative overpressure P _w/P _i, exceeds 1, but this only happens if the fracture is highly transmissive (in the hydrological sense). Fracture transmissivity, however, is extremely variable in nature. Low transmissive fractures, for instance tight fractures or fracture planes with many rock bridges, would dampen the fluctuating water pressures in the fracture, and hence lower the peak water pressures within the fracture (e.g. Neupane and others, Reference Neupane, Panthi and Vereide2020). The water pressure in such fractures will remain close to the long-term average (e.g. 80–95% of P_i, depending on the long-term glaciohydrology of the relevant sector in the ice sheet) but not rise above 100%, and thus not cause hydraulic jacking. Fracture transmissivity can vary in nature easily by many orders of magnitude. Instead, we introduce dimensionless transmissivity factor T _j to model the effects of variable transmissivity along basal fractures, shown here for a hemispherical obstacle (Eqn (15b)):

(16)$$F_{{\rm rj}} = \mu _{{\rm rr}}\;\left[{\displaystyle{2 \over 3}\pi r^3\;( {\rho_{\rm r}-\rho_{\rm w}} ) g + gh_{\rm i}\rho_{\rm i}\left({1-\displaystyle{{T_{\rm j}\;P_{\rm w}} \over {P_{{\rm i\;}}}}} \right)\pi r^2} \right]$$

Resisting forces with discontinuous basal fracture and some intact rock

Fracture networks in the basement rocks of eastern Sweden (and elsewhere) are highly variable, and continuous subhorizontal fractures at the base of an obstacle are not ubiquitous: in various quarries and natural outcrops it was observed that subhorizontal fractures occur at different levels and are hence discontinuous (or: ‘stepped'); in other sections subvertical fractures are dominant over subhorizontal fractures (Krabbendam and others, Reference Krabbendam, Palamakumbura, Arnhardt and Hall2021). This variability is modelled as a hemispherical obstacle where the base constitutes part intact rock and part fracture (Fig. 3c). The total resisting force F _r is then the resisting force of the fracture plus that of intact rock, proportional to the area occupied by the basal fracture A _j and intact rock A _r respectively:

(17)$$A_{xy} = A_{{\rm r\;}} + A_{\rm j}\;\;and\;F_{\rm r} = F_{{\rm rr\;}} + F_{{\rm rj}}$$

In proportion to area, this becomes:

(18)$$F_{\rm r} = F_{{\rm rr\;}}A_{{\rm r}\;} + F_{{\rm rj}}\;A_{\rm j} = F_{{\rm rr\;}}A_{xy}\;\;\displaystyle{{A_{{\rm r\;}}} \over {A_{xy}}} + F_{{\rm rj}}A_{xy}\;\displaystyle{{A_{{\rm j\;}}} \over {A_{xy}}}$$

So that, using Eqns (10) and (14):

(19a)$$F_{\rm r} = A_{xy}\;\tau _{{\rm rk\;}}\;\displaystyle{{A_{{\rm r\;}}} \over {A_{xy}}} + \displaystyle{{A_{{\rm j\;}}} \over {A_{xy}}}\;\mu _{{\rm rr}}\;\left[{V\;( {\rho_{\rm r}-\rho_{\rm w}\;} ) g + gh_{\rm i}\rho_{\rm i}\left({1-\displaystyle{{P_{\rm w}} \over {P_{{\rm i\;}}}}} \right)A_{xy}{\rm \;}} \right]$$

and hence for a hemispherical obstacle:

(19b)$$F_{\rm r} = \pi r^2\tau _{{\rm rock}}\;\displaystyle{{A_{{\rm r\;}}} \over {A_{xy}}} + \displaystyle{{A_{{\rm j\;}}} \over {A_{xy}}}\mu _{{\rm rr}}\;\left[{\displaystyle{2 \over 3}\pi r^3\;( {\rho_{\rm r}-\rho_{\rm w}\;} ) g + gh_{\rm i}\rho_{\rm i}\left({1-\displaystyle{{P_{\rm w}} \over {P_{{\rm i\;}}}}} \right)\pi r^2{\rm \;}} \right]$$