Maring correction for dust settling

[pcolarco]

The implementation of the Maring correction to the dust settling velocity is changed in GOCART2G from Legacy code. Previously the settling correction was applied as:
if(present(correctionMaring)) then
if ((correctionMaring) .and. (radiusInp(n) .le. (0.5*diameterMaring))) then
vsettle = max(1.0e-9, vsettle - v_upwardMaring)
endif
endif
where the correction is only applied to size bins less than a certain size (i.e., the diameterMaring = 7.3 um). In GOCART2G that second conditional is modified so the code reads as:
if(present(correctionMaring)) then
if (correctionMaring) then
vsettle = max(1.0e-9, vsettle - v_upwardMaring)
endif
endif
and the correction is applied to all dust bins regardless of size.

All GOCART2G runs appear to have DUSD001, DUSD002, DUSD003 resolving to negative values (which is unphysical and undesirable). The values are small and are almost certainly running up against the minimum allowed value (1e-9 m/s) fall velocity which results in essentially zero settling loss. This is not found in legacy GOCART runs.

It is not clear the change is causative, or logically why it would be. There is possibly some other numerical change that has happened.

1. Is the change in the application of the correction desired or just accidental (@adarmenov @amdasilva )?
2. What is the solution desired to avoid the incorrect negative fluxes (@dkim6git @pcolarco @amdasilva )?

[adarmenov]

I cannot comment on what is in G2G as I am not directly involved in this effort. That said the new G2G settling code was probably copied from my GOCART dev branch. I stand behind the implementation of the Maring correction in my branch.

Below are excerpts from Marring's paper that I can summarize in a single sentence -- they did not find evidence that the (unknown) mechanism counteracting the settling is size dependent:
[Section 3, paragraph 19] ...the shift from virtually no preferential loss to essentially total removal versus particle size appears more like a step function than the gradual transition predicted by Stokes settling. ...we hypothesized an upward velocity inside the air parcel partially counteracting Stokes settling. This upward velocity reduces the dust removal rate and maintains a relatively uniform dust aerosol size distribution in the vertical. We estimated the magnitude of the upward velocity to be equal to the Stokes settling velocity of the particles located at the point where preferential losses become observable (i.e., 7.3 mm diameter with Stokes settling velocity of 0.33 cm s-1). In Model 2 we hypothesized a vertical velocity for dust particles equal to the Stokes gravitational settling velocity minus 0.33 cm s-1. Calculated vertical velocities less than zero were assumed to be zero.
...
[Section 4, paragraph 21] ...A simple model with a vertical velocity for dust particles defined by the Stokes gravitational settling velocity minus 0.33 cm s-1 accurately reproduces observed losses of dust particles as a function of aerosol size.

In other words, codewise, the settling velocity correction should be applied to all dust bins. To some degree, this also makes sense physically if the hypothesized mechanism indeed acts similarly to an updraft - dust particles of all sizes should be affected by it.

It would be interesting to investigate the effect of relaxing the requirement for positive settling velocity, as this would allow for actual upward transport. I would not expect excessive accumulation of dust in the stratosphere (lower air density -> higher settling velocity).

I am not sure what could link the Maring correction to the negative DUSD001|2|3 observed in G2G. The small positive value in the vsettle = max(1.0e-9, vsettle - v_upwardMaring) expression should guarantee positive loss flux. Moreover, the 'older' and 'newer' codes should be working in the same manner for the small size bins. Why would only 1,2,3 be affected by this change? I would expect to see differences in the larger bins not the small bins.

Hope this clarifies the rationale for my changes to the Maring correction in GOCART.

[acollow]

I did not realize this discussion was here and stumbled upon the negative values for the first 3 bins from my own evaluation. I toggled the Maring flag on and off and can corroborate Pete's findings. I ended up putting in a conditional for my replay test run that the bin number must be greater than 3 but in hindsight that was not the best option for flexibility of the code. So, I am tempted to reinstate the radius conditional that Pete copied from the legacy code earlier in this thread. I am open to other suggestions if that doesn't seem like the best approach.

[pcolarco]

Allie, thank you for revisiting this, it evidently slipped off the radar. The dust code uses Chem_Settling function in the process library (as opposed to others). Could you test what an upward settling velocity does to the DUSDxxx diagnostics? I suggest right before the call to "SettlingSolver" just make vsettle a small (negative) number, e.g.,
vsettle = -0.001
which should equate to an upward velocity of ~0.1 km day, nothing dramatic, but meaningful. The diagnostics should be "zero" in value.

[acollow]

Yeah, I can do that test. Though I am going to continue my investigation of the settling solver first since I don't understand how column mass can increase from settling with a downward velocity.

[pcolarco]

By all means read the solver. Maybe there's an outright bug. Otherwise I suspect that numerical precision is the issue. You could increase that fall speed to be 1e-8 or 1e-7 and so on.