Opened 3 years ago

Last modified 17 months ago

#2777 accepted project-task

Create global Van Genuchten soil ancillaries

Reported by: vidale Owned by: pmcguire
Component: JULES Keywords: Van Genuchten, Brooks and Corey, soil ancillaries, JULES, pedotransfer functions, PTF, ANTS
Cc: imtiaz.dharssi@…, ass98av, mtodt, omuller Platform:
UM Version:


Can you create global Van Genuchten soil ancillaries for JULES?

Here are the instructions:
Copy the Montzka files from
Read the VG theta_s, theta_r, n, alpha and Ks parameters.
Convert to theta_s-theta_r [m3/m3], 1/(n-1), 1/alpha [m] and Ks [mm/s] in Jules units.
Calculate theta_c and theta_w using VG equation.
Interpolate to 0.5 x 0.5 grid.
Add the thermal ancillaries and albedo from existing ancillary file.

Attachments (1)

UM_soil.pdf (795.8 KB) - added by pmcguire 3 years ago.
Should you dine at ROSETTA's, white paper presentation, by Imtiaz Dharssi

Download all attachments as: .zip

Change History (38)

comment:1 Changed 3 years ago by pmcguire

  • Status changed from new to accepted

comment:2 Changed 3 years ago by pmcguire

Here're the equations and code that I have working in Python2.7 now. It's just some essential excerpts from the code.
I will do the regridding with CDO and I will then merge the old soils WFDEI ancillary file (0.5 deg) with these new fields.
Can you look at the code below to see if I did things right or not?
Also, the original JULES WFDEI ancillary for soils does not have depth information (const_z = TRUE). But the Montzka database has the parameters at multiple depths.
Should I average these parameters over depth somehow? A weighted average?
Patrick McGuire?

vname_alpha = 'alpha_fit_0cm'
vname_n = 'n_fit_0cm'
vname_Ks = 'mean_Ks_0cm'
vname_theta_s = 'mean_theta_s_0cm'
vname_theta_r = 'mean_theta_r_0cm'
h_c = -300 #cm ; pressure head at field capacity, defined after Equation 3 of Montzka et al, 2017
h_w = -1500 #cm ; pressure head at wilting point, defined after Equation 3 of Montzka et al, 2017
fname = 'Hydraul_Param_SoilGrids_Schaap_0/'

data_alpha, lon, lat, nx, ny, fill_value = ReadData?(fname,vname_alpha)
data_n, lon, lat, nx, ny, fill_value = ReadData?(fname,vname_n)
data_Ks, lon, lat, nx, ny, fill_value = ReadData?(fname,vname_Ks)
data_theta_s, lon, lat, nx, ny, fill_value = ReadData?(fname,vname_theta_s)
data_theta_r, lon, lat, nx, ny, fill_value = ReadData?(fname,vname_theta_r)

data_m = data_n - 1.0/data_n # defined in between equations 1 & 2 of Montzka et al, 2017
data_SatE_c = (1.0 + (data_alpha*abs(h_c))data_n)(-data_m) # equation 1 of Montzka et al, 2017, with h = h_c
data_SatE_w = (1.0 + (data_alpha*abs(h_w))data_n)(-data_m) # equation 1 of Montzka et al, 2017, with h = h_w
data_theta_c = data_theta_r + (data_theta_s - data_theta_r) *data_SatE_c # inverting equation 1 of Montzka et al, 2017, with SatE = SatE_c
data_theta_w = data_theta_r + (data_theta_s - data_theta_r) *data_SatE_w # inverting equation 1 of Montzka et al, 2017, with SatE = SatE_w

theta_r = dataset.createVariable('theta_r', np.float32, ('lat','lon'))
theta_s = dataset.createVariable('sm_sat', np.float32, ('lat','lon'))
theta_w = dataset.createVariable('sm_wilt', np.float32, ('lat','lon'))
theta_c = dataset.createVariable('sm_crit', np.float32, ('lat','lon'))
satcon = dataset.createVariable('satcon', np.float32, ('lat','lon'))
sathh = dataset.createVariable('sathh', np.float32, ('lat','lon'))

theta_w[::-1,:] = data_theta_w
theta_c[::-1,:] = data_theta_c
theta_r[::-1,:] = data_theta_r
theta_s[::-1,:] = data_theta_s
satcon[::-1,:] = data_Ks / 100.0 / 86400.0 # convert from cm/day to m/s
sathh[::-1,:] = (1.0/data_alpha) * 100.0 # convert from cm to m

comment:3 Changed 3 years ago by pmcguire

The code looks OK. I only spotted this

satcon[::-1,:] = data_Ks / 100.0 / 86400.0 # convert from cm/day to m/s

Jules uses units of kg/m2/s which equate to mm/s. So the 100.0 should be 10.0.

Also you don’t have the 1/(n-1) field.

Jules doesn't store theta_r explicitly but instead uses theta_s-theta_r, theta_c-theta_r and theta_w-theta_r.

comment:4 Changed 3 years ago by pmcguire

How should I handle the layering in the Montzka data? Should I average over the layers somehow?
I don't recall seeing theta_r subtracted off the other thetas in the JULES documentation.

comment:5 Changed 3 years ago by pmcguire

I don't have a good answer on how to handle soil parameter layering. It is possible for Jules to use ancillaries with vertical variations. I will have to dig out my notes to see how to generate those. For now you might consider using the parameter values from a depth of around 50cm? PLV may be able to offer advice on this. Maybe we can try both ways and see how much sensitivity there is.


Changed 3 years ago by pmcguire

Should you dine at ROSETTA's, white paper presentation, by Imtiaz Dharssi

comment:6 Changed 3 years ago by pmcguire

Hi Patrick

The VG (van Genuchten) or BC (Brooks and Corey) parameters are usually determined using Pedotransfer functions (PTFs). There are many PTFs published in the literature and they give very different estimates. The third slide of the attached presentation (UM_soil.pdf) compares several PTF estimates of CH b and VG 1/(n-1) parameter. I've also shown uncertainty estimates from papers that provided them.

You can see that ROSETTA tends to produce lower values of VG 1/(n-1) [range: ~1 to 4] while Toth et al and Woesten et al tend to provide much higher values for VG 1/(n-1) [range: ~1 to 10].
Cosby et al provide values of CH b that are in reasonable agreement with VG 1/(n-1) from Toth et al and Woesten et al.
Carsel and Parrish tend to agree with ROSETTA except for Clay.
Rawls et al produce values that are generally higher than ROSETTA but lower than Cosby, Woesten and Toth.

My conclusion from this is that PTF derived parameters have very large uncertainties; typically factor of 2 or 3 for the CH b and VG 1/(n-1) parameters. Rosetta VG 1/(n-1) values seem too low while Cosby CH b values seem too high, especially for clay. I think that the Toth et al PTF might be the better choice than either Rosetta or Cosby.

Let me know if you agree or not.


comment:7 Changed 3 years ago by pmcguire

  • Cc a.verhoef@… added

comment:8 Changed 3 years ago by pmcguire

Hi Imtiaz

Thanks for sending this to me. I like the analysis and think you have done a great job. I have gone through some of the related journal papers.

I agree with your conclusion based upon the data you compiled about n,b, etc.

One question I have then is why you suggested that I use the Montzka VG n map to try to make new ancillaries? Montzka uses ROSETTAv1 but you think Toth VG n would be better.
Another question is: can I post this attachment on the job ticket website?


comment:9 Changed 3 years ago by pmcguire

Hi Patrick

You can post the presentation onto the ticket.

There is a lot of uncertainty and confusion about soil properties so I think that it is still worth testing the Montzka VG ancillaries. If the results are poor then we can eliminate Rosetta from future use.

I have a python script that can generate ancillaries using soilgrid and toth and would like to test those. If the results are good then perhaps we can make a case to PLV and Anne. Maybe even the UKMO will listen.


comment:10 Changed 3 years ago by pmcguire

Dear Imtiaz

Great! Thank you!

Would you want me to run your Python script and make the ancillaries? Or do you want to do that yourself?


comment:11 Changed 3 years ago by pmcguire

Hi Patrick

My ticket & code for the soilgrids/Toth is here:

I can send you some ancillaries if the code won't run on your machine.


comment:12 Changed 3 years ago by pmcguire

Hi Imtiaz

I had already run your ANTS soil ancillary code in December 2017, with your default settings as described in your ticket. It doesn't look like there have been changes to your ticket or to the code since then. Maybe ANTS has changed, however. I am going to rerun with the same setup as in December 2017, but since my anaconda path on JASMIN has changed since the group workspaces have been moved, I need to reinstall conda on CEDA JASMIN now, am I am doing that now.

Due to my reading in the past few weeks and due to my correspondence with you and with Carsten Montzka, I do understand things about soil ancillaries much more than I did in December 2017. I hadn't tried to run your ANTS soil ancillary code since then, and I hadn't tried to make the ancillaries for the WFDEI 0.5 deg grid or anything. But now, since I understand things better, then maybe I can make more headway.

One thing I found out from Carsten a day or two ago is that the Brooks and Corey soil ancillaries that he made for us (Pier Luigi's group in Reading) were made using the Rawls and Brakensiek (1985) PTF using the ROSETTAv2 software from Schaap. This ROSETTAv2 software apparently was capable of doing 16 different PTFs.

Have you been using the Toth ancillaries in your modelling studies?

Last edited 3 years ago by pmcguire (previous) (diff)

comment:13 Changed 3 years ago by pmcguire

Hi Patrick

I was sure that the ancillaries PLV sent me weren’t using the Rosetta PTF. Good to know that the Rawls and Brakensiek PTFs were used.

I haven't had much time to test any soil ancillaries. [redacted] I'm hoping that you can test a few different PTFs to see which ones work best with Jules.


comment:14 Changed 3 years ago by pmcguire

Hi Imtiaz

I adapted your Toth/ANTS/N320 python program to make VG ancillaries with all the 11 fields that the Rawls BC ancillaries have.
I tried running in JULES, but the pixels are not matching up right with the other ancillary files, so it doesn't work.

The biggest problem I think is that the source data I chose (UM ancillary from Monsoon:xcslc0:/projects/um1/ancil/atmos/n512e)
doesn't have a fractional landmask enabled, so all the coastal pixels have 0 land fraction, giving a different shape to the continents and islands. I think most of the other um1 ancillaries are this way. I chose n512 since it was at a higher resolution than the ERA 0.5 degree regular grid that we use for WFDEI driving. My hope was to regrid (at least initially) from n512 to ERA 0.5 degrees.

Do you know the path or where I can get the UM/PP versions of the n512 ancillaries (or the ERA 0.5 degree regular grid ancillaries) with the fractional land mask enabled?

Maybe my next step (without the above) will be to try to get your Toth/ANTS/N320 python code to read in my NETCDF ERA 0.5 degree regular grid files directly instead of reading in the UM/PP versions of the n320 or n512 gridded files.

comment:15 Changed 3 years ago by pmcguire

Hi Patrick

Does your machine have a directory "ancil/atmos/n512e/orca025"? If so, try using those ancillaries. I think climate models use the orca land/sea mask that contains many more islands.

If this doesn't work, please send me your land/sea mask and I will try to generate the ancillaries you require.

comment:16 Changed 3 years ago by pmcguire

Hi Imtiaz

Thank you very much for your suggestion! I tried the orca025 gridded UM/PP input ancillary files as you suggested, with your Toth script and the versions in there of the input ancillaries from orca. And it indeed helped a lot. There were more of the islands than before, for example, and the grid matched up a lot better. But I don't think it will be a good enough match for JULES for WFDEI. For example, the Britsh Isles have a different shape with the ORCA grid than with the WFDEI landmask with coastal pixels.

So I am trying right now some input UM/PP ancillaries I made by converting my WFDEI NETCDF ancillaries to UM/PP format with Xancil.
Maybe those will work with your script and maybe the resulting output NETCDF ancillaries [w]ill have the same landmask as my original input WFDEI NETCDF files.

I also investigated today why India and Mexico are so spatially homogeneous for example in their 'b' exponent. The TEXMHT USDA Texture classes appear to put most of India and most of Mexico in textural categories that are 'clay rich' categories with similar 'b' exponents. By 'clay rich', I mean 'Cl', 'ClLo?', and 'SaClLo?', for example. They're probably not necessarily clay-rich; that's just how I am calling them.

comment:17 Changed 3 years ago by pmcguire

Hi Patrick

Would you be able to send me copies of your WFDEI ancillaries. I can also try generating Toth/soilgrids ancils on the WFDEI grid.

I agree that Mexico and India don't show much variability. This seems to be due to Soilgrids, Fig 5 of Montzka shows the same thing for VG n and alpha; .


comment:18 Changed 3 years ago by pmcguire

Hi Imtiaz

I think I have some WFDEI soil ancillary files that are almost there with the Toth PTF with ANTS. I might try them with a JULES gridded run tomorrow (Tuesday) (at least over a small grid like the UK, if not for a global grid).

I do note that some of the islands and a few of the coastal pixels are different in the resulting files than in the original WFDEI NETCDF files. I traced this down to the fact that the TEXMHT_M_sl1_5km_ll.tif file doesn't show the same islands that are present in the original WFDEI landmask. But most of the coastal pixels islands are the same as before.
If you know a good way to solve this problem, please let me know.

If you want to try to repeat or inspect what I have done, on CEDA JASMIN, I have put some files in
The files in the nancil_test_pljules-vn4.4_pcmTest1 subdirectory are the input WFDEI NETCDF ancillary files (from Pier Luigi Vidale from Carsten Montzka, using Brooks and Corey parameters from the Rawls & Brakensiek 1985 PTF (implemented in ROSETTAv2)) and the output UM ancillaries from Xancil on Monsoon.
The soil_bulk_density variable is messed up in the output file now, but maybe that variable is not needed.
I ran the script run5b.scr on CEDA JASMIN (bsub < run5b.scr), which calls a modified version of your Python program
to make .
Then I use cdo to put this on the WFDEI grid:
cdo -v selindexbox,361,360,1,360


comment:19 Changed 3 years ago by pmcguire

The latest status is that I now want to figure out how to get the right input TIF grid of USDA soil textures, so that when we combine this with the Land/Sea? Mask from the old WFDEI input ancillaries, there are no missing islands or missing coastal pixels in the new WFDEI output ancillaries.

comment:20 Changed 3 years ago by pmcguire

  • Type changed from task to project-task

comment:21 Changed 3 years ago by pmcguire

  • Cc lecSS13 added; a.verhoef@… removed

comment:22 Changed 3 years ago by pmcguire

  • Cc ass98av added; lecSS13 removed

comment:23 Changed 3 years ago by pmcguire

Hi Patrick

I've slightly modified the Python code as shown here to improve the handling of missing data. This change should allow all the small islands that were missing to appear.

I'm testing the change now, if you also want to test the change then please go ahead and send feedback.


comment:24 Changed 3 years ago by pmcguire

Thanks for the update to the code, Imtiaz!
I am running the Python code now.
Will defining the missing islands as NODATA really fix things?
I will try it with JULES soon. Do you think it will also work for the UM?

comment:25 Changed 3 years ago by pmcguire

Hi Imtiaz

I ran the Python script with your modification. And the output doesn't have the missing islands!
Thank you! I am trying to run it now in JULES for the UK only,


comment:26 Changed 3 years ago by pmcguire

Hi Imtiaz
I have been able to get a VG TOTH/PTF soil ancillary into JULES, when running over the UK only on the WFDEI grid. I am using the one with the bugfix from you […]. But I still have a mismatch to sort out between the number of soil points and the number of land points.

I had subdivided the UK over 8 processors (mostly to check the parallel processing part of the code). And for each processor, I get some errors like:

{MPI Task 0} [FATAL ERROR] init_ic: All points should be either soil or land ice points - have land_pts = 71 and soil_pts + lice_pts = 65
{MPI Task 5} [FATAL ERROR] init_ic: All points should be either soil or land ice points - have land_pts = 29 and soil_pts + lice_pts = 23
{MPI Task 2} [FATAL ERROR] init_ic: All points should be either soil or land ice points - have land_pts = 98 and soil_pts + lice_pts = 87
{MPI Task 1} [WARNING] init_ic: Provided variable 'rgrain' is not required, so will be ignored
{MPI Task 1} [FATAL ERROR] init_ic: All points should be either soil or land ice points - have land_pts = 88 and soil_pts + lice_pts = 82
{MPI Task 3} [WARNING] init_ic: Provided variable 'rgrain' is not required, so will be ignored
{MPI Task 3} [FATAL ERROR] init_ic: All points should be either soil or land ice points - have land_pts = 96 and soil_pts + lice_pts = 89

I will keep at it.

comment:27 Changed 3 years ago by pmcguire

Hi Patrick

I was setting sea points to np.nan which might not be interpreted correctly. I've modified the code so that sea points are set to missing. The changeset is here

Try generating ancils using this new code.


comment:28 Changed 3 years ago by pmcguire

Hi Imtiaz
I am trying that change now. But so far, it seems that JULES freezes or something in the first spinup year in the UK. It doesn't freeze and it gets to other spinup years OK with the old WFDEI soils ancillary. I don't know if I have any error messages, since the log files don't get written until the end of the run. And there is no explicit failure yet, so the end of the run hasn't happened.

I note that one difference between the old soil files and the new ones is that for the new ones,
the sea points are 1e20 for the recalculated layers and -1073741824. for the layers of the soil ancillary
that were just pulled from the old soil ancillary.

So I just tried to make a new ancillary with[isea] = -1073741824.
for the new layers. But that doesn't seem to work properly. In ncview, the NAN is now messed
up for the different layers,

Maybe I will try to set all the layers including the ones that are not recomputed to have 1e20 for the mask value.


comment:29 Changed 3 years ago by pmcguire

Hi Imtiaz
I didn't try the convsh. But thanks for the suggestions.
I did try today to modify further your original Python code instead.
In addition to your suggestion of[isea] = MASKED_VALUE
I added another line that does this for the variables that I modified:
MASKED_VALUE = -1073741824.

The output NETCDF soil ancillary looked fine for me in ncdump and ncview, with masked values all the same for all the variables.

Then I tried to run it over the UK with JULES.
In the file src/initialisation/standalone/initial_conditions/
It's failing now this condition on sathh for a small number of the land_pts pixels in the UK grid region:

 IF ( ANY( sathh_soilt(:,:,:) < 0.0 ) ) THEN 

I put some print statements, and it looks like those 'small number of pixels that fail' have the mask value, which is large and negative.
Most of the land_pts are ok, it looks like.

None of the land_pts should be masked.

I have to sort our why this is true.
I'll keep at it.


comment:30 Changed 3 years ago by pmcguire

Hi Patrick

It took a lot of tries before I managed to get me code to use a different missing data indicator. In the end, I only needed to add a line, -32768.0*32768.0)

Your attempt was pretty close to what I ended up doing.

The changeset is here

The code is here


comment:31 Changed 3 years ago by pmcguire

Hi Imtiaz
Thanks for checking the landmask vs. the output soil ancillary!

JULES is running now for the UK with the new soil ancillary. It made it through spin up pretty quickly
and now it's doing the main run from 1979-2012 ([four] years of the main run so far).

The problem I had last night with a few bad pixels was caused by a shift by half a pixel in the longitude grid which I hadn't noticed before. It is fixed now. The half a pixel shift was caused by the postprocessing I was doing on the output NETCDF soil ancillary. After running the Python program, I was doing two cdo commands to get the right format for the NETCDF file for JULES.
The first cdo command was the problem since it changed the grid, and I have replaced the first cdo command with this command
ncrename -h -d latitude,lat -v latitude,lat -d longitude,lon -v longitude,lon
The second cdo command that I am currently using remaps from the UM grid to the JULES grid, shifting the longitude by 180 degrees for this half degree grid:
cdo -v selindexbox,361,360,1,360


comment:32 Changed 3 years ago by pmcguire

I have put the code I have developed from Imtiaz's code in the MOSRS repository:

This code is used by calling scripts to generate from a choice of Toth et al.'s 3 different PTFs.
The latest calling scripts are at:

The 3 PTF calling scripts compute on CEDA JASMIN the soil ancillaries for all 7 soil layers, and they call this script for a single layer:

The Conda ANTS environment needs to be activated prior to running the 3 PTF calling scripts.

I have computed the hydraulic conductivity at the wilting point for the whole globe on a WFDEI grid for each of the three Toth (Van Genuchten) PTFs. I am doing this to have some idea of the minimum value of hydraulic conductivity. Once I have done the same for the Rawls&Brakensiek (Brooks&Corey) PTF, I will send details of the comparison of the 4 PTFs' hydraulic conductivity at the wilting point.

Last edited 3 years ago by pmcguire (previous) (diff)

comment:33 Changed 3 years ago by pmcguire

Here is the poster we presented on the soil ancillaries at the NCAS staff conference earlier this week.
There is also a slides version of the poster.

PDF format:
PPTX format:
slides version: PPTX format:

comment:34 Changed 2 years ago by pmcguire

  • Cc mtodt, omuller added
  • Keywords JULES, pedotransfer functions, PTF, ANTS added; JULES removed

Hi everyone
It looks like we made a decision during the wider University of Reading Land Surface Processes (LSP) group meeting (also with Martin) today to go ahead and try running JULES in global offline mode with the Tóth et al soil-physical-property pedotransfer functions (PTFs) 17and20. The PTF17 is for the Ksat saturated hydraulic conductivity, and the PTF20 is for the b=1/(n-1) exponent, the saturated soil moisture sm_sat, and the saturated psi sathh. The second choice was the (older) OLAM-Vereecken option with Soilgrids using the Weynants et al 2009 PTF.

Martin asked about the ANTS-based code for generating the ancillaries with Toth 17+20. The slightly out-of-date code for this is listed above in this ticket, in the comment 32, where there are links to my branch of Imtiaz's branch of the MO ANTS Python code, as well as to links in my branch to scripts that call this Python code. I will try to check in my more recent changes to this code later today if I can.

After the talk during the LSP group meeting today, I edited the version of the presentation a bit further, getting rid of some typos, and improving the sub-figure arrangement for Ksat a bit. I am putting links to this at:
It is version '*v8b_VideoWall1c' and '*v8b_VideoWall2c'. The '*v8b_VideoWall1c' file is the normal talk, but is has the option halfway through to switch to the '*v8b_VideoWall2c' version of the talk that is better-optimized for the video wall. If you have a small screen, it might be better not to use that option.

I will try to start running the offline JULES WFDEI-0.5degrees runs with the Toth17and20 PTF today if I can.

comment:35 Changed 2 years ago by pmcguire

I have started runs with the suite u-bb422UKTQ for the UK region and u-bb422CNTQ for the S+E Asia region on JASMIN now, with the latest version of the TothContin17and20 soil pedotransfer function (PTF). This is using WFDEI driving data from 1979-2012 and the Medlyn branch from JULES 5.2.
If these work without crashing, and if they produce reasonable data output, I will start global runs and other runs with JULES 5.2 trunk and JULES 5.6 trunk, etc.

comment:36 Changed 2 years ago by pmcguire

The runs with the suite u-bb422CNTQ over the S+E Asia region didn't work. They had some missing islands in the ocean off the SE coast of Asia.

I fixed those for now by defining NODATA values for the soil parameters, to correspond to the TothDiscretePTF19 Silt Clay Loam discrete-texture class. The new suite for the S+E Asia region (u-bb422CNTS) ran to completion over the main run period of 1979-2012.
I have also run to completion a global run with 35 years of spin up u-bb422GLTS. I am trying to evaluate these runs in ILAMB now.

These NODATA values need to be updated to correspond perhaps to the b=1/(n-1) exponent, the Ksat, and the theta_sat computed from TothContin17+20 PTF for the global averages of Sand,Silt,Clay,pH, Cation-Exchange Capacity (CEC), and Organic Carbon.

comment:37 Changed 17 months ago by pmcguire

  • Cc p.l.vidale@… removed
  • Reporter changed from pmcguire to vidale
Note: See TracTickets for help on using tickets.