Energy leak?

[AndrewILWilliams]

Hi there.

Too long, didn't read: I appear to be getting steady global-mean temperatures, even when the TOA fluxes are extremely out-of-balance.

Model set-up: aquaplanet, perpetual equinox, two-stream gray rad (frierson), prescribed q-flux, no solar absorption in atmosphere, fixed ocean albedo=0.22

As far as I can tell, in this configuration the net TOA radiation is: ds['swdn_toa']*(1-0.22)-ds['olr'], because the albedo=0.22 (I've confirmed this by outputting the albedo variable).

My issue is that I end up with situations where the global mean net TOA radiation is $\sim\mathcal{O}(10)$ W/m2 out of balance, but the surface energy budget is satisfied to within less than 1W/m2, and the global-mean SST is steady to within 0.01K variations.

I suspect that I'm misunderstanding something about how to calculate the net TOA radiation? I can't see how, but everything else is quite simple so I don't think I'm messing up any other calculations...

Appreciate the help!

Andrew

=======================

Namelist params below:

exp.namelist = namelist = Namelist({
    'main_nml': {
        'days'   : 360,
        'hours'  : 0,
        'minutes': 0,
        'seconds': 0,
        'dt_atmos':200,# cut timestep in half when T0290_T2-60 runs crash
        'current_date' : [1,1,1,0,0,0],
        'calendar' : 'thirty_day'
    },

    'idealized_moist_phys_nml': {
        'two_stream_gray': True,
        'do_rrtm_radiation': False,   
        'convection_scheme': 'SIMPLE_BETTS_MILLER', # vs 'NONE' for no convective scheme
        'do_damping': True,
        'turb':True,
        'mixed_layer_bc':True,
        'do_virtual' :False,
        'do_simple': True,
        'roughness_mom':3.21e-05,
        'roughness_heat':3.21e-05,
        'roughness_moist':3.21e-05,                
    },

    'vert_turb_driver_nml': {
        'do_mellor_yamada': False,     # default: True
        'do_diffusivity': True,        # default: False
        'do_simple': True,             # default: False
        'constant_gust': 0.0,          # default: 1.0
        'use_tau': False
    },
    
    'diffusivity_nml': {
        'do_entrain':False,
        'do_simple': True,        #Affects computation of PBL depth on the unstable side
    },

    'surface_flux_nml': {
        'use_virtual_temp': False,
        'do_simple': True,
        'old_dtaudv': True    
    },

    'atmosphere_nml': {
        'idealized_moist_model': True              #Use False for Held-Suarez
    },

    'mixed_layer_nml': {
        'depth': 20,
        'albedo_value': 0.22,
        'tconst' : 285.,
        'prescribe_initial_dist':True,
        'evaporation':True,
        'do_qflux':False, #Don't use analytic formula for q-fluxes 
        'load_qflux': True, #Do load q-flux field from an input file
        'time_varying_qflux': False, #q-flux will be time-varying
        'qflux_file_name': f'{T0}_{T2}', #Name of q-flux input file
        'do_read_sst' : False, #Read in sst values from input file
        'do_sc_sst' : False, #Do specified ssts (need both to be true)
        'sst_file' : 'sst_ideal', #Set name of sst input file
    },

    'qe_moist_convection_nml': {     
        'rhbm':0.7,
        'Tmin':160.,
        'Tmax':350.   
    },
    
    'lscale_cond_nml': {
        'do_simple':True,
        'do_evap':False              #Reevaporation of rain
    },
    
    'sat_vapor_pres_nml': {
        'do_simple':True,
        'tcmin_simple':-250          #I changed that so thar C-K computations works in very cold climates
    },
    
    'damping_driver_nml': {
        'do_rayleigh': True,
        'trayfric': -0.5,              # neg. value: time in *days*
        'sponge_pbottom':  50.,
        'do_conserve_energy': True,         
    },

    'astronomy_nml': {
        'obliq': 0.0 # eccentricity zero by default
    },


    # FMS Framework configuration
    'diag_manager_nml': {
        'mix_snapshot_average_fields': False  # time avg fields are labelled with time in middle of window
    },

    'fms_nml': {
        'stack_size': 6000000,                               # default: 0
        'domains_stack_size': 6000000                        # default: 0
    },

    'fms_io_nml': {
        'threading_write': 'single',                         # default: multi
        'fileset_write': 'single',                           # default: multi
    },

     'spectral_dynamics_nml': {
        'damping_order': 4,             
        'water_correction_limit': 200.e2,
        'reference_sea_level_press':1.0e5,
        'num_levels':30,
        'valid_range_t':[50.,800.],
        'initial_sphum':[2.e-6],
        'vert_coord_option':'uneven_sigma',
        'surf_res':0.1,
        'scale_heights' : 5.0,
        'exponent':3.0,
        'robert_coeff':0.03
    }
    
})

[AndrewILWilliams]

A clarification, the Q-flux has a non-zero global mean value—I'm trying to force the slab model into some slightly weird climate states. One can think of this non-zero value as being a geothermal flux from the planet's interior. So really, the global-mean energy balance (in equilibrium) should be:

$SW_{\mathrm{TOA}}^{\downarrow}(1-\alpha)-OLR+Q_{\mathrm{flux}} = 0$

However, this equation has a residual that increases with the equator-to-pole temperature gradient! T2 is the second coefficient in the legendre polynomial expansion, and is related to the eq-Pole temp gradient.

I am now really confused...is there a state-dependent energy leak in the model? Or am I doing something silly?

[AndrewILWilliams]

I'm currently finding a similar magnitude TOA energy imbalance in the frierson test case (no qflux), where the only adjustment I've made is to set atm_abs=0

@ntlewis @sit23 and co, can you reproduce this?

[AndrewILWilliams]

Ran this out for ≈4000 days and the 2.5W/m2 persists

[AndrewILWilliams]

This implies an effective albedo of ≈0.318, as opposed to the prescribed albedo of 0.31 in the Frierson test case's mixed_layer_nml

[sit23]

Hi @AndrewILWilliams - I have checked energy conservation in the model on quite a few occasions and been reasonably satisfied, but I haven't got the data to hand just now. I'll try and take a look next week and see what I can find. In the meantime, could you share the code you're using to calculate the global average of these quantities? Thanks.

[AndrewILWilliams]

Thanks! Appreciate it.

def weighted_global_avg(da):
    weights = np.cos(np.deg2rad(da.lat))
    weighted = da.weighted(weights)
    weighted_mean = weighted.mean(['lat', 'lon'])
    return weighted_mean

[AndrewILWilliams]

For the frierson plot:

ds_tmp=xr.open_mfdataset(
            f"/scratch/cimes/aw2143/Isca/frierson_test_experiment/run*/atmos_monthly.nc", 
            decode_times=False
        )

fig, ax = plt.subplots()
weighted_global_avg(ds_tmp['swdn_toa']*(1-ds_tmp['albedo'])-ds_tmp['olr']).plot(ax=ax)
ax.set_ylim(-10,10)
ax.set_xlim(0,4000)
ax.axhline(y=0, c='grey', linestyle='--')
ax.set_title(r"Global mean of ''$SW^{\downarrow}(1-\alpha)-OLR$''"+"\n"+"Frierson test case, atm_abs=0")
ax.set_ylabel(r"W m$^{-2}$")
ax.set_xlabel("Time (days since run start)")
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)

[ntlewis]

@AndrewILWilliams, forgive me if you've included this and I somehow missed it, but you need weights = np.cos(lat)*np.diff(latb). This is because the grid points are not evenly spaced in latitude. latb should be in your Isca ouput. It contains the latitude points for the cell edges (so if you have N points in latitude, then you have N+1 edges).

[AndrewILWilliams]

Hi Neil, thanks! I just checked this and it only made a small difference—the imbalance in the Frierson run is still ≈2.5W/m2

Here's a dropbox link to the time average of the frierson run (once it reached steady state).

[AndrewILWilliams]

I mentioned this to Isaac and he suggested checking if all the do_conserve_energy flags were on, because 2.5W/m2 is about the magnitude of KE dissipation. It turns out I had it switched on in the namelist, but a postdoc at Princeton had hard-coded it to .false. in the call to vert_diff_init in idealized_moist_phys.F90 (for reasons which I don't understand).

I reran the frierson test case, fixing this, and now I have energy conservation! Yay!

Note to self: this also accounts for why the energy leak is state-dependent in my runs with differing eq-to-pole temperature gradients, because the storm track intensity (and presumably KE dissipation) is strengthening

Thanks both! Sorry for the false alarm. I'll pass on this bug fix to the princeton folks

[sit23]

Glad to hear you found the problem! Thanks @AndrewILWilliams!