The Theoretical Diver – Page 5 – Theorizing about scuba diving

VPM-B: How to compute your deco

This is the first technical post about the inner workings of the VPM-B deco model. Here I will just present the algorithm (including the relevant formulas, some of them maybe in an unusual form) while the physical motivation or derivation will be the subject of a future post.

I try not to parrot the original explanations (which, honestly, I do not really understand) but rather put in words what I learned by staring long enough at the subsurface implementation of the model. As I explained earlier, we tried to make this an independent implementation but of course in the end we made sure it produces the identical deco schedules to the original FORTRAN code (that we took as the definition of the model).

First of all, it is not that different from the Bühlmann model, it rather builds on it. For the sake of presentation, I will assume you are familiar with how that model works (if not, I wrote a brief summary with some background in German, to be translated to English soon). In particular, it computes the same tissue loadings with inert gases following a simple diffusion equation

$$\frac{dp_i}{dt} = \alpha_i(p_{amb}-p_i)$$

where the time constant $\alpha_i$ is essentially the inverse tissue half-time.

The only difference is how to compute the current ceiling for a snapshot of tissue loadings (“tensions” in the deco theory speak) $p_i$: Bühlmann proposes an affine dependence on the ambient pressure (with the affine relation given in terms of the famous $a$ and $b$ coefficients. If you use gradient factors, this affine relation is further modified (depending on your lowest ceiling), but when you plot the ceiling as a function of $p_i$ it is still a straight line, i.e. affine.

This is different for VPM-B: Only in the zero’s approximation it is still a straight line. In fact, there the ceiling is simply given by a maximal gradient (difference between tissue pressure and ambient pressure) that is independent of the depth:

$$p_{amb} – p_i \le G_{0,i}.$$

This is called the initial gradient $G_0$. In this approximation, VPM-B is like Buehlmann, only with different a and b coefficients.

For the next approximation, the “Boyle compensation” we need to know a little bit about bubbles. in fact, the only thing one needs to know is that they have a surface and in the total energy balance of the bubble the surface contributes proportional to its area:

$$E_{surf}= \gamma r^2$$

for a bubble of radius r (we can absorb all kinds of numerical constants like $4\pi$ in the proportionality constant $\gamma$. Then the force on the bubble surface that wants to shrink it is

$$F_{surf} = \frac{\partial E_{surf}}{\partial r} = 2\gamma r.$$

Finally, this force acts on all of the surface, by dividing by its area again, we get the surface’s contribution to the pressure

$$G=\frac{2\gamma}r.$$

This type of formula, the a bubble surface contributing to the bubble’s internal pressure by $2\gamma/r$ you will see over and over again.

In particular, that the deepest point in the dive, or actually when reaching the first ceiling, the gradient $G_0$ is exactly this. And the Boyle compensation means that the (now depth dependent) maximal gradient (i.e. ambient pressure minus tissue pressure) is given by letting such a bubble expand according to Boyle’s law, i.e. by assuming that $pV$ is constant. $p$ inside the bubble (which is supposedly the same as the tissue pressure) is then given by the ambient pressure plus the momentary maximal gradient while the volume is proportional to $r^3$ which according to the above formula is proportional to $1/G^3$. So in total, the allowed difference between the tissue and the ambient pressure at depth $d$ is given by solving

$$\frac{G+p(d)}{G^3} = \frac {G_0+p_{1st\ ceiling}}{G_0^3}$$

So, rather than being a straight line, the maximal gradient is now the solution of a cubic equation (for which, in Subsurface, we use an analytic expression rather than running a Newton root finder as in the original implementation).

So far, I did not tell you what the initial gradient $G_0$ is. The VPM-B expression for it is

$$G_0 = 2 \frac\gamma{\gamma_c} \frac {\gamma_c – \gamma}{r_{reg}}.$$

Here, $r_{reg}$ is the crushing radius which is again determined from an expression of a bubble surface tension:

$$\frac{2(\gamma_c-\gamma)}{r_{reg}} = p_{max} + \frac{2(\gamma_c-\gamma)}{r_{crit}}$$

Here $\gamma$ and $\gamma_c$ are some surface tensions (effectively parameters of the model) as is $r_{crit}$. The latter is where one can tune the conservatism of the model: Each level of conservatism (like three levels in VPM-B+3) increases $r_{crit}$ by 5 percent. $p_{max}$ is the maximal ambient pressure encountered during the dive. In principle, it is also assumed that over time, $r_{reg}$ also exponentially approaches $r_{crit}$ again, but it is totally safe to ignore this since the decay time is supposed two weeks, far too long to play any role at the typical time scales of the dive.

Now, you have everything to compute a deco schedule: Like in the Buehlmann model, you know for each depth what the maximal gradient is.

Turns out, these schedules take far too long. This is where the “critical volume algorithm” comes into play. We now slightly increase the first gradient (the role so far played by $G_0$) which then is also subject to the Boyle compensation.

The idea here is that the body can tolerate an additional amount of gas going into bubbles. The maximal volume, again, is a model parameter and it is compared to the excess volume as follows: One assumes there is some distribution of bubble sizes. By enlarging the gradient $G_0$ to $G_0′$, the number of additional bubbles will be (this is one of the places where I do not understand the derivation, but let’s take the expression anyway) proportional to

$$\frac{G_0′-G_0}{G_0′-\alpha G_{crush}}.$$

Here, $\alpha$ is another constant and $G_{crush}$ is the maximal excess of the ambient pressure over the tissue pressure over the dive which also happens to be the pressure gradient over the surface of the bubble (if that is too big, we are in the “impermeable regime” and then there is another Boyle type compensation).

Furthermore, rate of additional out gasing is assumed to be proportional to the gradient $G_0′$ (totally ignoring the first Boyle compensation) and the whole thing lasts (assuming a constant gradient) over all of the decompression time (defined as the time from the first tissue starting to release gas again until we reach the surface). This decompression time $t_D$ is what we computed above where I said we can now compute a deco plan from the gradients.

Once we reach the surface, there is still an overpressure in the tissue compared to the ambient pressure, it starts with $G_0′$ (again ignoring the Boyle compensation) and then, according to the diffusion equation above, approaches 0 exponentially. so the time integral of the gradient yields $G_0’/\alpha_i$. Throwing all together, we get as a criterium for the maximal $G_0′$

$$V\ge \frac{G_0′-G_0}{G_0′-\alpha G_{crush}} G_0′ \left(t_D +\frac 1{\alpha_i}\right)$$

This we then solve for $G_0′$ which we can use for computing a new deco schedule which gives a new $t_D$ and we repeat this until $G_0′$ does not change anymore.

Apart from some minor fine print (like how in detail the deco schedules are computed) is all there is to VPM-B. With this information, you should, in principle be able to implement it yourself.

In the next post, I will look into the physics that is invoked to argue for some of the expressions that I only quoted here.

VPM-B: How it works (the no non-sense version)

One of the features that Subsurface users kept asking for were other deco models besides Bühlmann (including gradient factors). In particular, they wanted to see a bubble model implemented. I was reluctant to do so for some time since I believed (like many other divers) that you can produce fine decompression profiles that look like they are coming from a bubble model by using gradient factors, in particular by using settings like ${GF}_{low}=20-30$. But more important: I did not understand how they worked.

And this is not because I cannot understand how you can derive this or that formula from physical laws (I do that for a living), but because the text that are supposed to describe these models appeared to me to be confused. There were huge gaps in my understanding of the original texts, like this or this, those seemed to me like spending a lot of space on absolutely trivial calculations (like integrating an exponential function) while on the other hand not clearly saying what the actual strategy is and what all those symbols meant. Plus some of the derivations seems to use some slightly unorthodox algebra… There are other texts by other authors that are clearly derived from these originals but I keep having the impression that those secondary authors did not understand the difficult parts either, since many of them are just paraphrasing the originals.

For RGBM there is no hope of understanding it since its inner workings are a trade secret. But the situation of VPM-B is actually better: There is a very concise definition of the model. But it comes in the form of a FORTRAN program by Eric Baker. It is about 75 pages when printed out on paper and it is FORTRAN, so again, the core of the algorithm does not directly present itself to the reader.

Somebody had run this FORTRAN program through f2c but of course the resulting C code is not any more readable. But at least, that allowed be to single step through it with a debugger (rather ran running my finger through 75 pages of spaghetti code). Since then, there are a few more implementation, Since then there are a few more implementations including the C version by Bryan Waite which is at least idiomatic C. But still it follows pretty much in a 1:1 fashion the original FORTRAN code including the separation into sub routines and variable names.

But last summer, there was the opportunity to have a Slagvi, a Google Summer of Code student work on this and do a VPM-B implementation from scratch. I was his mentor and this here is what I learned. Slagvi’s task was not just to link the existing code to Subsurface but really do a reimplementation so that in the process we could understand what is actually going on (and what part of the “explanation” is just a smoke screen).

You can find the result in the subsurface source code, in particular in the two files deco.c and planner.c.

This post here is a start of a series of (to be written) posts that aim at a fresh presentation of how VPM-B works. The next looks at the actual algorithm.

What is this?

Welcome to my new theoretical scuba diving blog!

What is “theoretical diving” you will probably ask. Let me explain. By day, I work as a physicist at Ludwig-Maximilians-Universität München, as a theoretical physicist more specifically. But I also like a lot to go scuba diving. Unfortunately, I do not have as much time anymore to go diving as I used to have (for a variety of reasons). So instead, I like at least to think about diving. And I think about it as a theoretical physicist. In this blog, I would like to share what I came up with thinking about scuba diving.

Besides this blog, another theoretical diving activity of mine is contributing to Subsurface, the open source dive logging software started by Linus Torvalds and Dirk Hohndel. In particular, there I try to contribute to decompression calculations and the dive planner. In the course of working on that program, I also keep learning things (besides how to work in an open source project) that I would like to share here.

The plan to start this blog (besides my regular which means mainly physics blog) has been there for almost two years. So let’s see at which rate I manage to write posts. I hope I can also find some other contributors.

So, here we go…