proposed fix for LSM negative W and sigmaS

[LenkaNovak]

Both snow, σS, and water, W, can become negative in the bucket, if their reserves are depleted and evaporation is nonzero. Simple test - set all fluxes to zero except evaporation. Evaporation should be attenuated by the beta_factor but small negative values of W and σS can still arise if

• fluxes are explicit and are calculated based on the state of the previous timestep
• when using CombinedStateFluxes, where ocean and land states from the same grid point are combined (this can/will be avoided when using PartitionedStateFluxes)

To ameliorate this, we could do one of the following :

• adapt beta - e.g. nonlinear decay to zero before zero W or σS are reached (in the Bucket)
• implement a limiter - e.g. something along the lines of the following (in the Bucket / Coupler)

        # crop W and σS for negative values
        @. dY.bucket.W = max(dY.bucket.W, 0)
        @. dY.bucket.σS = max(dY.bucket.σS, 0)

        # collect negative W and σS
        W_res = @. -min(dY.bucket.W, 0)
        σS_res = @. -min(dY.bucket.σS, 0)

        dY.bucket.Ws = @. (
            (
                liquid_precip + snow_melt -
                (1 - snow_cover_fraction) * evaporation
            ) - infiltration - W_res - σS_res
        ) # Equation (5) of the text, snow melt already multipied by snow_cover_fraction

The energy from sigmaS needs to be tracked for conservation (also check F_into_snow - this seems to be unused in the Bucket).

[juliasloan25]

Not sure if this will affect this fix, but the sign convention for precipitation/infiltration in ClimaLSM has changed as of v0.6.0 (see diff here)

[kmdeck]

@LenkaNovak I like the idea of adapting the beta function, but currently it is only used for W, and not for sublimation of snow. We could introduce a beta function for snow which only damps nonlinearly to zero, though, and tends to 1 otherwise. There could still be issues for large enough timesteps and large enough evaporation rates even with a beta that nonlinearly decays to zero with water content. But, the system would be mathematically well behaved in the limit of W->0.

A limiter would work too. The only complication there is that if we run out of water in the bucket during a timestep due to some evaporation, i.e. W(t) < E(t)*dt, we would have W(t+dt) = 0, and would need to then "condense" water from the atmosphere in the next step (since the atmosphere took the entire E(t)*dt the previous step, but the land did not have that much to give). So something like -(E(t)*dt - W(t))/dt would be added to the vapor flux atmos BC in the next step. I guess this would need to be tracked in the coupler. Im not sure we want the coupler to have to handle things like this, though.

Are the negative values causing things to break? Can we tolerate small negative values in Y since they conserve water exactly, and interpret them as zero water in reality?

[LenkaNovak]

(discussed offline to implement beta snow and check that small negative values can be tolerated within the bucket parameterizations, thank you @kmdeck ! )