The Theoretical Diver – Page 3 – Theorizing about scuba diving

Setting Gradient Factors based on published probability of DCS

This is a contributed post by Doug Fraedrich in which he reports on a recent paper by him and which is also mentioned in Doolette’s blog post I just commented on.

Recently I published a methodology on how to set conservatism factors on commercial dive computers based on published probability of decompression sickness (DCS) data and statistical summaries/models thereof, (Reference 1.) That paper described the general methodology on how to adjust the conservatism of several different algorithms; this note will focus specifically on setting gradient factor levels for the Bühlmann ZHL-16C algorithm.

The basic methodology is to use the output of probabilistic models from the literature that were derived from many well-documented experimental man-trials of known dive profiles. Then pick a probability of DCS isopleth curve and use those to set gradient factors (or conservatism factors) such that the algorithm will output the No-stop time or total decompression time (TDT) specified by that isopleth. This was done for three metrics: no-stop times for short recreational dives profiles, total decompression time for longer/deeper decompression dives and “first stop depth” for decompression dives. There were no probabilistic models for first stop depths/decompression profiles, so for this metric, we used data from a NEDU controlled experiment, which compared P_DCSoutcomes of dives of the same depth, bottom time and same TDT but with different stop profiles (Reference 2.) This NEDU study was performed specifically to assess the effect of deep versus shallow profiles as dictated by dual-phase vs tissue-loading algorithms.

So now to apply the methodology from that paper to the task of setting Gradient Factor levels for the Bühlmann ZHL-16C algorithm. From previous publications on the topic of gradient factors, we know that the value GF-Hi mainly affects total decompression time and GF-Lo affects the depth of the first stop. For GF-Hi, we will use the probabilistic model from a new study on No Stop Times by Howie et al. (Reference 3.) Using their criteria of 0.1 % serious DCS, we pick several point-pairs of depth and No-stop time and iteratively run theBühlmann ZHL-16C algorithm (specifically MultiDeco version 4.12)with different values of GF-Hi and match at what level the algorithm allows for a direct ascent at each individual depth-bottom time pair. The corresponding value of GF-Hi is 80.

As mentioned above, for GF-Lo, we use the US Navy NEDU study on deep stops, Reference 2. In the NEDU study, the maximum depth was 170 fsw, the bottom time was 30 minutes, the ascent rate was 30 fsw/min and the TDT for both profiles was 180 minutes. The group of divers who started their first stop at 40 fsw (for 9 minutes) had significantly lower P_DCS(P= 0.0489 one-sided, Fisher Exact Test) than the group that stopped first at 70 fsw (for 12 minutes). Since there is not sufficient information to know exactly where between the two tested depths is optimal, this test case used the two depths as a maximum and minimum for the first stop criterion. For that dive profile, MultiDeco was iteratively run varying GF-Lo to see which levels of GF-Lo yield the two limiting first stop depths. This procedure indicates that GF-Lo should be higher than 55 (this was previously reported in Reference 1.) If you pick the mid-point of the two first stop depths between those profiles (assuming the optimal point is somewhere in the middle and not at either extreme) the resulted in GF-Lo being on the order of 70.

Of course, the above setting for GF-Hi was based on no-decompression dives which had a relatively low value of PrT of ~10-12 (Pressure Root Time is an indicator of the severity of the dive exposure where P = pressure in bar, T = dive time in minutes, Reference 4.)We know from Reference 4 that historically, the incidence of DCS is significantly higher for dives with a PrT > 25 and we know from Reference 1 that the required value of GF-Hi for ZHL-16C needs to decrease as the PrT of the dives increases. So to set GF-Hi for higher PrT decompression dives, we re-visited methods and models presented by Van Liew and Flynn, where they were specifically assessing the suitability of the US Navy’s decompression schedule (used at that time) by fitting data on single-level, non-repetitive, nitrogen-oxygen dives from the US Navy Decompression Database to a logistic regression that resulted in P_DCSisopleths as a function of bottom time and TDT (Reference 5).We limited the domain of applicability of the current study to PRT< 40, so we used the “StandAir” Model which Van Liew andFlynn based ondata from standard air dives which had depths of less than 190 fsw and bottom times of less than 720 minutes. They assessed the StandAir statistical model to be reasonable except at the two depth extremes (nominally < 60 fsw and > 190 fsw). Van Liew and Flynn compared TDT required by the algorithm-under-test to the P_DCSisopleths from their statistical model, and found that the TDT required by the algorithm-under-test lay between the 2% and 3% P_DCSisopleths and thus deemed the algorithm acceptable for US Navy use. In the current study, we selected the 3% P_DCSisopleth as an initial standard of comparison as a compromise between managing DCS risk while not requiring excessive total decompression times. Note that in their analysis, they combined data from all dives that had the same depth and bottom time, thus their PDCS results are averaging over many different (non-optimal) decompression profiles.

The process here is similar to what was described above using the No Stop Data, only with the Deco dives, it was a trio of data values of: depth-bottom time-TDT (at the 3% PDCS isopleth curve.) MultiDeco was iteratively run to match TDT for the input values of Depth and bottom time. The value of GF-Hi was noted and plotted against the PrT of the dive profiles used. These results are summarized in Figure 1.

Note the curve in Figure 1 is of sigmoid shape: It is flat at both low and high levels of PrT with a transition in between. These results are based on dive profiles with observed DCS symptoms, and there are many more dives at the lower end of the PrT scale than the higher; so generally speaking the “error bars” start off small on the left side of the graph and get bigger as PrT increases. The point at which it levels out is considered somewhat uncertain.

To summarize:

GF-Lo is recommended to be >= 55
The recommended level of GF-Hi depends on the PrT of the dive, but is never greater than 80. It decreases with increasing PrT of the dive (to a point)
While the uncertainty at the low end (say PrT < 25) is low, it increases with PrT (i..e. more research is needed in this part of the domain.)

References:

Fraedrich DS, Validation of algorithms used in commercial off-the-shelf dive computers, Diving and Hyperbaric Medicine, V48 No 4 , 2018
Doolette DJ, Gerth WA, Gault KA. Redistribution of decompression stop time from shallow to deep stops increases incidence of decompression sickness in air decompression dives. NEDU TR 11-06, Panama City FL; 2011.
Howie, LE, Weber PW, Hada E, Vann RD, Denoble PJ, The probability and severity of decompression sickness, PLoS ONE 12(3): e0172665,2017
Balestra C. Dive computer use in recreational diving: Insights from the DAN-DSC database. In: Blogg SL, Lang MA, Møllerløken A, editors. Proceedings of validation of dive computers workshop. 2011 Aug 24 Gdansk. Trondheim: Norwegian University of Science and Technology; 2012. p. 99–102.
Van Liew HD, Flynn ET. A simple probabilistic model for standard air dives that is focused on total decompression time. Undersea Hyperb Med.2005;32:199–213.

Short comment on Doolette’s “GRADIENT FACTORS IN A POST-DEEP STOPS WORLD”

David Doolette has written a blog post “GRADIENT FACTORS IN A POST-DEEP STOPS WORLD” on the GUE blog that has caught some attention on the inter webs. He summarises the state of affairs about the current sentiment against too excessive deep stops, a trend that got momentum in the wake of the NEDU study.

It’s a nice summary but for those who followed the debates about this there are not too many news except in the last paragraph where he describes his personal take away: He says that newer results indicate that the increase in allowed over-pressures with depth (in the Bühlmann model expressed by the fact that the b-parameters are smaller than one and thus the M-value grows faster than the ambient pressure) is doubtful and that in fact the algorithms used by the Navy have the allowed over-pressure independent of depth, i.e. an effective b=1 for all compartments.

As with standard dive computers you are not at liberty to change the a and b parameters of your deco model. Rather the tuneable parameters are usually the gradient factors. So he proposes to set

\(GFlow = b GFhigh\)

citing 0.83 being an average b parameter amongst tissues and thus justifying his personal choice of gradient factors to be 70/85 saying “Although the algebra is not exact, this roughly counteracts the slope of the “b” values.”

This caught my attention as many divers lack a good motivation for some specific setting of their gradient factors (besides quoting that “they always felt good with these settings” and thus ignoring that this is as best subjective anecdotal evidence for a statement the would need a much higher number of tests under controlled circumstances to have any significance.

We could add this as a feature to Subsurface where you could turn on “Doolette’s rule” and then your GFlow would automatically be set as 0.83 of your GFhigh (we could even use the actual b value for the compartment and make the GFlow depend on the compartment). Of course, since we have access to the source code, we could directly set the all the b’s to 1 by hand and get rid of the gradient factors in that mode entirely. That might involve also modifying the a’s (as the now depth independent M-value) which made me worry that I might be pulling numbers out of thin air which when implemented might cause other divers to actually use them in their diving while being without any empirical testing, something I wouldn’t like to be responsible for.

So, rather being the theoretical diver, I fired off mathematica to better understand the “non-exact algebra” and to see if it could be improved. Turns out “non-exact” is somewhat of a euphemism.

Allowed over pressure: Green is 100/100 plain vanilla Bühlmann, blue 70/85 as of Doolette, orange is actually b=1 depth independent and red is ambient pressure. All numbers are bar.

Of course, how bad this approximation is depends on the depth at which GFlow applies (i.e. the first stop depth) but from the plot it is clear that a constant M-value corresponds to smaller and smaller GFlow (as a fraction between the red and green line).

Required GFlow/GFhigh for constant M-value as a function of ambient pressure (in bar) at first stop

As you can see, for somewhat greater fist stop depths (corresponding to deeper and/or longer dives), to keep a depth independent M-value one needs such a small GFlow that I would consider this to be very well in deep stop territory, opposed to what the blog post started out with.

If you want to play around with the numbers yourself, here is the mathematica notebook.

Equalizing real gases

As for filling cylinders real gas corrections do matter, we try to take the compressibility factor into account in Subsurface. As the theoretical physicist I thought (as I was taught) this is handled by the van der Waals equations. Turns out, in the pressure/temperature range relevant for diving cylinders is very far from getting the corrections to the ideal gas equations right. There are errors at least in the 30% range. This is probably since we are quite far away from the tri-critical point where it is supposed to be a good approximation.

In the end, we found a cubic approximation with coefficients read off from tables of compressibility factors to be an economic way to do these computations. You can find the source here.

When today somebody asked in a diving forum for a spreadsheet to calculate the resulting pressure when streaming between two cylinders, I thought I could use these interpolations to do the calculations for real gases with different compositions. Here is what I came up with.

Intro to Bühlmann/GF by H&W

In my posts here at the theoretical diver, I usually assume that besides some mathematics and physics literacy (which allows me to use formulas) my readers already know about the workings at least of the Bühlmann model and its generalisation using gradient factors.

Of course, that assumption is not always justified and for some time I wanted to write a foundational post on these topics (maybe with some personal twist, let’s see) but I never got around doing that. There is my old text (in German) but that covers only the very basics and in particular it is very short on gradient factors.

Some consolation was that there are a number of such descriptions available on the inter-webs. But now, on the occasion of the Boot trade show, there is a new one written by Ralph Lembcke und Matthias Heinrichs of Heinrichs Weikamp, the manufacturer of OSTC dive computers (and friends of the blog), Matthias and Ralph are now the principal developers of the OSTC firmware. It has a broader audience (and is in German as well) in mind and thus does a good job using diagrams instead of formulas. I like it a lot and give a strong reading recommendation (even though I might phrase the justification for using smaller GFlow differently, and there is of course my private pet peeve of not thinking of tissue half times as part of the model)!

Update September 2020: Now also available in English.

NDL and Gradient Factors

There is an interesting discussion over at ScubaBoard how the non-decompression limit (NDL) is affected by the settings for gradient factors (GFlow and GFhigh), in particular which of the two is relevant. My initial reaction was: It’s clearly GFlow, as that sets the depth of the first stop and the criterion for NDL is that the (theoretical) first stop is at a depth of 0.

Others argued that it must be GFhigh as that applies, by definition, at the surface.

And indeed, in this limit (of the first stop is at the surface) the idea of gradient factors degenerates: In that limit, the rate of change of effective gradient factor as a function of depth diverges. So, like in the last post, there is an interesting point that involves taking limits.

What I had in mind is the implementation in Subsurface: As you can see, as long as there has not been a stop yet (because the ceiling is still above the surface), the effective gradient factor is GFlow. So, in a no decompression dive, you would think, you never see anything else.

But to show that in the recreational mode in the planner turned out so be quite hard: I had to play a lot with the parameters to find a dive where GFlow has any influence on the total dive time of a recreational dive (defined as a dive without mandatory stops and without running out of gas). Eventually I found one: For an air dive to 20m (with an ascent rate of 20m/min for the last segment, see below why this is important), you get a maximal run time of 49min with GF settings 20/100 while you can stay for 50min with GF settings of 100/100. But changing GFhigh has much stronger influence:

What did I miss? In the end, it’s the fine-print of the definition of “first stop depth” that I already talked about in an earlier post: The problem is that for real world dives there is no clear distinction between ascent and stop. So you need to come up with some definition which depth one wants to use to actually anchor GFlow. Subsurface uses the lowest ceiling encountered during the dive so far. But in particular for dives with very little (or none at all) deco obligation that is not exactly what others might consider the first stop depth. The difference is that the diver first need to get to that depth of the ceiling before the ceiling actually becomes a stop depth. And during that ascent, there is already off-gassing going on which can eliminate the ceiling during the time it takes the diver to get there.

As an example, you could have a first ceiling (which as I explained above is determined by GFlow) at say 1m of depth. But then, in this last meter of water, the effective gradient factor has to vary from GFlow to GFhigh. Given that we are talking about dives that are only marginally deco dives, it is likely that this first ceiling comes form a very fast tissue so it is likely that much of it goes away during the short time of ascent to that depth. Then, to find the NDL, the remaining question is if there is ceiling left below the surface. But then the GFlow is already anchored at 1m so for the surface it’s really GFhigh (and GFlow is no longer relevant as there is no ceiling left at 1m where it applies).

So the challenge to find a dive where GFlow makes any difference at all for the NDL was to produce a dive where there is something left of the initial ceiling at the time when the diver gets there in the marginal case of staying a little bit shorter not occurring any stops at all. So the dive must not be too deep (otherwise the ascent takes too long and there is a lot of on the way off-gassing). That’s why I had to increase the ascent rate.

So the upshot is: It is almost entirely GFhigh that sets the difference between a non-stop dive and a decompression stop dive. But if you stay a little bit longer the depth of your first stop (and also the duration) depends a lot on GFlow.

I should not end without pointing out that once more this discussion is quite academic: Gradient factors were invented for dives that have significant deco obligation to force deeper stops. Here, we are in the limit of recreational no-stop diving. So we are really not in the realm of gradient factors. And this manifests itself in the degeneracy of the model in the case of the ceiling being exactly at the surface that determines the NDL. But it was interesting anyway.

A few thoughts on Oxygen Toxicity

Recently, at Subsurface, we decided to look into the calculations of oxygen toxicity related to diving. Here is what I learned (admittedly, not too much). According to what we could find on the inter-webs, there are two things to consider: The effect on your brain (the CNS value) and on you lungs (OTU for oxygen toxicity unit). The former seems to be simply a table look-up, so there is not much to do about it. The best we could do (thanks Willem) was to do a linear interpolation of table values.

For OTU, there is some thing to calculate. The empirical basis seems to be a measurement of the reduction of effective lung capacity after breathing O2 at higher partial pressure for an extended period of time. As far as we could tell, the current theory goes back the the thesis of Clark and his supervisor Lambertsen.

They found that there is no effect as long as the partial pressure of O2 stays below 0.5bar and any effect is proportional to the excess of partial pressure over 0.5bar and, according to their fit, it is also proportional to the exposure time to the power 6/5, i.e.

\((p-0.5bar) t^{6/5}\)

The power slightly bigger than 1 sounds somewhat believable since a lung that already suffers could be more susceptible to further damage. This value was obtained by plotting their data on a log-log-scale, fitting a straight line and reading off the slope of that straight line as the exponent thanks to

\(\log((p-0.5bar) t^{6/5}) = \log(p-0.5bar) + \frac 65 \log(t).\)

You can see this in their figure 18A on page 164. There you will also spot that the linear regression is not really good. This was also pointed out later in a more comprehensive study of the Naval Medical Research Institute that found that later data but also the original data of Clark and Lambertsen was more compatible with an exponent of 1. But the 6/5 stuck in the literature and all that follows below is non-trivial only because the power is different from 1.

What people then did was to take the 5/6-th root of this expression as the definition of the oxygen toxicity unit (normalised to a partial pressure of 1bar), i.e.

\(OTU = \left(\frac{p-0.5bar}{0.5bar}\right)^{5/6}t.\)

For time varying O2 partial pressure this could then be integrated over time, as for example proposed by Erik Baker (yes, the Erik Baker of VPM-B). Again, he publishes an implementation in his favourite programming language, FORTRAN. Well except, that he uses the reciprocal of that fraction and raises it to the power -5/6, with the additional benefit of risking a division by zero at the boundary case of partial pressure 0.5bar.

And of course, integrating this 5/6-th root makes no sense: You still get something linear in time whereas originally it was found that the damage progresses faster than time!

As you are integrating, you can also write a closed formula for a time segment where the partial pressure changes linearly in time (like during an open circuit dive during an ascent or descent). You need to compute

\(\frac 1{p_b-p_a}\int_{p_a}^{p_b} ((p-0.5bar)/0.5bar)^{5/6}dp = \frac 3{11(p_b-p_a)} \left.\left(\frac{p-0.5bar}{0.5bar}\right)^{11/6}\right|_{p_a}^{p_b}\)

This is Baker’s equation 2. Computationally, this formula has the disadvantage that a constant partial pressure is no longer a special case for this formula as one encounters a 0/0 floating point exception. Of course, taking a proper limit yields the above equation for the OTU but this is not convenient in a computer program.

So, what we did was to introduce better variables \(p_m=(p_a+p_b)/2\) and \(\delta = (p_b-p_a)/p_m\) in the integral expression above and then expand in powers of \(\delta\). By symmetry, only even powers appear and so a two term quadratic expression if good up to \(O(\delta^4)\), by far good enough for our purposes. This yields the improved expression

\(OTU = \left(\frac{p_m-0.5bar}{0.5bar}\right)^{5/6}t\left(1- \frac{5(p_b-p_a)^2}{216((p_m-0.5bar)/0.5bar)^2}\right)+O(\delta^4).\)

that can be calculated easily without treating the case \(p_a=p_b\) separately.

Once more, all this is very likely purely academic as it is not so easy to do dives that get into the regime where OTU matters at all. That is probably also why the empirical data is so poor.

Even more recently, there was a report on new results in this area on rebreather.org. Their study produced data summarised in these plots:

Obviously, this data clearly suggests to fit it by some convoluted formula that yields a line that I manually erased from the graphs. Guess what it is and then check out the original link.

Can we calculate no-fly-times?

There are various recommendations how long one should not fly after surfacing from a dive. DAN has recently done a study where they did Doppler measurements on a plane and recommends an interval of 12 hours after a single dive and 24 hours after repetitive diving.

But one would think that it should be possible to use a decompression model that works well under water to compute such a time. So let’s do this in this post or at least compute a conservative estimate. To be specific, we will use the Buhlmann model. We do that calculation compartment by compartment and assume that when leaving the water, the tissue has the maximally allowed partial pressure for surfacing (i.e. this tissue was the guiding tissue for the final part of the ascent). Clearly, this is a conservative assumption. Then, according to the model

\(p_s = (p_i -a)b.\)

Here, ps is the surface pressure and a and b are the usual Bühlmann coefficients.

Then we do a surface interval (whose length we wish to determine in the end) during which the partial pressure decays exponentially:

\(p_i(t) = f_{N_2} p_s+ (p_i(0) – f_{N_2}p_s)e^{-\gamma t}\)

where f is the N2 fraction the the breathing gas (0.79 for air which we are probably breathing while waiting for the plane to board). Finally, we don’t want the cabin pressure to be less than the minimal ambient pressure that Dr. Bühlmann recommends. We want to parametrize the cabin pressure using the barometric formula which asserts (assuming constant temperature) that the pressure drops exponentially at height with decay constant that the pressure is 1/e at about hs=8500m above sea level. So, we set

\(e^{-h/hs} p_s= (p_i(t) -a)b.\)

This we can then solve for t. Actually, being conservative, we want to throw in gradient factors. Again, being conservative, we don’t further linearly extrapolate gradient factors, but will use GFhi everywhere on the surface. Plugging everything in (with the help of Mathematica and some manual massaging) we find

\(e^{\gamma t} = \frac{a GF/p_s + 1-f_{N_2} + GF(1/b -1)}{a GF/p_s – f_{N_2} +e^{-h/h_s}(1+GF(1/b-1))}.\)

Now, we have to plug in numbers. For the cabin pressure, according to Wikipedia, we will assume h=2400m for older aircrafts. The wait times (in hours) for the different compartments are then shown in this plot:

No-fly-times in hours for different compartments plotted against gradient factor. The labels are the tissue half-times in minutes.

You are supposed to see a number of things:

The waiting time depends strongly on which tissue we are dealing with. For reasonably large gradient factors, only the tissues with half-times of several hours contribute significant waiting times. Remember, for this calculation, we assumed the loading is at its maximal value when you get out of the water. For realistic sports/tec-diving scenarios (as opposed to saturation diving) that should be quite hard even on a weeklong liveaboard with at least five dives a day. If slightly faster tissues are leading, the inferred no-fly times are much shorter, probably shorter than the queue at check-in. I looked at some data from real dive trips where people got everything out of their booked liveaboard but they got nowhere close to exciting the slow tissues. In the Subsurface planner I had to do five consecutive 2-hour dives with less than two hour surface interval to see at least some ceiling for the 239 minute compartment. In the plots, this has the blue line and leads to a no-fly time of much less than 5 hours.
The plot ends on the left (small gradient factors) with diverging values. These divergences move to higher gradient factors when you increase the cabin pressurisation equivalent altitude (for example by assuming a loss of cabin pressure, remember this is when the oxygen masks are supposed to drop from the ceiling)
Same plot with 4500m of altitude

This comes about as the no-fly time becomes infinite or even complex as according to the Bühlmann-with-gradient-factors limit, your are not allowed to be at the ambient pressure at that altitude even when saturated with the nitrogen that you experience at sea level. We can compute this limiting altitude by solving
\(f_{N_2} p_s = e^{-h/h_s} p_s (GF/b-GF +1) + a GF /p_s\)
for the altitude via
\(e^{-h/h_s} = \frac{f- aGF/p_s}{GF/b -GF +1}.\)
This is shown here, the maximum altitude after waiting an infinite amount of time:

Maximal altitudes in meters for the different compartments as a function of the gradient factor.
All these calculations are for air (or nitrox underwater since all we used was the assumption that the nitrogen saturation is at the limit). In particular, in view of an earlier post (N2 vs. He, what’s the difference?) there should not be large differences.
You could try to repeat the same argument for VPM-B but according to that model, if you followed it during the ascent, the no fly time would always be infinite: The ascent is determined such that when you surface and stay at ambient pressure, you will just create the maximum amount of free gas that is barely allowed. So going to any altitude and lowering the ambient pressure further would release more free gas than allowed, no matter how long you waited. The only way out would be to produce fewer bubbles on the earlier ascent while still in the water, then you would have some reserves to go to altitude.

What are we supposed to conclude from this? One takeaway message is waiting for the recommended 24 hours is not totally off, in particular if there is a chance that your very slow tissues have loaded a significant amount of nitrogen.

On the other hand, for realistic dives, 24 hours is likely on the very conservative side. At least from the perspective of decompression theory. From this perspective it is a total mystery to my what kind of reasoning dive computers use to determine a no-fly time.

Apart from these model considerations, DCS symptoms often enough do not show up immediately after a dive but up to several hours later. And when you are in that situation that you find yourself with DCS symptoms (even those that would have occurred irrespective of flying or not) your chances for immediate proper treatment are probably much higher if you are not confined to an aircraft above the middle of the ocean. So even from that perspective it makes sense to wait a bit more to make sure you will not need a chamber in the next few hours.

PS: If you want to play around with the formulas, here is my Mathematica notebook.

More confusion from Isobaric Counter Diffusion

After I expressed my worries about the commonly used criterium for the occurrence of Isobaric Counter Diffusion (ICD), namely a gas change in which the N2 fraction increases by more than a fifth of the drop in the He fraction, I learned (suggested by the OSTC source code, around lines 2220) about another (and much better motivated) criterium: Simply check if the leading tissue is taking up more N2 than it dumps He to the environment. That would mean that the net effect is an on-gassing of inert gases.

\(\Delta p_{He}<0,\quad \Delta p_{N_2}>0,\quad \Delta p_{He} + \Delta p_{N2} >0\)

This sounds reasonable: If this is true for the leading tissue, the decompression is ineffective as the total inert gas pressure goes up.

So I implemented it in Subsurface. Turns out, this does not trigger where you expect it. For example for a 60m dive with 20min of bottom time on 18/45, forcing a gas change to air (i.e. high N2) at 45m does not trigger it: Yes, He in the tissue goes down and N2 goes up but He is so much faster that the net effect is still off-gassing.

But for the same dive, it triggers at a different, unexpected place: At the beginning of the ascent (at 58m to 54m). How does this come about? At the end of the bottom time (and also during the start of the ascent), the leading tissue is the second (with 8min of N2 half-time). After 20min bottom time, He is almost saturated but N2 still has a way to go. Thus, at the start of the ascent, pretty much off the bottom, He starts off-gassing while for N2, the depth difference is not really noticeable, so it is still on-gassing with positive net. So there actually is counter diffusion even without a gas change!

I guess, nobody would suggest that leaving the bottom at 60m would be dangerous. But this seems to be the only place where counter diffusion actually happens!

I tend to believe more and more that this whole ICD story is either not explained at all by a diffusion model (maybe because it is only relevant in the inner ear that does not follow this simple tissue+environment model) or it is totally bogus.

So I would like to hear from you, dear readers, do you have any experiences with ICD or could you suggest a dive profile where it should be relevant?

Isobaric Counter Diffusion Criteria

I would like to hear your opinion about isobaric counter diffusion (ICD) and the criteria you apply to decide if a gas switch should be considered dangerous.

But first a bit of background to make sure we are all on the same page: ICD is the phenomenon that occurs when you switch from a mix with a lot of helium to a mix that contains less helium but more nitrogen. Then it can happen (depending on the tissue loading) that some tissue is off-gassing He but on-gassing N2, i.e. that the two inert gases move in opposite directions. And historically this is considered bad for your deco.

As the He atoms are much small than the N2 molecules the former diffuses faster than the latter so, in typical situations, in total the off-gassing should happen on a shorter timescale than the on-gassing and thus the net effect will be that the tissues inert gas loading goes down. Anyway, this effect should be covered by the usual diffusion based deco algorithm (as this is exactly what it simulates) and no additional care would need to be taken as long as the diver stays within the boundary of the deco algorithm.

There is, however, an argument due to Steve Burton that suggests to take the solubility of the gases (in typical tissue) into account to compute the absolute amount of inert gases in the body. And since that is about 5 times higher for N2 than for He, he argues that in order to have a net unloading of the absolute amount of gas one has to limit the change in N2 percentage in the breathing gas by 1/5 of the change in He percentage in the breathing gas. And at least one technical diving school of teaching seems to have adopted this criterium.

We are currently discussing if we implement this check into Subsurface and warn the user if a gas switch violates this “rule of fifth”. I am not sure though if I buy into this line of thought. After all, everything we do in deco planning is based on partial pressures, we never consider the absolute amount of gas. It is the partial pressure that plays the role of the fugacity determining if a particle moves across the diffusion boarder (in or out of the tissue say) and the rate is proportional to the differences (this is the mantra of diffusion based decompression models). So the solubility and with it the total amount of gas should play no direct role. I wonder if this rule of fifth (as is seems it comes out of a theoretical consideration with questionable initial assumptions) has received any empirical evaluation. Note that for example it forbids to change from Tx18/45 directly to EAN50 (as He goes down by 45% of which a fifth is 9% while the N2 fraction increases by 13 percentage points, more than 9%) even though I understand this is rather commonly done.

There is another ICD theory investigated by Doolette and Mitchell that focusses on inner ear DCS. They argue that in the inner ear the common diffusion model assumptions are violated since there is a relatively large amount of inner ear liquid that is not in direct (diffusive) contact with them ambient pressure (blood in practice) but only indirectly via inner ear tissue. So all off-gassing of the liquid goes first into the tissue. So what happens in an ICD caused inner ear DCS accident is that while that tissue is already on-gassing N2 from the environment (=blood) it is still receiving the He that comes out of the liquid and therefore experiences an over all uptake in inert gas which eventually causes DCS.

But i have not seen any hard “you should avoid doing the following” criterion that is derived from this line of thought beyond a general “shouldn’t really be a problem when using a CCR and otherwise be careful and always maximise the O2 fraction within the boundaries set by MOD”.

So, for the deeper trimix divers here: How do you decide which gas switches are ok and do you want your dive planning software to warn you about those.

PS: A last warning: It makes no sense to think about ICD in a way like “stuff moves in different directions, so the in-moving particles clog the out-going ones” as this little thought experiment shows: There are two stable nitrogen isotopes N-14 and N-15 (the former much more common in nature) that are chemically not distinguishable (only for example in a mass spectrometer). Image you are breathing a nitro mix with only N-14. Then you switch to the same mix but with all the N-14 replaced by N-15. Then there will be a counter diffusion of N-14 vs N-15. But of course, since both are chemically equivalent there are absolutely zero physiological consequences even, so the argument that the in-moving N-15 clogs up the out-moving N-14 cannot hold (and so it cannot hold in the case above).

VPM-B for real dives (or not)

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!