Evaporation does not affect soil concentration

[dfyx]

Note: this report is mostly based on the blog post. Might be that the code already does something that I missed

The way you have implemented evaporation applies it as a constant factor to the overall volume of your drop. Shouldn't it only affect the liquid parts, leaving the sediment unaffected while also updating the sediment concentration? Something like this:

float sedimentVolume = drop.volume * drop.sediment;

// Only evaporate liquid
float newLiquidValue = drop.volume - sedimentVolume * (1.0-dt*evapRate);
drop.volume = sedimentVolume + newLiquidVolume;

// Calculate new sediment percentage
drop.sediment = sedimentVolume / drop.volume;

Otherwise you'd effectively evaporate sediment.

[weigert]

That is a correct observation.

If you think of the sediment value as a dissolved concentration instead of a volume in its own right (a suspension), then this is easier to handle.

The volume is lowered by the evaporation rate factor:

volume = (1.0 - dt evapRate)

The sediment would be concentrated accordingly:

sediment = (sediment volume) / (volume (1.0 - dt*evapRate))

or

sediment /= (1.0 - dt*evapRate)

This makes the system mass conservative.

I have added this minor change, but it made little difference. This is because evapRate has a very small value. For larger timesteps or evaporation rates, it would have a more significant impact.

Good catch!

[dfyx]

Your approach still does something different from mine. Imagine three drops with the following parameters:

Name	volume	sediment
A	1.0	0.5
B	2.0	0.75
C	1.0	0.75

A and B have the same amount of liquid while A and C have the same overall volume. Taking a high dt * evapRate of 0.1 to illustrate the difference, we get the follwing values.

Name	Your volume	Your sediment	My volume	My sediment
A	0.9	0.56	0.95	0.53
B	1.8	0.83	1.95	0.77
C	0.9	0.83	0.98	0.77

In your model, the amount that gets evaporated depends on the total volume while in mine it depends on the liquid volume. In extreme cases your model could still lead to evaporation of sediment and a sediment percentage of over 100%. I'm not sure if this can ever happen within realistic parameters but it's something that could eventually lead to serious problems.

[dfyx]

Also note that your death condition for particles might be affected by this. A particle that's close to 100% sediment probably won't move in the same way that a suspension does.

[weigert]

We are both modeling a mass conservative system, but making different assumptions in the beginning.

I say that mass is conserved, where sediment represents a concentration (or mass fraction), you say that mass is conserved where sediment represents a volume fraction. That makes the assumption that the sediment has its own volume (which is a fine assumption to make). It's the same deal though for mass conservation, and makes no effective difference. I do not evaporate sediment.

If the volume approaches zero for my system and sediment mass is conserved, then of course the sediment value will diverge to infinity. This makes sense for concentrations and for mass / volume fractions. Besides, what is the mass or volume fraction of sediment in a control volume of zero volume and finite mass?

Both our our code implementations conserve mass and are correct, just in a different way.

N represents the current time-step.

My case:

V_TOT(N) = V_L(N) V_L(N+1) = V_L(N)*(1-dtF) V_TOT(N+1) = V_TOT(N)*(1-dtF) M(N) = V_TOT(N)*X(N) //is conserved M(N+1) = M(N) = V_TOT(N+1)*X(N+1) X(N+1) = X(N)/(1-dtF)

Your case from your code posted above:

V_TOT(N) = V_S(N) + V_L(N) V_S(N) = V_TOT(N)*X(N) //is conserved V_L(N+1) = V_L(N)*(1-dtF) = (V_TOT(N) - V_S(N))*(1-dtF) //I assume you forgot a bracket above V_L(N+1) = V_TOT(N)*(1 - X(N))*(1-dtF) V_S(N+1) = V_S(N) V_TOT(N+1) = V_S(N+1) + V_L(N+1) = V_TOT(N)*X(N) + V_TOT(N)*(1-X(N))*(1-dtF) V_TOT(N+1) = V_TOT(N)*(1 + X(N) - X(N) -dtF + X(N)dtF) V_TOT(N+1) = V_TOT(N)*(1 + X(N)dtF - dtF)

therefore

X(N+1) = X(N) / (1-dtF+X(N)dtF))

Conserving the volume is entirely valid for a constant density.

Note that for this case, as the evaporation rate dtF approaches 1.0, the volume fraction approaches one. This physically makes sense when thinking about volume fractions.

Here is your calculation again: If dtF = 0.1, for an initial sediment value X (concentration for me, volume fraction for you) of 0.5:

A: Me: X(N+1) = X(N) 1 / (1-dtF) = 0.56 You: X(N+1) = X(N)(1-dtF+X(N)dtF) = 0.53

B: Me: 0.75 -> 0.83 You: 0.75 -> 0.77

and so on.

These formulas correspond to exactly what you put in the table, so you know my derivations are correct. Basically you can see the sediment value in my code represents something different than you think.

The mass transfer in my method is based on a concentration difference. This works, because my expression for the equilibrium concentration can have values outside of the range [0, 1]. Your mass transfer method, if applied directly to the sediment volume fraction, would have to be expressed in volume terms, and be bounded in that range. How do you do this?

In the end it is the exact same thing. Both systems conserve mass, it is just represented differently. Thanks for the discussion though!

Edit: I made a math mistake. Our implementations are not identical. My other points still hold.