( previous, home, next )

A dribble castle at the beach

During summers as a little kid, my friends and I liked to build sand castles at the beach. We would sometimes decorate the castles by putting little dribble towers on the top. We did this using a plastic funnel that would drip out sand and water from the bottom. However, getting the funnel to work right was a little tricky. If the sand I used was too fine and closely packed, the water would not drain down from the top through the funnel and no sand or water would come out the bottom.

The funnel of sand and water is one small example of a very general question – how does liquid flow through a porous medium? Solutions of this question can be useful for fields from the design of sponges and insulation to well water, and oil and gas extraction. Solutions can even help us understand the spread of malware over the internet and the Ebola virus in Africa. But to answer it, we have to develop new kinds of models that we haven’t studied before – models that include space.

How to define “space”

The concept of “space” is a rather flexible one in mathematics. There are many kinds of spaces: vector spaces, Euclidian spaces, polish spaces, … all of which share some aspect of our intuition about what a space is.

Intuitively, we interact and conceptualize space though are ability to move through it. Space is effectively one-dimensional for a tight-rope walker who can only move forward or backward. We live most of our lives in two dimensions, being able to move forward and backward and being able to turn left and right. But nature itself is concretely 3 dimensional, many animals like fish and birds can exploit this and we ourselves have invented devices from stools to rockets that allow us to exploit this third dimension.

In a spatial model, we have to explicitly define how things can move. A good example of how we do this is found in cellular automata models. A cellular automata model treats space like a tiled floor. Suppose you stand on the floor and occupy one tile, or “cell”. From your cell, you can take a step into each of the neighboring cells, and by repeatedly stepping from one cell to the other, you can move long distances.

But exactly how you move depends on the layout of tiles on the floor. You will have different numbers of neighbors to choose from, depending on the layout. We do not worry about this in real life – the idea of moving from tile-to-tile is only a rough approximation of real movement. But when building a cellular automata model, these details are crucial. So, how do we define a “neighborhood” in two dimensions?

There are 4 conventional ways to define neighborhoods in a rectangular array. Each of these can be described with a “stencil” connecting the origin to each of its “neighboring” nodes in the Cartesian grid.

Square lattice
Square lattice
Triangular lattice
Triangular lattice
Moore neighborhood
Moore neighborhood
Hexagonal lattice
Hexagonal lattice

Site Percolation on a random square lattice

In any real setting, percolation phenomena depend on the chemical and physical composition of the medium, the viscosity of the fluid, temperature, pressure, and even erosion and contamination. But the core concepts of percolation can be qualitatively modeled in very simple and illuminating ways. Along the way, we’ll get an example of how something interesting can emerge from something very simple.

Suppose we start with a rectangular array of square tiles, randomly black and white. Let \(p\) be the probability that a given site is empty and marked white while \(1-p\) is the probability that a given site is occupied and marked black. The variable \(p\) is called the “void fraction”, and \(1-p\) the “fill fraction”. The state of each site is independent of all the other sites. Such a lattice can be generated with just a few lines of python code.

[Show code]
p = 0.9
n = 16
A = ceil(rand(n,n) - p)

We can use this kind of code to create a menagerie of random media by tuning the void fraction \(p\) and the size \(n\). Let’s start by looking at some examples of random arrays. Starting at the top, can we move through white-space and reach the bottom?

[Show code]
#!/usr/bin/env python2

from scipy import *
from matplotlib.pyplot import *

for n in [16,32]:
    i = 1
    for p in [0.1, 0.3, 0.4, 0.5, 0.6, 0.9]:
        A = floor(rand(n,n)+p)
        subplot(2,3, i)
        i += 1
        xlabel('p=%.2f'%(1-p), fontsize=18)
    savefig('siteperc_n=%02d.png'%n, bbox_inches = 'tight', transparent=True)

for p in [0.3, 0.4, 0.5]:
    A = floor(rand(n,n)+p)
    subplot(1,3, i)
    i += 1
    xlabel('p=%.2f'%(1-p), fontsize=18)

savefig('siteperc_n=%02d.png'%n, bbox_inches = 'tight', transparent=True)

n = 16 cells

Well, if \(p\) is close to 1, it is easy – there is lots of open white space to move through. And for \(p\) close to zero, the answer is clearly no. But for values of \(p\) in the middle, things are a little trickier.

what directions can we move in? Can we go diagonal? By convention, let us assume can only travel left, right, up, or down, in a von Neumann neighborhood around each cite.

n = 32 cells

n = 64 cells

We can compute how far down we go – the “depth” – by hand, but that can get hard when \(n\) gets large! We can also compute the depth algorithmically, and python (or any language of your choice) to implement that algorithm.

The algorithm uses recursion.

[Show code]
#!/usr/bin/env python

from scipy import *

def printbinarymatrix(A):
    m,n = A.shape
    f = lambda i : (" " if i < 0.5 else "#")
    s = newline.join([ "|%s|"%"".join([f(i) for i in row]) for row in A])
    s = newline.join(["+%s+"%("-"*n),s,"+%s+"%("-"*n)])

def calc_percolation_depth(A, i0=None):
    Given an array A of {0,1}'s, find the farthest across the
    matrix we can walk on 0-entries, when
    using a von Neumann neighborhood, if we start
    at the left side in row i0.

    If i0 == None, return the maximum over all initial rows.

    upper_bound = A.shape[1]
    if i0 == None:
        return max([calc_percolation_depth(A,i) for i in range(upper_bound)])

    neighbors = [(1,0), (-1,0), (0,1), (0,-1)]

    # to track which sites we have visited ...
    visited = 0*A
    def next(i,j):
        if i < 0 or i >= A.shape[0] or j < 0:
            return 0
        if j >= upper_bound:
            return upper_bound
        if visited[i,j] == 1:
            return 0
        visited[i,j] = 1
        if A[i,j] > 0.5:
            return j
        return max([next(i+dx,j+dy) for (dx,dy) in neighbors])
    return next(i0,0)

def test():
    p = 0.45
    N = 12
    A = floor(rand(N,N) + p)
    assert set([0,1]) == set(map(int,A.flatten()))
    print("Here we have a %dx%d matrix, with fill fraction %f"%(N,N,p))
    print("If we start on the top, how far toward the bottom can we move on blanks?")
    for i in range(N):
        d = calc_percolation_depth(A,i)
        print("From column %2d, we go %2d steps down"%(i,d))
    print("Max depth: ", calc_percolation_depth(A))
    print("Now, double-check with your own eyes.  Are all of these correct?")


Here we have a 12x12 matrix, with fill fraction 0.450000

If we start on the top, how far toward the bottom can we move on blanks?

From column  0, we go  2 steps down
From column  1, we go  2 steps down
From column  2, we go  2 steps down
From column  3, we go  0 steps down
From column  4, we go  0 steps down
From column  5, we go  0 steps down
From column  6, we go 12 steps down
From column  7, we go 12 steps down
From column  8, we go  0 steps down
From column  9, we go  3 steps down
From column 10, we go  3 steps down
From column 11, we go  3 steps down
Max depth:  12

Now, double-check with your own eyes. Are all of these correct?

Note above that water will get through for small \(p\), and that it will not for large \(p\). And interestingly, it seems that as the we look at larger matrices, the chance of getting through becomes less-random and more predictable based on \(p\). Very important question: For very large systems, is there a specific value \(p=p_c\) when the water stop flowing?

Percolation threshold for site-percolation on a von-Neumann neighborhood, \(p _ c \approx 0.5927\) although this is based on crude estimates and it remains a great open problem of statistical mechanics to find an efficient formula for this number.

But we should also emphasize that this result would change if we changed the neighborhood of allowed movement. If we could move in diagonal directions as well, like a Moore neighborhood, then it would be easier to move through the medium, and the percolation threshold would be smaller.

Bond percolation with von Neumann neighborhoods

Bond percolation is a different geometric approach to a similar problem. In bond percolation, the parameter \(p\) controls the probability that two adjacent nodes are connected along an edge. So, it’s not the cells themselves that are important, but rather the edges connecting the cells.

Instead of each cell being empty or filled, in bond percolation, each edge is present with some probability \(p\) or absent with probability \(1-p\). Here is an example of an \(80 x 80\) matrix with \(p=0.45\).
A bond model of percolation

Can you travel from top to bottom along the bonds?

Despite the similarities to site percolation, the difference in geometry has a huge impact on our ability to answer this question. For von Neumann bond percolation, we can compute \(p_c\) exactly as the lattice becomes infinitely large. The threshold below which you can’t travel from top to bottom \(p_c=1/2\)!

To explain why \(p_c\) must be \(1/2\), we need to introduce the concept of our bond lattice’s dual lattice. Imagine that all bonds were present, so we had the full square lattice. At the center of each elemental square in the full square lattice, create a new node. We will call this set of new nodes the “dual nodes”. The neighbors of the dual nodes are also given by the von Neumann neighborhood. Now, for any two neighboring dual nodes, connect the nodes with an edge if and only if the new edge will not cross any of the edges in our original bond lattice. The lattice made up of the dual nodes and their edges is called the dual lattice. Because of the conditions for edge existence, the dual lattice will have edges present with probability \(1-p\). Here is an example.

Example of a random-bonds von Neumann lattice with dual lattice filled in

Now, let’s consider percolation on our lattice from top-to-bottom, simultaneously with percolation of the dual lattice from left-to-right. If our lattice percolates, then there will be atleast one unbroken path of bonds from top to bottom. This implies that the dual lattice can not have an unbroken path of bond from left-to-right, and hence, the dual lattice can not percolate. On the other hand, if the dual lattice percolates from left to right, then our original lattice will not be able to percolate.

Imagine that \(p _ c < 1/2\) for bond percolation with von Neumann neighborhoods. Pick \(p\) so \(p _ c< p < 1/2\). Then both the original and the dual lattices must have percolate, but this is a contradiction to our observation that only one lattice can percolate at a time. Thus, \(p _ c \geq 1/2\).

On the other hand, if \(p _ c > 1/2\), then if we pick \(1/2 < p < p _ c\), the lattice will be below the percolation threshold, and thus will not contain a connected path from top to bottom. This implies there is a connected path in the dual lattice separating top from bottom. Thus, the dual lattice percolates. But since \(1/2 < p < p _ c\), \(1-p < 1/2 < p _ c\), so the dual lattice is below percolation threshold and should not percolate either. Thus, a contradiction.

It follows then that the only possibility is \(p _ c = 1/2\).


  1. The “neighborhood” in a cellular automata model defines the geometry of the space. In class, we discussed 4 different kinds of neighborhoods for 2-dimensional cellular automata on a Cartesian grid – 8-neighbor Moore, the 6-neighbor, the 4-neighbor von Neumann, and the 3-neighbor.

    1. Let \(N_k(x,y)\) be the set of neighbors of site \((x,y)\) where \(x\) and \(y\) are integers. If \(k = 3\), then \[ N_3(x,y) = \begin{cases} \{ (x+1,y), (x,y+1), (x,y-1) \} & \text{ if $x+y$ is even}, \\ \{ (x-1,y), (x,y+1), (x,y-1) \} & \text{ if $x+y$ is odd}. \end{cases} \] Find similar formulas for \(N_4(x,y)\), \(N_6(x,y)\), and \(N_8(x,y)\).

    2. In a regular lattice, the atomic loop length is the smallest number of neighboring edges in loop from \((0,0)\) back to itself where the same edge between two neighbors is never traversed more than once. Each such loop with this minimal loop length is called an “atomic loop”. For each \(k \in \{ 3,4,6,8 \}\), find the atomic loop length and the number of atomic loops containing \((0,0)\) for lattices with neighborhoods \(N_k(x,y)\).

    3. Given a neighborhood \(N_k(x,y)\), we can define a metric \(d_k( (x,y), (u,v))\) to measure the distance between points \((x,y)\) and \((u,v)\) recursively as follows: \[d_k( (x,y), (u,v) ) = \begin{cases} 0 & \text{if $(x,y) = (u,v)$}, \\ 1 + \min \{ d_k( (w,z), (u,v) ) : (w,z) \in N_k( (x,y) ) \} & \text{otherwise}. \end{cases}\] For each \(k \in \{ 3,4,6,8 \}\), draw the set of points \(\{ (u,v) : d_k( (0,0), (u,v) ) \leq 2\}\) on the appropriate lattice.

    4. For each \(k \in \{ 3,4,6,8 \}\), find \(d_k( ( 0, 0), (3,3) )\).

  2. In class, we were able to visually calculate percolation depth because of our awesome pattern analysis wet-ware. Calculating percolation depth algorithmically is a little more complicated, but can be done recursively. See here

    1. Let \(A(p,N)\) be a random \(0/1\) matrix with shape \(N\times N\) where sites are filled (\(1\)) with probability \(p\), and empty (\(0)\) with probability \(1-p\). Recall that in class, we generated example matrices like this with the python code \[ A = \text{floor(rand(N,N) + p)} \] For each value of \(p \in \{ 0.3, 0.4, 0.45, 0.5, 0.6 \}\), simulate \(1000\) \(40 \times 40\) matrices. Let \(g(x,p)\) be the fraction of these matrices with percolation depth less than or equal to \(x\). Plot \(g(x,p)\) for each value of \(p\), all in one plot. Remember to label your plot axes and include a legend.

    2. What particular feature of your plot change when \(p\) is between \(0.45\) and \(0.5\)? What does this change mean?

    3. Suppose we define \(h(p,N)\) as the fraction of \(N\times N\) \(A(p,N)\) matrices which percolate all the way through. Plot \(h(p,N)\) as a function of \(p \in [0.3,0.6]\) (this is an INTERVAL, not a set!) for \(N \in \{ 5, 10, 20, 50, 100 \}\). Use atleast 1000 matrices for each. (Warning: This calculation may take you a long time.)

    4. Extrapolating from your plot, what do you think will happen to \(\lim_{N\rightarrow \infty} h(p,N)\)?

    5. Reality, of course, is often three-dimensional, rather than two-dimensional. Describe how you think the percolation threshold will change when we switch from two dimensions to three dimensions and why it will change.

  3. Estimate the percolation threshold

( previous, home, next )