Viscous Dissipation in the NASA/GISS ModelE Code is Wrong …
it is not even viscous dissipation. There is also an energy imbalance in the NASA/GISS ModelE code due to the error in the ‘viscous dissipation’. The energy imbalance is about the same magnitude as the imbalance associated with increased concentrations of CO2 in the atmosphere.
Updated January 14, 2009, down near the end.
This is another continuation post about viscous dissipation. A followup on the discussions over at Luica’s and additional information related to both the general issue and some specific aspects of the NASA / GISS ModelE model / code.
I originally ran across a ModelE code segment that purported to be the viscous dissipation as mentioned in this post. There is a some discussion with NASA / GISS employee Gavin Schmidt there in the comments. That led to a more general discussion in this post, and given the estimated magnitude of the viscous dissipation, I posted this thread about the general issues of conservation of mass and energy in model equations.
At the time that this was ongoing, my preliminary conclusions were that the quantity purported to be viscous dissipation in the ModelE code was in fact not viscous dissipation. I also decided that the ModelE model / code had the potential to be missing a not-insignificant contribution to the thermal-energy form of the energy balance equation.
Based on additional reading, I will now say that both of my preliminary conclusions are correct. The ModelE model / code does not include accounting of the viscous dissipation, and the magnitude of the missing contribution is about the same value as the energy imbalance introduced by increased concentration of CO2 in the atmosphere. The papers by Becker and colleagues that I cited in this post, plus some others cited in these, provide important information. The papers that I find useful are Becker 2001, Becker 2003, and Burkhardt and Becker 2006. There is an earlier article by Fielder; Fiedler B. H., 2000: Dissipative Heating in Climate Models. Quart. J. Roy. Meteor. Soc., 126, 925–939, but I haven’t paid the bucks to get it.
Becker cites the following to conclude that the physical processes in the atmosphere correspond to viscous dissipation of about 2 W/m^2:
It is noteworthy that a simulated mean dissipation rate of about 1.9 W m–2 fits well within the range of existing estimates based on observational analyses of the general circulation. Oort (1964) finds a value of 2.3 W m–2 (see also Lorenz 1967, chapter V). Some modern textbooks suggest somewhat lower dissipation rates between 1.8 and 1.9 W m–2 (e.g., James 1994, chapter 5.3; Prandtl et al. 1989, chapter 8). Thus, with regard to the total kinetic energy typically contained in the atmosphere (~7 × 1020 J), neglecting dissipation in conventional GCMs means that within about 8 days an equivalent amount of kinetic energy is removed from the flow without being dissipated into heat (see also Pichler 1986, chapter 9).
Becker and colleagues base their developments on two fundamental aspects of fluid motions as follows; (1) conservation of angular momentum requires the stress tensor to be symmetrical, and (2) the viscous dissipation must always be a non-negative positive-definite contribution to the thermal-energy forms of the energy conservation equation.
I think Becker 2003 discusses explicitly the Total-Potential-Energy Change-Kinetic-Energy approach mentioned in passing by Gavin Schmidt:
It is worth mentioning that an attempt to include frictional heating in GCMs can also be found in Hamilton (1996) or Kiehl et al. (1996, sections 3b and 4d). With this alternative method, no knowledge of the stress tensor is required. One rather assumes that, locally, the rate of change of kinetic energy due to friction enters the thermodynamic equation of motion with negative sign. This assumption may be interpreted as the local counterpart of the fact that, in the global mean, the frictional loss of kinetic energy balances the frictional heating (e.g., Pichler 1986, chapter 9; Smagorinsky 1993; see also section 2b). However, assuming local equivalence implies that there is no turbulent stress at all acting at the resolved scales. Also, the local frictional heating rate has arbitrary signs.
I’ll have to track down the references he cites to see if this is what Schmidt is saying.
Based on the information in these papers and the coding in NASA / GISS ModelE code I now conclude
1. The quantity called ‘dissipation’ in the ModelE model / code is in fact not viscous dissipation. As I noted in the original post, I suspect the approach in the code implicitly contains numerical artifacts in addition to any ‘physical’ accounting.
2. The above error in the ‘dissipation’ leads to an energy imbalance of about 2 W/m^2 in the NASA / GISS ModelE calculations. This imbalance is about the same as that associated with increased concentrations of CO2 in the atmosphere.
3. The potential errors in the calculation of viscous dissipation seem to have been recognized only recently now early in the 21st Century. This situation leads me to conclude that the results of many GCM calculations used by the IPCC, especially in the first ‘Annual Reports’, very likely contain the same error as ModelE ( speculation on my part ).
4. The Total-Potential-Energy Change-in-Kinetic-Energy approach can indeed produce negative numbers for the ‘viscous dissipation’; a most unusual outcome, in my opinion.
In regard to 3. above, Boville and Bretherton 2003 state:
The spectral Eulerian core in CAM2 includes a biharmonic horizontal diffusion operator that cannot be represented by a symmetric stress tensor and therefore the kinetic energy dissipation cannot be correctly defined, as noted by Becker (2001).
Where Becker (2001) is the paper cited above. The authors note that the Change-in-Kinetic-Energy approach has the potential to give cooling of the fluid. These authors also give the viscous dissipation to be about 2 W/m^2.
Reverse engineering model equations from coding is a process that has very high probabilities for errors; especially when even the most rudimentary documentation is not available. I think there is more light than heat here, but under the present circumstances I can not be certain. Additionally, maybe I’m biased because these papers are in agreement with preliminary conclusion that I obtained 🙂
Update:
Here is the citation and abstract of the Fiedler (2000) paper:
B. H. Fiedler, “Dissipative Heating in Climate Models”, The Quarterly Journal of the Royal Meteorological Society, Volume 126 Issue 564, Pages 925 – 939, 2000.
Abstract
The various strategies towards dissipative heating in meteorological models are reviewed. The strategies are implicit formulations, explicit formulations, and exclusion. A thermodynamic formulation based on dry static energy, which has long been used in hydrostatic models, implicitly includes dissipative heating. Many modem non-hydrostatic mesoscale models, being based on potential temperature, do not include any dissipative heating. The National Center for Atmospheric Research community climate model (CCM3) explicitly formulates the dissipative heating for the resolved motion, which is based upon potential temperature, but, uses the implicit formulation for the sub-grid motion.
Those who want to use the mesoscale-model formulation as a cloud-resolving climate model are advised to consider the effect of dissipative heating. We show a dry radiative-convective model, formulated without dissipative heating, in which climate equilibrium occurs with a net inward radiative flux of about 18 W m-2 at the top of the atmosphere, and a surface temperature that is 8 K too low.
Models Methods Software 
I just searched Chapter 8 of AR4, Climate Models and Their Evaluation, for Becker. Got no hits.
This chapter should be cited as:
Randall, D.A., R.A. Wood, S. Bony, R. Colman, T. Fichefet, J. Fyfe, V. Kattsov, A. Pitman, J. Shukla, J. Srinivasan, R.J. Stouffer, A. Sumi and K.E. Taylor, 2007: Cilmate Models and Their Evaluation. In: Climate Change 2007: The Physical Science Basis. Contribution of Working Group I to the Fourth Assessment Report of the Intergovernmental Panel on Climate Change [Solomon, S., D. Qin, M. Manning, Z. Chen, M. Marquis, K.B. Averyt, M.Tignor and H.L. Miller (eds.)]. Cambridge University Press, Cambridge, United Kingdom and New York, NY, USA.
The NCAR paper Nick linked similarly suggests doing this right is “difficult”. It’s a fairly new paper, and seems to be discussing something people only recently got a handle on.
I still think we can obtain quasi-steady solution. However, there will be errors if this is not dealt with properly. How much? Beats me. It’s good you brought up the topic.
My background is in hard and medium-soft rocks, so please ignore this question if it has no significance to the discussion.
We had viscosity problems modelling wet sediment movement. The boss used to place an empty, wet beer glass on the bar and heat it with a lighter. The escaping gas facilitated “floating” the glass and its flow from one part of the bar to another took far less energy than when it was dry.
Is there an equivalent concept in atmospheric turbulence where thin shears facilitate the bulk movement of large gas volumes? If not this mechanism, is it risky to assume that viscous dissipation can be modelled on smoothly-varying properties without allowance for formation of bundles (that might look like cells in clouds, for example)? Do zones of slippery turbulent flow bound parcels of sluggish non-turbulent atmosphere?
Geoff :
No , no and no .
Liquids and gases don’t behave in the same way even if they obey both to Navier Stokes equations .
The turbulence is much higher in gases .
Besides your example is not even fluids – it’s solid/solid friction .
In Kolmogorov theory you can consider that energy dissipation (that happens at very small scales)happens through interaction of small eddies with isotropic and homogenous distribution .
There is a hierarchy of scales .
At large scales inertial effects dominate (that’s where the energy of the fluid is) .
At small scales viscous efects dominate (that’s where the dissipation of fluids is) .
And the whole dynamics is a process where energy casacades from large scales that are the entrance door of the system to small scales that are the exit door of the system .
That’s why the spectral methods are almost always used when people speak about dissipation (obviously it is important to know the energy dependence on the wave number : E(k))
Hi Dan,
I decided to download the Model E source code and check this business out for myself. You are right – ediff has nothing to do with viscous dissipation. It appears to be correction to the temperature based on changes in kinetic energy (as stored in the array DKE).
Consider this code snippet from ATMDYN.f
C**** Add in dissipiated KE as heat locally
!$OMP PARALLEL DO PRIVATE(I,J,L,ediff,K)
DO L=1,LM
DO J=J_0,J_1
DO I=1,IMAXJ(J)
ediff=0.
DO K=1,KMAXJ(J) ! loop over surrounding vel points
ediff=ediff+DKE(IDIJ(K,I,J),IDJJ(K,J),L)*RAPJ(K,J)
END DO
T(I,J,L)=T(I,J,L)-ediff/(SHA*PK(L,I,J))
END DO
END DO
END DO
!$OMP END PARALLEL DO
END IF
You will notice that DKE is defined above this code fragment as
DKE(I,J,L)=0.5*(U(I,J,L)*U(I,J,L)+V(I,J,L)*V(I,J,L)
* -USAVE(I,J,L)*USAVE(I,J,L)-VSAVE(I,J,L)*VSAVE(I,J,L))
As far as I can tell, USAVE, VSAVE are the flow velocities prior to the Shapiro Filtering that is done to them (this filtering is probably worthy of its own post!). So, the delta KE here is merely the difference between the KE of the filtered velocities with the KE of the unfiltered velocities.
Why do this? The filtering must alter the velocity field enough to warrant the energy modification – that is, it is purely a numerical artifact and NOT something based on physics (again, as far as I can tell). I suppose that if you don’t correct the energy, these small differences (actually numerical errors) will, over the course of the time-marching process, become large and corrupt the solution. Remember that climate models are run over *decades* worth of simulation time using relatively small time steps to maintain stability.
By the way, if you look at the temperature modification that is done (T(I,J,L) = …), the parameter SHA is the specific heat at constant pressure of dry air (a constant in the model), and the array PK(…) is a pressure function which appears to be
PK = pressure^((gamma-1)/gamma)
where gamma is the specific heat ratio. I was trying to figure out where they got this relation from – it appears to be related an isentropic relations for compressible fluids
(T1/T2) = (P1/P1)^((gamma-1)/gamma)
but how it is used here is a little weird.
If you go further down in ATMDYM.f, you’ll come across the subroutine DISSIP. Here the same ediff calculations and temperature corrections are repeated. This is only called from the main program, as contained in the file MODELE.f, and it appears that this time DKE(…) comes from a surface-based calculation:
C**** ADD DISSIPATED KE FROM SURFACE CALCULATION BACK AS LOCAL HEAT
CALL DISSIP
I think DKE for this instance comes from the following code in ATMDYM.f
C**** WL is restricted to Wmax by adjusting X, if necessary;
C**** the following is equivalent to first reducing (U,V), if necessary,
C**** then finding the drag and applying it to the reduced winds
CDN=CSDRAGl(l)*xjud
IF (cd_lin) CDN=(X_SDRAG(1)+X_SDRAG(2)*min(WL,wmaxj))*xjud
X=DT1*RHO*CDN*min(WL,wmaxj)*GRAV*BYDSIG(L)*BYPIJU
if (wl.gt.wmaxj) X = 1. – (1.-X)*wmaxj/wl
C**** adjust diags for possible difference between DT1 and DTSRC
AJL(J,L,JL_DUDTSDRG) = AJL(J,L,JL_DUDTSDRG)-U(I,J,L)*X
DP=0.5*((PDSIG(L,IP1,J-1)+PDSIG(L,I,J-1))*DXYN(J-1)
* +(PDSIG(L,IP1,J )+PDSIG(L,I,J ))*DXYS(J ))
ang_mom(i,j) = ang_mom(i,j)+U(I,J,L)*X*DP
DUT(I,J,L)=-X*U(I,J,L)*DP
DVT(I,J,L)=-X*V(I,J,L)*DP
DKE(I,J,L)=0.5*(U(I,J,L)*U(I,J,L)+V(I,J,L)*V(I,J,L))*(X*X-2.*X)
Again, this has nothing to do with actual viscous dissipation, but appears to be a correction designed to prevent numerical errors (in the form of energy errors) from accumulating over long periods of simulation time.
What makes this all somewhat tricky is that I don’t know if any non-dimensionalization was employed in the equations, so trying to relate the code back to fundamental equations of fluid dynamics and thermodynamics requires figuring out the basic equations that are being solved and the forms employed for each of the dependent variables.
Hope this helps (a little).
cheers,
Frank
Frank – No, I think the DKE used by DISSIP() (which is indeed where the viscous dissipation is converted to heat) comes from the turbulence model in ATURB.f, line 470. This makes sense – they are calculating there the vertical momentum diffusion, which is effectively the viscous dissipation, and according to the comment, they compute the associated KE loss and store it in DKE for future heat adjustment, which I think is the DISSIP call.
Hi Nick,
Thanks for correction. Here is the code:
C**** Save additional changes in KE for addition as heat later
!$OMP PARALLEL DO PRIVATE (L,I,J)
DO L=1,LM
DO J=2,JM
DO I=1,IM
DKE(I,J,L)=DKE(I,J,L)+0.5*(U_3d(I,J,L)*U_3d(I,J,L)
* +V_3d(I,J,L)*V_3d(I,J,L)-U_3d_old(L,I,J)*U_3d_old(L,I,J)
* -V_3d_old(L,I,J)*V_3d_old(L,I,J))
END DO
END DO
END DO
!$OMP END PARALLEL DO
Again, this is NOT viscous dissipation in the sense of the dissipation function defined in the energy equation, which is:
Phi = tau(i,j) dui/dxj
where tau(i,j) is the viscous stress tensor, and dui/dxj represents derivatives of velocity components wrt the Cartesian coordinate directions. If you were to write this out in Cartesian component form, it would look like this:
Phi = mu*[2*(du/dx)^2 + 2*(dv/dy)^2 + 2*(dw/dz)^2 + (dv/dx + du/dy)^2 + (dw/dy + dv/dz)^2 + (du/dz + dw/dx)^2 – 2/3*(du/dx + dv/dy + dw/dz)^2]
where mu is the dynamic viscosity, and u,v,w are the Cartesian velocity components.
I don’t see that here. I do see that they are taking a difference between the kinetic energy of the flow velocities altered by the surface friction and the unaltered velocities at the near-surface cells. They hence must assume that the difference in kinetic energy is due to viscous effects, but I’d like to see a derivation showing a connection between delta KE and the viscous dissipation function above. Is there a paper which shows this?
Thanks,
Frank
Frank and Nick, Thanks for diving into the coding. I can do only so much of that when there are no documented equations to use as reference. I had decided that DKE was some kind of difference between velocity-components-squared, but had not found the details. I had also run into units problems with that PK variable. Can either of you get the units to work out?
I have yet to find any papers that develop this concept. It would seem to me that the Available Potential Energy, in contrast to the Total Potential Energy, should enter somehow. It is also noteworthy that the no-slip condition applied at a stationary surface would seem to imply that no work can be done by that surface.
Thanks again for your assistance.
Frank,
As you may know, there’s a parallel discussion of this at Lucia’s blog, and in my comment 8659 I responded to Dan on this topic. I think that the misapprehension that model E omits viscous heating is based on laminar flow thinking, which you’ve made explicit, whereas these are very high Re flows. You will not find a viscous stress term based on velocity gradient in the dynamic equations there, and rightly so. If you tried to calculate viscous stress, or viscous heating according to such formulae, the answer would be extremely small. Viscosity represents the diffusion of momentum, and if you rely on molecular diffusion, it is negligible. Instead momentum mixing occurs entirely through turbulent transport, and you need to compute an eddy diffusivity. This they do in ATURB.f. When you have done that, you can compute the deceleration, but also, directly, the KE reduction, which they do.
Just to reinforce this point, suppose you did try to calculate your Phi. You might have a grid cell 100 km square and a velocity difference of 10 m/s -> gradient 10^-4. Viscosity about 10^-5 (SI units). Result, about 10^-13 W/m3, or maybe 10^-10 W/m2 integrated over height.
Nick, the turbulent phenomena is given by the same form of the equations. The presence of turbulence is handled by a variety of methods; too many to review here.
I’ll put it on my to-do list to track down how turbulence is modeled in GCMs. I strongly suspect that a mixing-length approach is used along with parameterizations for a ‘turbulent viscosity’ model. Zeroth-order approach is to use a large ‘viscosity’ and carry out the calculations with equations that have the exact form as the laminar-flow equations.
The size of the spatial increments is, in my opinion, an issue that impacts all the discrete algebraic approximations; especially as convergence is still not checked or attempted in GCM calculations.
Hi Nick,
Thanks for the note. You can indeed add a viscous dissipation term which includes turbulence effects – this is done all the time in CFD! If you go through the derivation of the Reynolds Averaged NS equations, the equation for energy will contain the viscous dissipation term – in conservative form it is:
Div(tau * V + Reynolds-stress-related-term)
where tau is the stress tensor and V is the velocity vector and * indicates dot product and Div is the divergence operator. The Reynolds stress related term comes from the averaging and as usual involves averages of fluctuating terms. The standard approximation is to relate fluctuating terms to the mean flow velocities, velocity gradients, and a turbulent viscosity. And the turbulent viscosity can be many orders of magnitude greater than the laminar viscosity (even up to 10^4 – 10^5).
Now, in your estimate, I think you were considering the horizontal grid dimensions. Actually, the gradients I was thinking of here are ** vertical ** gradients in the planetary boundary layer. I doubt that the height of the first cell above the earth’s surface is 100 km!! I would guess that the cell height in the near-surface layer would be on the order of 100 – 1000 m, and gradually increase as you go up though the atmosphere. So, the estimate may be something like
dissipation –> mu-turb * (vel / dy)^2 = 10^3 * (10 / 100)^2 = 10
By the way, the units here are W / m^3 (N-s/m^2 * m^2/s^2 / m^2)
Now, having said all this, there’s nothing wrong with using a drag ocrrelation to apply a friction force and add that to an otherwise inviscid calculation (CFDers were doing that back in the 70s and 80s).
My only puzzle at this stage is how they relate the delta KE to temperature. My feeling is that the T(…) array is not physical temperature in Kelvin but some derived variable (probably a non-dimensional temperature). Same for the pressure. That would be the only way you could employ a thermodynamic relation involving P^((gamma-1)/gamma) – P must be a normalized pressure. Some real documentation would help here … sigh … ;^)
cheers,
Frank
Dan and Frank,
OK, yes, formally the process is similar, as long as we’re clear that its not actually viscous dissipation in the laminar sense that is being sought. Then, as I said, you compute an eddy diffusivity and proceed. GISS do this in ATURB.f, based on mixing length lscale(). The diffusivity for momentum is km() and for KE is ke().
They are reasonably only concerned about vertical diffusion, and so can solve just a 1D diffusion equation. For horizontal velocities U() and V(), this is done starting at line 431. They call a tridiag ode solver de_solver_main() to solve separately in U and V, with a bottom boundary condition based on friction, and slip at the top.
Having figured the new velocities after this diffusion solve, at line 470 they store DKE, the resulting KE increment. This is the DKE that is then used in the call to DISSIP() to compute the heat loss.
It’s all there.
Hi Nick and Dan,
Nick – Thanks for the information on the calculation of U and V.
Dan and Nick – I still don’t see the relationship between delta KE and viscous dissipation (frictional heating) in the boundary layer. Let me know if you find a derivation for this. Their formulation looks very ad hoc to me.
In conclusion, I don’t think we can call what they do “viscous dissipation”. After all, why would they adjust temperature due to changes in the velocity resulting from the Shapiro filtering in ATMDYN.f, as I outlined above? There is no physical basis for that – it’s simply a fix for numerical errors.
Frank
Frank, the Shapiro smoothing issue is quite different. Whatever the justification for the smoothing, it does make a separate change to the KE. They have chosen to use a separate bit of code (which you quoted) to make a balancing change in heat energy. That doesn’t affect the process whereby they compute a heat gain corresponding to the momentum diffusion KE loss. That’s real and clear physics, and is what people here were saying should be done. They’ve done it.
The boundary layer fluxes, including uflux and vflux, are computed in SURFACE.f. If you want to have some fun with that, you can note the comment on authorship. Anyway, these then become the boundary condition for the vertical diffusion calculation. Momentum diffuses down to the ground, and is lost to the air. KE is lost throughout the profile. The process is frictional, and so heat is created, and is added in where the KE was lost.
Hi Nick,
Thanks. Again it would be good to derive the equation which equates change in kinetic energy to a change in temperature. Anyone up for that? I looked that the kinetic energy equation (derivable by multiplying the momentum equation by the velocity vector and manipulating terms), and there may be a connection there.
So we agree that the KE change and subsequent temperature change due to Shapiro filtering is unphysical? It really is nothing more than a numerical fix. It would be interesting to visit this whole filtering business as it strikes me as being very adhoc way of fixing a bad time marching algorithm…
Frank
Frank,
I don’t have a view on the use of Shapiro filtering here, except a general predisposition to believe that smoothing is more often harmless than beneficial. But it may be addressing a stability issue. Some pro’s and con’s are discussed here</a”.
However this is rather a change of subject. This post has a standing assertion that Model E is wrong in not accounting for viscous dissipation, and I think it would be useful if we could settle that.
Hi Nick,
Thanks for the link to the paper which discusses the Shapiro filter. My opinion is that the smoothing is harmless, except where it’s not :^) Seriously, instabilities in numerically computed velocity or pressure fields are a sign of problems with the algorithm, and this kind of filtering is simply an adhoc fix with no physical basis (and as the authors mention, it also increases numerical diffusion of the solution). But we can save that for another day…
As for the current issue, I think we can establish that Model E does not have a proper viscous dissipation term in its energy equation. And whether of not the current method of equating differences in kinetic energy with changes in temperature is equivalent to viscous dissipation remains an open question. I for one would like to see a derivation of the mathematical relations upon which their method is based.
Frank