Bugfix layer thickness

[rcnijzink]

Water table zw does not align properly with the surface level (i_cz). Solution is implemented (see https://github.com/schymans/VOM/pull/11), but we need to check.

[rcnijzink]

Comment Stan:

Thanks for fixing this, Remko! I just merge it, although I am still surprised that it makes a difference in the simulations, as I thought the bottom soil layer thickness was adjusted again elsewhere. How big is the difference?

[rcnijzink]

See attached:

[rcnijzink]

Comment Stan:

Thanks, Remko! This is strange, as I thought this affected only the initial conditions. Could you verify? There could be something else wrong as well.

[rcnijzink]

I think the initial conditions change slightly as well of course (will check), but the correction is only the thickness of the soil layers. It basically comes from this statement in transmodel.f90:

! * number of soil layers s_maxlayer assuming same thickness everywhere if (s_maxlayer .eq. 0) s_maxlayer = ceiling(i_cz / i_delz)

So the ceiling statements leads to a total soil layer thickness that can be slightly larger than i_cz, which happens by this:

    s_delz(:)   = i_delz

So now the lowest layer has a different thickness. Need to check maybe where generally the layer thickness (i_delz) is used, which should be now a loop or vector of s_delz.

[schymans]

Something is strange, it could be that the problem was fixed in the wrong place. In watbal.f90, the location of the water table is calculated in Lines 335-342, as zwnew. Could you check again why s_delz(s_maxlayer) needs a different value that the other layers, and how the error propagated into zwnew? The problem is that this one sets just one layer in the profile to a different thickness, which stays the same even if the water table moves away from there.

[rcnijzink]

I've looked a bit into it and tried to better understand what it actually does, but I'm actually not really sure what you mean.

So we have a vector s_delz which is completely static and stays the same during all timesteps. Condition here should be that the sum of s_delz equals i_cz minus i_zr. This was not the case, but now it is with the bugfix. So in the old case, the water table zw could reach higher levels then i_cz-i_zr. This happened in the lines you refer to:

zwnew = SUM(s_delz(wlayernew:s_maxlayer)) - 2.d0 & & * s_delz(wlayernew) * pcapnew(wlayernew) & & / (s_delz(wlayernew+1) + s_delz(wlayernew))

So, if pcapnew is not very big (i.e. we are close to (the supposed to be) surface), zwnew could be bigger then i_cz-i_zr, due to the use of sum(s_delz), now sum(s_delz) = i_cz-i_zr, so that cannot happen anymore.

The initial conditions are similarly defined:

do jj = wlayer_, 1, -1 pcap_(jj) = 0.5d0 * s_delz(jj+1) + 0.5d0 * s_delz(jj) + pcap_(jj+1) enddo

So here s_delz is used as a vector in a loop, so the correct values are used and pcap still increases linearly from the bottom layer.

Of course, we can also fix the problem in lines 335-342, by enforcing zw to be less or equal to i_cz-i_zr. But it all comes from the definition of s_delz and I think it's better to fix the problem at the root.

I will also check the code for using a fixed layer thickness (i_delz) or using all different values of the vector.

[schymans]

OK, here is how it should behave. The condition is that the water table (zw_) cannot be higher than i_cz, but it can indeed be higher than i_cz - izr, which is a pre-requisite for any seepage face flow. The initial condition is, however that zw = i_cz - izr. In your plot above, there are occasions where zw > i_cz, which I suspect can only happen if wlayernew = 1 and the second term in the zwnew equation above is negative. So I suspect that pcapnew has become negative there, which is not physically possible. It could happen if soil saturation degree becomes > 1. I don't know why the change of s_delz would have avoided this, as it only defines the static boundaries for the soil layers.

[rcnijzink]

Ok, I agree with your last point, with a negative second term we have values that are too high. I will go over the code to see if there are any checks preventing this and else I'll add it.

But in the occasion the second term is small and positive, we end up with a value just less than SUM(s_delz). When SUM(s_delz) is not correctly defined and in total more then i_cz, we also still may end up with values larger then i_cz. So I still support my first bugfix...

Just as an example: i_cz = 36.87 m i_delz = 0.5

This originally lead to 74 layers of 0.5m, with a sum of 37 m. So if the second term is less then 0.13m (37-36.87), we end up with water tables higher than i_cz. In cases we are close too saturation I think this is very likely to happen.

[schymans]

Ah, sorry, I did not realise that i_cz is the problem, not the initial conditions! It does not make sense to define i_cz with a centimeter precision, and in fact when I coded it, I exptected i_cz to be a multiple of the layer thickness, i.e. of 37 m in this case. You are right, to avoid such problems, we can make the bottom layer thinner than the others. We can still run into problems on the other end of the spectrum, when the water table vanishes but the bottom layer turns out to be e.g. 0.01 m. Such a thin layer would lead to quite a numerical burden. Perhaps it would be good to mention this in the documentation that i_cz and i_zr should ideally be a multiple of i_delz for smooth computations. OK, so there is no problem with the calculation of soil saturation degree or the position of the water table, just the fact that the soil profile is deeper than i_cz, but the calculations are based on the assumption that i_cz coincides with the position of the top layer. We can close this issue then.

[rcnijzink]

Ok, thats a good point, about the numerics. Let's just leave it open for a bit, cause I think there are two actions:

-add to documentation that i_cz must be a multiple of i_delz -maybe add an error/stop statement in the code when this is not the case, or a warning

[schymans]

Actually, we should have a general test that the stack of soil layers has a total thickness of i_cz and a layer transition at i_zr. If not, an error message should be raised. This will be particularly important if someone supplies a file with layer thicknesses.

[rcnijzink]

Closed with https://github.com/schymans/VOM/pull/21