Apparent Water Loss in All-Land Experiments

[mcthreems]

Hello Isca team,

We recently realized that our experiments are not conserving water from their starting reservoir. We aren't sure if this behavior is known/expected, especially since we are running the bucket model with only land grid cells. We measured our starting reservoir by taking the starting bucket depth (same for every grid cell) and calculating the total volume/mass (accounting for change in grid area with latitude/longitude). We then compared this mass to the same calculation, but including the specific humidity at each level, for all subsequent years of output. The result is a significant decrease, over 80%, from the starting reservoir amount by the 25th year. The main confusion is where the water would go in the model. Previously, we had the issue where our maximum bucket depth was too small, and lost water when the buckets overflowed. But currently the maximum depth is set to be larger than the theoretical depth if all the starting water were to be placed in a single grid cell, so that shouldn't be where it disappears to. If anyone has any advice or suggestions on where to look in the code, it would be appreciated.

[gkvallis]

That's a very large water loss so it should presumably be diagnosable on a much shorter timescale, for example one year. The two places it could be lost are the buckets or in the atmosphere. Have you looked to see if the evaporation equals the precipitation when averaged globally and over time. If E > P then it looks like a sink in the atmosphere and if P > E then there is a sink in the ocean. That diagnostic requires a statistically steady state, so might not work in your case unless you can replenish the empty buckets with water from the full ones, but it might be a place to start.

Jamie, a student here, did some calculations like that quite some time ago to see where water ended up (the poles eventually), so will be curious to see your results when it is working.

[mcthreems]

Regarding the evaporation metric, we've been using the "flux_lhe" output, but wanted to clarify if that's the intended output for calculating the E for a P - E metric? Also, we don't appear to be in steady-state since the water is continuing to decrease over time, and we can't replenish buckets since the decrease is a decrease of the entire global reservoir rather than individual buckets. We would expect certain buckets to dry out over time and the water to move to other buckets, but instead we only see the drying in the tropics and a minor increase in depth at the poles.

[sit23]

In answer to your most recent query about flux_lhe this is the latent heat flux from the surface. So it's in Wm^-2, as here:
https://github.com/sit23/Isca/blob/master/src/atmos_spectral/driver/solo/mixed_layer.F90#L351. It's the flux of moisture multiplied by the latent heat of evaporation, as here:
https://github.com/sit23/Isca/blob/master/src/atmos_spectral/driver/solo/mixed_layer.F90#L742
To compare E and P you'll want them in the same units, so you could either get P into Wm^-2 or get flux_lhe into kg m^-2 s^-1. Up to you how you do that.
I've been doing some recent 'all buckets' simulations for Titan. I'll check if my total mass of vapour is conserved or not.

[gkvallis]

Once you have E and P properly computed, you should be able to calculate separate budgets for the atmospheric and bucket amounts, even when not in a statistically steady state. You can see if the increase, or loss, of water content in the atmosphere or bucket corresponds to the cumulative difference in E - P, and that might give you a sense of where the leak is occurring.

[RuthG]

Just to add, in case you've not seen, tendencies for the bucket depth can be output from idealised_moist_phys.F90, which might be useful in diagnosing what is happening. In my experience the atmospheric water budget is closed with a reasonably small residual, so I'd imagine the issue relates to the bucket, but under the conditions you describe I wouldn't expect a huge residual there either. Interested to hear what you find.

[jamesp]

If it's helpful, the code that I ran all-land buckets is here: https://github.com/jamesp/Isca/tree/jpbucket. Experiment file here: https://gist.github.com/jamesp/4b94941620cec9aa1c2546b6e8a536eb.

The difference from master is in idealised_moist_phys and sufrace_flux where I prevented buckets from overflowing with a namelist parameter. As far as I know this has not been merged into the master as I never finished the work of validating the results. It sounds like you may have already discounted overflowing buckets as an issue.

To make this non-overflowing change, I also had to change some of the logic around bucket saturation in the surface flux module. Evaporation is scaled by (bucket_depth/max_bucket_depth) so if you set the max bucket to be massive, you may have switched off evaporation. Have you considered this?

See master...jamesp:all_buckets to get the gist of the changes I made to run a non-overflowing bucket model. I don't recall water loss being a big problem (total bucket mass stayed constant over the decade timescales I looked at), but the atmosphere did dry out as all the water in the buckets eventually migrated to the poles from where it never returned to the atmosphere.

[mcthreems]

Thank you all for the suggestions. We'll look into the E-P budgets and the tendencies to see what we find. I'll get back to you with any new info.

[mcthreems]

I have looked into the E-P budget, and it looks like there is a consistent imbalance in favor of E up until the very end of the 25-year run. I have included a figure showing the cumulative difference over each day of the simulation.

[mcthreems]

Complicating things, I have found different values for water loss from the buckets based on the bucket_depth output itself and the tendencies. Using the base bucket_depth output, and calculating the total reservoir for each day of the experiment, I made a parallel figure to the one above showing the amount of water seemingly lost over time. This figure shows a loss of about one order of magnitude larger, as well as no leveling off at the end.

[mcthreems]

The bucket depth tendencies also seem to show a different amount of loss when adjusted to be kg/day and summed over time. The bucket_depth_lh tendency shows a much larger loss than the other two tendencies show in gain, as shown in the following figures.

[mcthreems]

Based on these figures, it appears that the E-P budget does not fully explain the loss seen in the bucket depth output (too small), but neither do the bucket depth tendencies (too large). I'm not sure if there is some other source/sink I'm missing, or if I've done the calculations wrong somewhere. Any further advice would be appreciated.

[sit23]

I'm afraid I'm not able to look at this in detail at the moment, but here are some things I would try:

1. Try changing the initial conditions - is all of your water initially at the surface or in the atmosphere? If it's all initially at the surface then you will see E win over P as the water evaporates, but it should eventually become steady. Try it the other way around (put all the water in the atmosphere instead) and see what's different.
2. Look at your water conservation options - in the spectral_dynamics_nml there is the option do_water_correction, and when set to true this will evaluate the total water amount before and after advection has been done, and adjust the total amount of water to make sure it's conserved. Have you got this option on or off? Also, there is the related namelist variable water_correction_limit, which sets the minimum pressure that water correction takes place over. If it's 0.0, then it will do this correction over the whole atmosphere. We often set it to 200hPa, because otherwise the water correction ends up placing a source of moisture in the stratosphere, which can be unhelpful in some experiments. If you've got it set to something non-zero, try making it zero.
3. If I were you I'd want to check that the dynamical core itself is conserving water OK. To do this cleanly, I'd turn off the bucket model, then turn off evaporation in the mixed-layer namelist, then turn off moist convection and turn off large scale condensation. That way any water in the atmosphere becomes totally passive, and you can check whether it is conserved when it's just being advected around. That's one way to rule things out.
4. It might be that you're having troubles in the boundary layer. I'd check how close to the ground your bottom pressure level is. Is it comfortably within the boundary layer? If it isn't, then I'd make sure there are a few levels down there to make sure things are resolved. If it is, then perhaps try using the neutral option for the boundary layer fluxes. That should simplify their calculation a bit, and so might reveal some complications you're having.

Hope that's enough ideas to get you started.

[mcthreems]

Regarding point 1 of your suggestions, the vast majority of our water starts at the surface, and we do expect a spin-up period where it evaporates out. Below is a figure showing the daily total water inventory of the atmosphere, calculated from specific humidity. The very early part of the run shows a quick increase from ~0, but this is done very quickly and the curve flattens out likely within the first year. In terms of magnitude it's very close to the loss calculated directly from the buckets, but that curve shows a steady loss over the entire run rather than a quick loss right at the beginning.

[mcthreems]

Regarding point 2, the do_water_correction option is set to true, while the water_correction_limit is set to 200 hPa. I can try setting water_correction_limit to 0, but might try some other stuff first since we also wouldn't want moisture in the stratosphere.

For point 3, I'll give this a try and get back to you on how it goes.

For point 4, our lowest pressure level is about 930 hPa, which I think would be in the boundary layer. The next level is about 800 hPa, which is probably above it. Would there be a way to add more levels near the ground without having to increase the entire vertical resolution? We can also just increase the vertical resolution generally, I'm just a bit concerned with increased runtimes.

[sit23]

OK - sounds sensible. Regarding your latest figure, is this consistent with your bucket water loss figures you sent earlier? From the look of the atmospheric water vapour total, you'd expect a big water loss initially from the buckets, and then it gradually creeping up from the first year, but that isn't what your bucket depth analysis shows, with its linear increase over time?

[mcthreems]

The fact that it is inconsistent with the bucket loss figure is likely part of the problem. The various figures I've made don't really match up with each other with regards to the shapes of the curves or the magnitudes, but I don't know why. My guess is water is being lost or going somewhere where it's not being included in these plots, but I need to figure out where that is.

[mcthreems]

I attempted to do a run as you described without any moist processes, but have been having trouble with errors. I'm not sure if there is a setting I need to adjust somewhere, but I turned off the buckets, turned off evaporation, and changed the convection to dry, which should also turn off the large-scale condensation. I decided to use the "initial_sphum" parameter to make sure there was water in the air, but when the model tried to start it crashed due to an overflow of the saturation vapor pressure table. I would have expected that table to not be in use, so I'm a bit confused what went wrong. The error message is copied below:

FATAL from PE 1: lookup_es_3d: saturation vapor pressure table overflow, nbad= 3840

[rosscastle]

Hi, that error is misleading, it usually is a symptom not a cause. The usual causes are unstable set ups (which I am afraid I can't help you with as this is not really my area, hopefully someone else can!). Was there any more error messages? You can try lowering the timestep as that can sometimes lead to things numerically blowing up and causing crazy values.

[sit23]

I agree with @rosscastle that this error message is often misleading. Once you've turned off those moist processes, the other places moisture is still used are :

1. In the calculation of relative humidity if you have asked for this as an output (easy thing is just to no longer ask for rh as an output)
2. In the surface flux calculation here: https://github.com/ExeClim/Isca/blob/master/src/coupler/surface_flux.F90#L416
   It's relatively straight forward to hack the latter calculation temporarily to stop this from happening. Just feed it artificial temperatures and sphum values so that it doesn't crash.

The alternative would be to use a very low value for intial_sphum. That way you should be within the lookup table's range of values.

In terms of the consistency of your graphs - I agree that they are inconsistent. In terms of interpretability, I think your total atmospheric moisture is the one that looks the most sensible, so I would believe that one. The bucket ones look suspiciously linear to me. I'm hoping to take a look at my all-bucket Titan runs today. Hopefully that will help in this investigation. Did you manage to look at some alternative initial condition cases?

[mcthreems]

Regarding the alternate initial conditions, I ran into the issue above when I attempted to increase the initial sphum to account for the same global bucket depths we normally have on the surface. I can try that again, but for a lower starting depth, and see if the model still crashes.

[sit23]

Hi @mcthreems - I've had a look at our Titan runs in terms of vapour conservation, and I've attached two related plots. The first atmospheric_vs_land_bucket_tracer_mass.pdf has the atmospheric and surface amount of methane in kg plotted over time (x axis units are Titan years). As you can see, we initially start this simulation with all the moisture in the atmosphere, and none at the surface, and the model gradually equilibrates. You'll see though, that the total vapour mass looks to be increasing. If you add them together and plot the total over time, you get this plot total_tracer_mass.pdf. This shows that, indeed, the total mass does increase slightly over time (by about 16% of the initial total). It is, however, to be noted that in this run I have set do_water_conservation:False. So it is to be expected that I may not be conserving vapour entirely here. Despite this, there isn't the kind of gradual vapour loss over time that you report, so this run in particular doesn't reproduce your result. If I were you, I'd think about trying a run with do_water_conservation:False, just to check there's nothing going wrong for you there.
The other thing to note about these runs are that I have a small amount of hyperdiffusion acting on the tracer field so smooth out some of the very fine scales (I've only added this because I had massive difficulty with small scales in the tracer field, and I noted that Tapio Schneider's moist Titan model has this hyperdiffusion). I doubt this is the problem you're having, but may be something you want to consider.
I would keep trying to see if you can get the model stable with some alternative initial conditions. It can really help with perspective! If you're still really struggling, please send us a minimum working example experiment script that shows your problem, and I'll try and take a look.

[mcthreems]

I have been unsuccessful in getting the passive water run working, but will continue to work on it. I did a run similar to what you described above where I tried to start with all water in the atmosphere using the "initial_sphum" nml parameter. I did this for an experiment intended to start with 1m of water in each grid cell, and made the initial_sphum value such that each grid should start with 1m of column water vapor. Below is the daily reservoir of the atmosphere and the surface, calculated the same as the plots I shared previously for other experiments. The atmosphere reservoir quickly drops from the start and eventually levels out about 2.5e18 kg below its initial value. The surface reservoir shows a quick increase after the start, but it only gets to a max of about 2e15 kg, a full 1e3 times lower than the observed drop in the atmospheric reservoir. It then drops off again back to near the starting value, without any corresponding increase in the atmosphere. So compared to your plots above these two curves seem pretty disconnected from each other.

[mcthreems]

Based on these two plots, it seems like water is leaving the atmosphere in the early part of the run without ending up in the buckets. What are the processes in the model that deal with removing water from the atmosphere? I know there is the convective and condensation precipitation, and both of those should connect to the buckets via their corresponding bucket depth tendencies. Is there another way for water to leave the atmosphere that I'm not accounting for here?

[sit23]

The latent heat flux of evaporation can sometimes be negative (i.e. can pull moisture from the air into the surface) but that doesn't happen very often. I think it's going to be difficult to make much more progress here without having a look at your namelist. Could you please send through your experiment script or the input.nml file for this run and I'll take a look.

[mcthreems]

This is the experiment file for the experiment with all starting water in the atmosphere, same as the one that produced the last two plots I shared above. It's written as a python script, but I made it a txt file to upload here.

run_rastest_bd1m_es01_rot1d.txt

[mcthreems]

I ran a comparison experiment equivalent to the one that produced the previous two plots, with the only difference being I turned off the do_water_correction option. The plots are pretty much identical though, which suggests that option isn't affecting much in my experiments.

[sit23]

Thanks for the experiment script - I'm now running some tests using some setups close to yours. Will let you know what I find.