Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Possible typo in topsoil function #13

Open
blimlim opened this issue Oct 12, 2023 · 11 comments
Open

Possible typo in topsoil function #13

blimlim opened this issue Oct 12, 2023 · 11 comments

Comments

@blimlim
Copy link

blimlim commented Oct 12, 2023

Hi everyone,

I'm hoping to replicate the calculation of the top 10cm soil moisture mrsos that APP4 does for some runs of ESM1.5, which I believe is done using the topsoil function in app_functions.py:

def topsoil(var):
    return var[:,0,:,:]+var[:,1,:,:]+var[:,2,:,:]*.012987

Based on the soil depths of 0.022, 0.08, 0.234, 0.643, 1.728, and 4.6m from the cableSoilLevels function, I'm wondering whether the .012987 should instead be 0.12987?

Many thanks,
Spencer
@paolap
Copy link
Collaborator

paolap commented Oct 12, 2023

I'm wondering if these depth levels might be different from what the model outputs. Based on what you say and assuming these values are defined at the centre of the level, then you wouldn't even need to use the 3rd level at all.
Does your output has the depth_bnds?

@blimlim
Copy link
Author

blimlim commented Oct 12, 2023

I've been basing the depths on this function in APP4:

#level values and bounds for soil depths in cable
def cableSoilLevels():
    levels=np.array([ 0.011 ,  0.051 ,  0.157 ,  0.4385,  1.1855,  2.872 ],dtype=np.float32)
    boundmin=np.array([0,0.022,  0.08 ,  0.234,  0.643,  1.728],dtype=np.float32)
    boundmax=np.array([ 0.022,  0.08 ,  0.234,  0.643,  1.728,  4.6  ],dtype=np.float32)
    bounds=np.column_stack((boundmin,boundmax))
    return levels,bounds

I might not have interpreted it correctly, but my understanding from this is that 0.022, 0.08, and 0.234 would be the bottom of the first three levels. Good question about what depths the model is actually using though. All I've been able to find in the output is a soil_model_level_number dimension which is just 1,2,3,4,5,6, but no depth bounds variables.

I have heard some different values for the soil depths elsewhere (access hive and Kowalczyk et. al. 2013). I'd assumed that possibly APP4 was using a more recent version of the depths, but it would be good to confirm this.

@paolap
Copy link
Collaborator

paolap commented Oct 12, 2023

So if those values are the bottom of the last level you are correct that should be 0.1298.. as it's 2 cm in the first 10 cm divided by the 15.4 cm total thickness of the layer

@paolap paolap mentioned this issue Oct 12, 2023
Closed
@chloemackallah
Copy link
Member

Hi Spencer,
You are right, this calculation is wrong!
It would be good to get an idea of the difference it makes, I'll need to think about whether we should put in an errata as an alert, or fully retract the data.
Chloe

@blimlim
Copy link
Author

blimlim commented Oct 12, 2023

Thanks Chloe and Paola,

If it would be useful I could run the two versions of the calculation on the runs I'm working with (30 members of AMIP style runs) to see what the errors look like?

@chloemackallah
Copy link
Member

Yes that would be fantastic!

@paolap paolap closed this as completed Oct 12, 2023
@paolap paolap reopened this Oct 12, 2023
@blimlim
Copy link
Author

blimlim commented Oct 12, 2023

I'll try and upload the results early next week!

@blimlim
Copy link
Author

blimlim commented Oct 16, 2023

Hi Chloe,

Just uploading some plots of the mrsos errors from the model runs. For some context on the runs – I'd run 3 AMIP style ensembles from 1980 - 2014, each with 10 members. One used observed SSTs and the other two used artificial SSTs, but for the following plots I've bunched the three ensembles together.

The following shows mrsos_true - mrsos_false averaged over the 30 runs and also over 1980 - 2014:

The following shows the maximum monthly error occurring in each grid cell over the 30 members and the run duration:

I've also looked at the average percentage error (100(mrsos_true - mrsos_false)/mrsos_true for the following plot):

The average percentage error is often around ~20%, though it does vary spatially.

I wanted to get an idea of whether the percentage error varied over time and between the runs, or instead if it stayed pretty constant. The following shows the standard deviation of the percentage error at each grid cell taken over time and the ensemble members:

The variation looks quite low in higher latitude NH, but larger elsewhere.

Finally here are just a couple of timeseries at different gridpoints from a single run:

I'm wondering whether it would be good to look at whether linear trends are affected much or not. I can also share the nc files on gadi if that would be useful!

Best,
Spencer

@chloemackallah
Copy link
Member

Hi Spencer,
These are great, thank you! So it looks like we have a rather significant issue here.

At this point, I'm thinking it would be easiest to just retract mrsos (the only variable that uses topsoil()), and issue an errata with a description of the issue (perhaps some plots, or just a summary), and the correct calculation (I think mrsol could be used by data users to recalculate mrsos?).

Linear trends would be helpful and add to the errata if that's easy for you.

Thoughts?
Chloe

@blimlim
Copy link
Author

blimlim commented Oct 31, 2023 •
edited
Loading

Hi Chloe,

Apologies for the delay in getting back to you on this! I think that sounds like a good plan. Do you know if mrsol was provided with the rest of the cmip fields? In some of the runs I've looked in the past I haven't seen it come up.

I finally got around to looking at those trends. I thought it might make sense to look at each of the three SST forcings separately, each with 10 members.

The following show the ensemble mean annual trends for the corrected mrsos on the left, and the difference between the uncorrected and corrected versions on the right:

download
download-1
download-2

When the corrected trend is positive, the uncorrected trend seems to be not positive enough, and when the corrected trend is negative, the uncorrected one seems to be not negative enough. Checking the difference in magnitudes, the uncorrected trend nearly always has too small a magnitude:

It's looking like the uncorrected trends are generally scaled down compared to the corrected ones.

In terms of percentage errors, I tried calculating the differences in ensemble mean uncorrected and corrected trends, divided by the ensemble mean trend. The results look a little strange, so I'm hoping I haven't made any errors in my calculation:


There are scattered very high values (eg ~11000), which I think occur when the correct trend is very small (I'll try to check this a bit further), but apart from that, I'm surprised how flat these plots are. In a lot of places, it seems very much like the uncorrected trend ≈ 0.8 * the correct trend.

Let me know if you think there's any thing worth looking a bit more into, or making plots of!

Cheers,
Spencer

@chloemackallah
Copy link
Member

chloemackallah commented Nov 21, 2023 via email

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees
No one assigned
Labels
None yet
Projects
None yet
Milestone
No milestone
Development

No branches or pull requests
3 participants
@paolap @chloemackallah @blimlim