en/latest/examples/CloseHeatBudget_POP2

[utterances-bot]

Calculate POP2 heat budget using xgcm

https://pop-tools.readthedocs.io/en/latest/examples/CloseHeatBudget_POP2.html

[Zhangqs-hub]

Hello, I would like to ask, are DZU, DZT and dz the same? There is no DZU and DZT in my pop output.

[dcherian]

You'll have to compute DZU and DZT. I have a mostly complete PR in #44 if you want to test it out.

[Zhangqs-hub]

Thank you very much for your reply. I'm studying cesm recently and want to do heat budget analysis. I read your reply and tried to use pop_tools.grid.calc_dzu_dzt() to calculate 'DZT' and 'DZU', but 'DZBC' doesn't exist in my POP output file. The POP grid I used is gx1v6. May I ask if you have any files including the 'DZT' and 'DZU' corresponding to this grid? If yes, can you send me in the format of 'nc' file? Thank you.

[ALDepp]

Hi there, we don't have DZT and DZU readily available for each grid as far as I'm aware, and I haven't used that grid myself. I'm asking around for the files containing DZBC and will let you know when I find out more.

[ALDepp]

Hi again, you are looking at 1° model output, correct? In that case DZT[k,:,:] = dz[k], because the 1° model does not use partial bottom cells.

[Zhangqs-hub]

Thank you very much for your reply. It is very helpful to me. I found that gx1v6 is a 1° grid on average, which is also a three-pole grid. According to your reply, if DZT[k,:,:] = dz[k], the program can run normally, but it seems that 'SUM' and 'TEND_TEMP' are not completely equal (not closed). For your example code above, except for DZT and DZU, I did not make any other modifications. This may be incorrect. My calculation code can be found here : https://github.com/Zhangqs-hub/pop_tools-test. But the POP output file used is a bit big, more than 2G, I don’t know how to send it to you. If you have any good ideas, please let me know. Thanks again.

[dcherian]

@Zhangqs-hub if you put that code in a notebook, run it, then upload it; we'll at least be able to see the figures even if we cannot actually run it.

[Zhangqs-hub]

OK, thank you, let me try.

[Zhangqs-hub]

I have run and uploaded the notebook file, please check it out, thank you.
https://github.com/Zhangqs-hub/pop_tools-test

[ALDepp]

@Zhangqs-hub , I suggest plotting some maps of the fields around the depths where the budget doesn't close. It's not obvious to me from the line plot why it doesn't.

[Zhangqs-hub]

Thank you for your suggestions. I drew the maps at z_t=12500 and uploaded them to the previous github. I would like to ask, can the small error between 'SUM' and 'TEND_TEMP' be the reason for the grid conversion, such as converting the three-pole grid of POP to xgcm will produce a small deviation, which will cause it to not close? In addition, is there any problem with the calculation and creation of my DZT and DZU?

[klindsay28]

@Zhangqs-hub, are you looking at model output that you are generating? If so, you might consider writing out double precision output and seeing if the magnitude of your residuals goes down. If the residuals don't go down in magnitude, then you can probably conclude that a term is missing or mishandled.

Double precision output can be enabled by running the following command in CASEROOT

./xmlchange POP_TAVG_R8=TRUE
./case.build --clean ocn
./case.build

That said, it doesn't look like your computation includes contributions from HDIFB_TEMP and RF_TEND_TEMP. These terms are zero in high-res configurations, and are thus not needed in CloseHeatBudget_POP2.html, but they are non-zero in our 1 degree model.

[Zhangqs-hub]

Thanks for your reply. I followed your reply and carried out the output of double precision, and tested it, and found that it is indeed still not closed. I checked the output file and found that there are indeed two variables, RF_TEND_TEMP and HDIFB_TEMP. Among them, the global value of RF_TEND_TEMP is smaller, 1e-10 degC/s, and the magnitude of HDIFB_TEMP can reach 1e-6 degC/s.
Maybe it is not closed due to lack of HDIFB_TEMP. How to put HDIFB_TEMP into the POP2 heat budget analysis code? Is this part of VDIF? I'm so sorry I don't know much about these variables and their physical meaning.

[klindsay28]

HDIFB_TEMP is a tendency from lateral mixing through the bottom of the cell, so I suggest including it with budget["HDIF"]. It is non-zero because of mixing along isopycnal surfaces. The sign convention is positive down, and the downward flux is normalized into a tendency by dividing the flux by dz. So the equation for net tendency at layer k from this process is

(HDIFB_TEMP[k-1,:,:] * dz[k-1] - HDIFB_TEMP[k,:,:] * dz[k]) / dz[k]

I'm not familiar enough with xgcm to implement this equation with xgcm. That said, HDIFB_TEMP * dz has the same sign conventions and grid location as DIA_IMPVF_TEMP. So to get the net tendency from this process, apply what you do to DIA_IMPVF_TEMP to HDIFB_TEMP * dz, with a zero-flux condition at the surface.

RF_TEND_TEMP is a plain tendency and can be added as a separate term to your tendency reconstruction.

[Zhangqs-hub]

Thank you very much for your clear reply, but there seem to be some problems after I tried. I have uploaded the notebook file (https://github.com/Zhangqs-hub/pop_tools-test),and you may be able to check. Thanks again for your help.

[klindsay28]

Please try replacing

budget["HDIFB_TEMP"] = gridxgcm.diff(dsxgcm.HDIFB_TEMP * dsxgcm.VOL.values, axis="Z") / dsxgcm.VOL
budget["HDIFB_TEMP"][:, 0, :, :]=0.0

with

budget["HDIFB_TEMP"] = -gridxgcm.diff(dsxgcm.HDIFB_TEMP * dsxgcm.VOL.values, axis="Z") / dsxgcm.VOL
budget["HDIFB_TEMP"][:, 0, :, :] = -dsxgcm.HDIFB_TEMP.isel(z_w_bot=0)

I think this more closely mimics the treatment of DIA_IMPVF_TEMP.

[Zhangqs-hub]

According to your modification, the thermal balance analysis has been closed!This's so great!
But I don’t know the meaning of RF_TEND_TEMP. What is its physical meaning? Which term of the temperature equations corresponds to RF_TEND_TEMP?
Thank you very much!

[klindsay28]

That's great news that the budget is closing.

Regarding RF_TEND_TEMP...
POP uses a leap-frog time stepping scheme, which leads to numerical instabilities. We suppress these instabilities in POP with a Robert filter (RF is an abbreviation for Robert filter). RF_TEND_TEMP is the tendency from this filter. It is mentioned briefly in section 2.2 of the the CESM2 overview paper (Danabasoglu et al., JAMES, 2020, doi:10.1029/2019MS001916). There are additional references to descriptions of the filter there.

The contribution of RF_TEND_TEMP to the budget analysis decreases in magnitude as the time duration of your budget analysis increases.

Prior to CESM2, we used a time filtering approach in POP that didn't lead to a contribution in tracer budgets. We switched to the Robert filter in POP in CESM2 to enable hourly coupling of the ocean to the atmosphere. If we used hourly coupling with the previously used time filtering approach, the computational cost of POP would have gone up by a factor of 3, because of extra time steps required by the previous approach. (The previous approach requires at least 3 time steps per coupling interval.)

In the hi-res configuration analyzed in CloseHeatBudget_POP2.html, the ocean model takes time steps much smaller than 1 hour. So no extra time steps are needed to use hourly coupling and the previously time filtering approach. So we use the previous approach in the hi-res configuration and RF_TEND_TEMP doesn't appear in the budget.

[Zhangqs-hub]

Thank you very much for your detailed answers!
There is a small question. Before proceeding with the heat budget, I want to interpolate all the variables on the gx1v6 three-pole grid to the standard grid points, such as TLON=1, 2,..., 359, 360, TLAT=-90, -89,..., 89, 90, is there any function or method available in pop_tools or xgcm that allows me to conveniently rearrange the grid?

[jn533213]

Looking at CESM2, it doesn't seem to have a HDIFE_TEMP term like in CESM1. All I could find was HDIFT which was the vertically integrated horizontal mixing tendency. Did you have to calculate your HDIFE_TEMP and HDIFN_TEMP terms from other ones?

[ALDepp]

@Zhangqs-hub I do not recommend interpolating before closing the budget. Have a look here for interpolation methods: http://xarray.pydata.org/en/stable/user-guide/interpolation.html

@jn533213 It sounds like you might not have all the budget terms saved. I am not aware of all the CESM simulations that are available, but if you provide more info on which simulation you are looking at I am happy to ask around whether the budget terms are available.

[jn533213]

I was looking to recreate the heat budget that you have done using the data from CESM2 LENS2 linked here: https://www.cesm.ucar.edu/projects/community-projects/LENS2/data-sets.html only over a different time period and location. Is this the proper data set to go forward with? If not, would you be able to provide a link to the specific data sets used in the pop_tools.DATASETS.fetch function?

[Zhangqs-hub]

Thank you for your reply. I've had some new problems.
1. I found that the amount of heat budget terms I got might be a bit of a problem. I have uploaded the notebook file (https://github.com/Zhangqs-hub/pop_tools-test), which I converted from ℃/s to ℃/month, and adjusted the colorbar range based on the value of the data. I have found that the value of advection terms is very large, for example, WTT values are greater than 30 ℃/month near the equator, VNT also has a very small value near the equator, less than -30 ℃/month. However, the total advection seems normal, and in our validation, the heat budget terms are closed. I'm not sure if there was a mistake, I'd like to ask you for advice.
2. I'm not familiar with python, so I use a simple command, budget.to_netcdf('budget_test_ctl.nc'), output the heat budget terms as nc file, and then process them with matlab, interpolate them into the new grid, such as: [slon, slat]= meshgrid (30:389, -89:90); TEND_TEMP2 = griddata(TONG, TLAT, z_t, TEND_TEMP, slon, slat), I compare the interpolation before and after no significant change, do not know whether this will affect our analysis of the heat budget?
3. I would like to express the heat budget terms of POP2 output with tracer transport equation, so as to explain more clearly the meaning of each term. Does this equation that represents the heat budget terms in POP2 are right? (I'm sorry I don't know how to show the equation here, so I uploaded it to Github in pictures. The terms on the right-hand side, respectively, are horizontal advection, vertical advection, vertical diffusion, horizontal diffusion and penetrative shortwave radiation)

[ALDepp]

@jn533213 I am inquiring for you whether the full heat budget has been saved (in the monthly output there is a term HDIFT that I am not familiar with, maybe it is the sum of HDIFE and HDIFN).

@Zhangqs-hub
regarding 1: the large terms seem to cancel each other, this tells you that in this model there is a lot of temperature advection in the zonal, meridional and vertical, due to strong currents and large temperature gradients. I can't spot a mistake from having a brief look.
regarding 2: interpolating after calculating should be fine (and in fact it seems you have checked that this does not make a difference).
regarding 3: The equation you link is the heat budget, as it is calculated here. The sign for the vertical processes depends on the definition of the direction of the z coordinate (can be negative or positive, therefore there can be different signs in the terms containing d/dz). The vertical diffusion term includes a non-local term from the KPP parameterization that you have not explicitly included in your equation, you can compare the formulation here, Equation 6, and attribution in 9 (the arrows only work in the pdf).

[Zhangqs-hub]

Thank you for your reply and it helps me a lot.
Regarding question 1, I found that near 40W, the advection in the x direction is discontinuous, resulting in discontinuity in the total advection here (see fig1, 40W corresponds to the vicinity of horizontal coordinates 20). Checking the code for calculating the heat budget, I found that it may be due to the boundary in the POP2 gx1v6 grid (I mean, around 40W is the horizontal boundary of POP TLONG, that is, the first row of TLONG is 320.5 and the last row is 319.4). Curiously, for HDIF, no such discontinuity was found near 40W(see fig2). I have not been able to solve this problem at present, I hope to get your help.
Regarding question 3, I’m sorry, I can’t see the webpage you shared. This seems to require a login username and password. I checked POP reference manual and found that the flux operator is used in the transport equation(see fig4).
But what difference does it make between the heat budget equations I've seen in many papers, as shown in fig 3? Which one do I use to illustrate the physical significance of each item better?

[Zhangqs-hub]

Sorry, the picture can be found here, https://github.com/Zhangqs-hub/pop_tools-test

[ALDepp]

@jn533213 You actually told me what HDIFT was and I had forgotten, vertically integrated. I'm afraid there are no other heat budget terms saved than what seems to be available to you. So you won't be able to close the budget, sorry about that.

[ALDepp]

@Zhangqs-hub great that you investigated! Since you are not taking a gradient in x for HDIF it doesn't surprise me that you don't see the same behaviour there, if it is due to the boundary. For the zonal advection term you are ((gridxgcm.diff(dsxgcm.UET * dsxgcm.VOL.values, axis="X")), so if you hit a boundary there might be a problem. Try setting periodic=True up where you convert your dataset to an xgcm compatible dataset (and probably change boundary={"X": "extend", "Y": "extend", "Z": "extend"}, to boundary={"Z": "extend"},). Let me know how it goes!

[Zhangqs-hub]

Thank you for your reply!

1. According to your suggestion, I only set periodic=True and the problem has been solved. fig5 is the new result, it seems to be fine.
2. For HDIFE_TEMP, it seems that the calculation method is the same as that of UET, and the gradient of the x direction is calculated( budget["HDIFE_TEMP"] = gridxgcm.diff(dsxgcm.HDIFE_TEMP * dsxgcm.VOL.values, axis="X") / dsxgcm.VOL), but look at the previous results(fig2), it seems that HDIF at 45W doesn't show the same discontinuity as the total advection(fig1). However, in fig2, there is a band distribution along the latitude, and I'm not sure if this is the same mistake. So I divided HDIF into three terms, as shown in fig6, and in the HDIFN_TEMP(fig6 lower left), it was obvious, and I'm not sure if it's normal or why.

[andersy005]

(I'm sorry I don't know how to show the equation here, so I uploaded it to Github in pictures.

@Zhangqs-hub, using the following syntax, you can embed the images in the markdown text:

![](image_url)

for e.g.

![](https://github.com/Zhangqs-hub/pop_tools-test/blob/main/fig4.png?raw=true)

[Ryohei17]

@jn533213 @ALDepp Unfortunately, CESM2 LENS did not save all heat budget terms as monthly output. But, you can close monthly heat budget precisely by using alternative terms to HDIFE_TEMP, HDIFN_TEMP, and HDIFB_TEMP.
CESM2 LENS monthly output has,
(Tendency) TEND_TEMP
(Advection) UET, VNT, WTT
(Horizontal diffusion) ISOP_ADV_TEND_TEMP, Redi_TEND_TEMP, RF_TEND_TEMP, SUBM_ADV_TEND_TEMP
(Vertical diffusion) DIA_IMPVF_TEMP
(Air-Sea flux and SW penetration) SHF, SHF_QSW, QSW_3D
(Ocean internal heating) QFLUX (this term is missed in the pop_tool example page)
Sum of the four horizontal diffusion terms is equivalent to 'HDIF' estimated from HDIFE_TEMP, HDIFN_TEMP, and HDIFB_TEMP. And then missing term is only 'KPP_SRC_TEMP' in vertical diffusion terms, so, the residual of heat budget computed from above terms should be contribution from 'KPP_SRC_TEMP'.

[tlsw231]

Shouldn't the derivatives be computed using gridxgcm.derivative rather gridxgcm.diff? I see the latter being used to obtain the advection terms but the 'diff' operation just forms the numerator in the derivative calculation, right?

[tlsw231]

I am trying to plot the globally averaged terms in the heat budget. I am not sure I fully understand how the weighting is being done here. Here is how I think of the sequence of operations:

1. Multiply each budget term by the local cell volume.(degC * cm^3/s)
2. Sum the term in 1 over (nlon,nlat) for each vertical level. (degC * cm^3/s)
3. Add together the cell volumes for each vertical level. (cm^3)
4. Divide answer in 2 by answer by 3 to get the globally averaged vertical profile in degC/s units.

But in the code above, the KPP_SRC_TEMP term does not appear to be weighted. I cannot just sum this term as is over (nlon,nlat) without weighting first, am I correct?

I am also confused why each term is divided by the local cell volume (dsxgcm.VOL). I thought the weighted mean requires summing the numerator (weighted terms) and the denominator (weights) separately before dividing one by the other.

[tlsw231]

@ALDepp, I forgot to tag you in the previous two posts of mine.

[ALDepp]

Hi @tlsw231, the reason we use grid.diff and multiply and divide by the volumes ourselves is the way POP outputs the fluxes. It performs a division by the cell area before saving the terms, which you would not account for if you used grid.derivative. Note that we also multiply by dsxgcm.VOL.values and then divide by dsxgcm.VOL. This is due to the same issue, there is a mis-alignment in the grid if you wish in this output term that xgcm would not like, and we are getting around it this way. Aa far as I know this is a POP specific thing, likely for other models you can use grid.derivative.
You can sum the KPP_SRC_TEMP term because it is output as "degC/s", while the DIA_IMPFV_TEMP term is output as "degC cm/s", so area weighting is not necessary.

[dcherian]

@tlsw231 it'd be nice to update the notebook with the information in @ALDepp's comment. Would you mind sending in a pull request please?

[dcherian]

I agree that it would be much nicer to use grid.derivative instead of grid.diff.

One way to do this would be to rescale the output i.e. undo the division by cell area, and then use xgcm's metric-aware functionality. I bet that would be less confusing if we can get it to work. @tlsw231 are you for experimenting with this approach?

[tlsw231]

@cherian_deepak I am glad to contribute. I am somewhat new to github, just created an account only very recently. To create a PR, I just follow the instructions in the 'Preparing pull requests' here, right?

I can experiment with grid.diff vs grid.derivative, this might take a little while though. I have a couple of meetings coming up in the next 10-12 days, hence the flurry of questions:). One question I have is what dataset to use for validation? Based on the comments section for this notebook, it appears terms that vanish for 0.1 deg POP output can be quite significant for 1 deg POP output. In my case, the 0.1 deg output does not have the TEND_TEMP` term, so I have been trying to use the restart files at the beginning and end of the averaging interval to compare the LHS and RHS in the budget. Haven't been successful so far but that might be because I have not been doing the weighting correctly.

For creating a PR with just the extra documentation added in, I was wondering what is the standard protocol regarding initiation of the PR. Is it more appropriate for the original author of the notebook to do this? I feel like that person has done all the hard work, even supplying the necessary additional documentation. I would just be doing a copy-paste job :)

[tlsw231]

@dcherian I am glad to contribute. I am somewhat new to github, just created an account only very recently. To create a PR, I just follow the instructions in the "Preparing pull requests" here, right?

I can experiment with grid.diff vs grid.derivative, this might take a little while though. I have a couple of meetings coming up in the next 10-12 days, hence the flurry of questions:). One question I have is what dataset to use for validation? Based on the comments section for this notebook, it appears terms that vanish for 0.1 deg POP output can be quite significant for 1 deg POP output. In my case, the 0.1 deg output does not have the TEND_TEMP` term, so I have been trying to use the restart files at the beginning and end of the averaging interval to compare the LHS and RHS in the budget. Haven't been successful so far but that might be because I have not been doing the weighting correctly.

For creating a PR with just the extra documentation added in, I was wondering what is the standard protocol regarding initiation of the PR. Is it more appropriate for the original author of the notebook to do this? I feel like that person has done all the hard work, even supplying the necessary additional documentation. I would just be doing a copy-paste job :)

[tlsw231]

@ALDepp Thanks for the detailed explanation, really appreciate it. Once you assemble the terms in the budget array, they all have the units degC/s. So, I should be able to algebraically sum all these terms over (nlat,nlon) without any weighting to get the horizontally averaged vertical profile. Am I interpreting your comments correctly?

[tlsw231]

@ALDepp I tried working with 1-degree POP output from CESM version 1.3 but am still having a hard time balancing the budget. I have uploaded a notebook here. One difference between my notebook and yours is that I am first computing the time mean of all the terms but don't see why that should create problems. In the right panel of the plot in the notebook, there is a large discrepancy between the tendency obtained from adding the terms on the right side of the equation and that from the variable TEND_TEMP saved as part of the POP output. I have tried adding/removing: (i) GM terms (ii) submeso terms; (iii) HDIFB (but retaining HDIFE and HDIFN) but always get an order of magnitude difference in the right panel.

One of the main sources of my confusion is what all does HDIFB include. In your notebook you omit HDIFB but get an excellent balance. Based on an earlier comment by @Ryohei17, HDIFB must also be included. That is the case in the notebook uploaded by @Zhangqs-hub, who also gets a good balance. If I read through the comments within the NCL code generate_TEMP_budget.map.monthly.ncl provided by NCAR, they say HDIFB includes both resolved and subgrid. Does that mean I should not add TOT_ADV to the sum? (tried that too but to no avail).

[tlsw231]

@klindsay28 I have been trying to balance the POP heat budget with 1-deg cesm output (version 1.3, see previous comment). I read your replies earlier in this thread. Would you have any idea what am I doing wrong (link to notebook)? Thanks!

[jn533213]

Is there a way to calculate the individual components of the air-sea heat flux minus the shortwave radiation rather than just the net air-sea heat flux minus the shortwave radiation? I'm not sure how to handle the DIA_IMPVF_TEMP when using individual components ie. sensible, latent heat flux, longwave radiation.

[HXJ1008611]

@klindsay28 From your reply , I learn that in 1 degree model we need HDIFB_TEMP; in CEMS1 the time filtering approach used in POP didn't lead to a contribution in tracer budgets, so in CESM1, there's no RF_TEND_TEMP. I try to recreate the heat budget that you have done using the data from CESM1 LENS linked here: https://www.earthsystemgrid.org/dataset/ucar.cgd.ccsm4.cesmLE.ocn.proc.monthly_ave.html
After considering all the terms, the budget is still not closing. is there any other term I need to add?
In the website, they don't provide the TEND_TEMP term. I would like to know how to calculate the TEND_TEMP by myself.
You can find my work here: https://github.com/HXJ1008611/pop-tool-heat-budget
Thanks for your help.

[Yefee]

@HXJ1008611 Did you finally close the heat budget in the CESM1-LENS? Have the same problem here.

[HXJ1008611]

我已收到您的邮件。我看到后会及时回复，谢谢。

[Zhangqs-hub]

@ALDepp Hello, may I ask if the budget['UET'], budget['VNT'] and budget['WTT'] calculated by xgcm correspond to -udT/dy, -vdT/dy, and -wdT/dz, or -d(uT)/dy, -d(vT)/dy, and -d(wT)/dz?

[HXJ1008611]

我已收到您的邮件。我看到后会及时回复，谢谢。

[ALDepp]

@Zhangqs-hub UET, VNT, and WTT are outputs from the model, they are online accumulated over the 5day period that I have output for. I'm assuming you don't have those outputs and want to calculate the equivalents? Note that then you are calculating from means and will have a missing contribution from eddy effects, which are accounted for when you use the total transports output by the model.
UET is eastward transport of heat, and VNT northward transport. You have dy everywhere, so there needs to be dx in your u component. If I'm understanding your question right you wonder whether you can take the difference between the transport out of the east facing face of the box minus the transport into the west facing box (for the example of UET) or whether you use the velocity times the temperature gradient, is that what you are asking? Are you interested in doing your calculations per grid point or over a box?