```
import porespy as ps
import numpy as np
import matplotlib.pyplot as plt
from edt import edt
```

# Adding Gravity to Image Based Drainage Simulations

The PMEAL team recently developed a way to add gravitational effects to the standard drainage simulation algorithm based on sphere insertion, also sometimes called *morphological image opening* or *full morphology*. The paper is available in Water Resources Research. This post will give a brief overview of how to use this new algorithm, which is being included in a new function called `drainage`

located in the also newly added `porespy.simulations`

module. The new `drainage`

function can produce results with or without gravity effects, creating a more powerful yet flexible function for performing such simulations.

Start by generating a 2D image of blobs:

```
= ps.generators.blobs([600, 600], porosity=0.7, blobiness=1.5)
im ; ps.imshow(im)
```

First we’ll compute a reference capillary pressure curve with no gravity effects using the new `drainage`

function. This function requires us to supply an image of capillary pressure values associated with each voxel. Typicall one would use the Washburn equation:

\[ P_C = -2\sigma \cos (\theta) \bigg(\frac{1}{r}\bigg) \]

where `r`

is conveniently provided by the distance transfer. This can be computed as follows:

```
= 1e-4
voxel_size = edt(im)
dt = -2*0.072*np.cos(np.deg2rad(180))/(dt*voxel_size) pc
```

The reason for requiring `pc`

instead of just the boolean `im`

is to allow users to compute the capillary pressure using any method they wish. For instance a 2D image of a glass micromodel has a spacing between the plates which can be accounted for as:

\[ P_C = -\sigma \cos (\theta) \bigg(\frac{1}{r} + \frac{1}{s}\bigg) \]

Note that `s`

could be set to `np.inf`

for a 2D simulation. Computing `pc`

externally also reduces the number of arguments that must be sent to the function, by eliminating the surface tension and contact angles.

In any event, the function is called as follows:

```
= np.zeros_like(im)
inlets 0, ...] = True
inlets[= ps.simulations.drainage(pc=pc, im=im, inlets=inlets, voxel_size=voxel_size, g=0)
drn
ps.imshow(np.log10(drn.im_pc))= plt.colorbar()
cbar 'log10(Pc)'); cbar.set_label(
```

The above image is colored according to the capillary pressure at which a given voxel was invaded. Note that we have applied a base-10 logarithm to the values to improve visibility. This format is handy since it contains the entire invasion sequence in a single image. Next let’s plot the capillary pressure curve. There are some functions available in `porespy.metrics`

for this, but the `drainage`

function actually does it for us and attaches the result to the returned object:

```
'b-o')
plt.semilogx(drn.pc, drn.snwp, 'Capillary Pressure [Pa]')
plt.xlabel('Non-wetting Phase Saturation'); plt.ylabel(
```

Now let’s perform the same experiment but including the gravity effects:

```
= ps.simulations.drainage(pc=pc, im=im, inlets=inlets, voxel_size=voxel_size, g=9.81)
drn2
ps.imshow(np.log10(drn2.im_pc))= plt.colorbar()
cbar 'log10(Pc)'); cbar.set_label(
```

We can see from the colormap that the invasion pattern is different. Plotting the capillary pressure curve will make the differences more obvious:

```
'r-o')
plt.semilogx(drn2.pc, drn2.snwp, 'Capillary Pressure [Pa]')
plt.xlabel('Non-wetting Phase Saturation'); plt.ylabel(
```

Combining both plots together for a better comparison:

```
'b-o')
plt.semilogx(drn.pc, drn.snwp, 'r-o')
plt.semilogx(drn2.pc, drn2.snwp, 'Capillary Pressure [Pa]')
plt.xlabel('Non-wetting Phase Saturation'); plt.ylabel(
```

The red curve above is shifted to the right with a lower slope. This is the expected result since the impact of gravity makes it more difficult for the heavy non-wetting fluid to rise up the domain and penetrate into pores. The decresed slope is also expected since pores are invaded more gradually rather than in a larger percolation event.