Tag: VPM-B

The VPM-B model (of which I have described my understanding in previous posts) has a “derivation” in terms of bubble physics which I, as readers should have noticed, have some trouble to follow through completely as it is, to say the least, ad hoc and sketchy in places. I bet, only based on this derivation in prose plus formulas, nobody would be able to write a program to compute a VPM-B deco plan.

But we are lucky: There is a reference implementation in terms of a FORTRAN program. So, even though the VPM-B code in Subsurface is a complete rewrite (and solves many things algorithmically differently than the original code), we did a lot of testing to make sure that the plans we produce are identical to the plans computed using the FORTRAN program, even in places where we thought it makes little sense (like starting the ascent to the next stop when the ceiling depth equals the depth of the next stop rather than starting a bit earlier and making sure never to violate the ceiling during the ascent to the next stop (which is different since the ceiling goes up during the time it takes to get to the next stop). The latter is, at least to me, better motivated physically (“simply never violate the ceiling”) than basing the stop time on what happens at a different depth (namely the next stop depth) and we use the latter when computing Bühlmann schedules. But for VPM-B, we thought, given that we do not really understand the physics, we should not modify the model based on physics.

I have talked in the past about the Bühlmann model not being well defined to make all implementations come up with identical plans. In a sense, we are in a better situation with VPM-B, since there is the FORTRAN program which (at least for us) defines the model. But the problem is: This definition works only for the situation that you can compute with the FORTRAN program: You specify the bottom part of the dive and then let the program work out the ascent. Strictly speaking there is no definition for real dives: Dives where the distinction between bottom part and deco is blurred. But this is exactly what you have if you were to implement it in a dive computer or, as for Subsurface: Implement it in a dive log to show a ceiling for a logged dive (after all, Subsurface is mainly a divelog). For real dives, you cannot tell the exact point where bottom time ends and deco begins.

Unfortunately, the model depends on this: As explained in the VPM-B derivation post, one parameter that goes into the computation is the total deco time \(t_D\) (clearly that ends when you reach the surface, but where does it start?) but there are other parameters that you are supposed to evaluate at the beginning of the ascent like \(p_{1st ceiling}\) and the initial gradient for the Boyle compensation. All this depends on the time you call the end of bottom time and thus the computed ceilings also depend on that.

For Subsurface, we decided to take the point of time with the deepest ceiling (which, strictly speaking, is a circular definition but in practice is irrelevant) as the end of bottom time and base the ceiling computation in logbook mode on that. But that, to some degree, is arbitrary. And even for dives that we planned using VPM-B (and thus agree with the schedules computed by the FORTRAN program), applying this logic yields a slightly different ceiling. So it can appear that a VPM-B planned dive violates the VPM-B ceiling. But this is only due to the fact that the model is not really defined outside the planning situation.

Or put differently: You need to make arbitrary assumptions (not based on the depth profile and the gas you are breathing) to come up with a VPM-B ceiling. Which, at least for me, doesn’t strengthen my believe into this model.

PS: I am planning to be at the Boot dive show in Düsseldorf on Friday January 26th. If you happen to be there as well, please let me know and we can shake hands!

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.

In VPM-B you can shorten your deco by doing another dive before.  Sounds crazy? But this is what this model predicts as I recently learned from Subsurface user who thought he was reporting a bug. To be concrete, let’s plan an air dive to 40m with 30min bottom time (just for concreteness, maybe this is not the optimal gas choice but let’s keep things simple, the argument does not depend on this). With conservatism +2 we get the following plan
<
Runtime: 93min

But then add a 60m for 30m dive a bit less than two days (this is the limit that Subsurface considers two dives “repetitive diving”) before and you get:

Subsurface dive plan (surface interval 44:00)
Runtime: 78min

So, you get out of the water 15 minutes earlier! From the point of view of a Bühlmann-type model this sounds totally ridiculous: You would think that the effect of the previous dive is to add inert gas saturation to your tissues and indeed, even for VPM-B that is the case, but the consequences are opposite. The reason is that the previous dive sets the “maximal crushing pressure” and thus the “crushing radius” at the deepest point of the dive. And, as you know from a previous post, this enters the maximal gradient (at least the initial one) that VPM-B will allow you: Translated into gradient factors, without previous dive, Subsurface computes effective GFs as 24/71 versus as a second dive 26/92.

I was even more surprised that in the VPM-B literature this is considered a feature rather than a bug of the model: The previous dive pre-compresses the bubbles and thus they are less dangerous (or something like this). See the very last paragraph of this article.

Actually, to see this behaviour in Subsurface, you cannot use the latest version (or you have to comment out  the lines
for (ci = 0; ci < 16; ci++) {
max_n2_crushing_pressure[ci] = 0.0;
max_he_crushing_pressure[ci] = 0.0;
}


in the function clear_vpmb_state() in deco.c as I did)  as we decided this is too odd and rather opted for the more conservative option to reset the crushing pressure between dives. But that of course poses the question “What is exactly ‘between dives’?”. Is a short visit to the surface starting a second dive and thus magically resetting the state of bubbles in your body?

And of course the underlying question remains (until somebody does an empirical study): Can you really shorten the deco time by pre-compressing the diver?

It doesn’t even have to be a separate dive. According to VPM-B, I can bounce to 60m and immediately return back to 40m and then I get the benefit of shorter deco (compared to the version of the dive where I stayed at 40m):

Runtime: 79min