Table of Contents

Introduction

Obukhov Length is a well-known length scale in boundary layer meteorology. It characterizes stability effects (buoyancy versus shear turbulence) in the turbulent layer of the atmosphere near the surface. It is used in air pollution modelling and to characterize the shape of wind profiles in wind energy applications (among many other applications).

The following procedure describes how the Obukhov length (L) can be calculated from ERA5 data.

The formula used in the code, computes the inverse of the Obukhov length (1/L [m-1]) as eq (5.7c) from Stull (1988): 

$$1/L = \frac{vk \ast g \ast tvst}{(tv2 \ast ust^2)}$$

where vk = 0.4 is the Von Karman constant, g is the gravitational acceleration, tvst is the turbulent temperature scale, tv2 is the virtual temperature and ust is the friction velocity.

The output is a 2D field from input data taken from the Climate Data Store (CDS) ERA5 hourly data on single levels from 1959 to present dataset. 

You will need:

Step 1: Download input data from the CDS

The variables that are needed from the ERA5 data, are:

  • 'surface_pressure',
  • '2m_dewpoint_temperature',
  • '2m_temperature',
  • 'instantaneous_eastward_turbulent_surface_stress',
  • 'instantaneous_moisture_flux',
  • 'instantaneous_northward_turbulent_surface_stress',
  • 'instantaneous_surface_sensible_heat_flux',
  • 'standard_deviation_of_filtered_subgrid_orography'


Please be careful about the ECMWF ERA5 sign convention for fluxes. All downward fluxes are positive, so during daytime the surface sensible heat flux (sshf) is typically negative.

The Python script used the CDS API to download the ERA5 data from the CDS. The suggested procedure is:

  1. Copy the script below to a text editor on your computer
  2. Edit the 'year', 'month', 'day', 'time', and output filename ('data_2020010112.grib' in this example) in the script to meet your requirements. Note: 'grid' and 'area' keywords can also be included.
  3. Save the script (for example with the filename as 'get_data.py')
  4. Run the script: python3 get_data.py

get_data
# **************************** LICENSE START ***********************************
#
# Copyright 2022 ECMWF. This software is distributed under the terms
# of the Apache License version 2.0. In applying this license, ECMWF does not
# waive the privileges and immunities granted to it by virtue of its status as
# an Intergovernmental Organization or submit itself to any jurisdiction.
#
# ***************************** LICENSE END ************************************
import cdsapi

c = cdsapi.Client()

c.retrieve('reanalysis-era5-single-levels', {
    'product_type': 'reanalysis',
    'format': 'grib',
    'variable': [
        'surface_pressure',
        '2m_dewpoint_temperature',
        '2m_temperature',
        'instantaneous_eastward_turbulent_surface_stress',
        'instantaneous_moisture_flux',
        'instantaneous_northward_turbulent_surface_stress',
        'instantaneous_surface_sensible_heat_flux',
        'standard_deviation_of_filtered_subgrid_orography'
    ],
    'year': '2020',
    'month': '01',
    'day': '01',
    'time': '12:00',
}, 'data_2020010112.grib')

Running the script produces a file in the current working directory called 'data_2020010112.grib' . This is a GRIB file containing the ERA5 variables needed. 

Step 2: Compute inverse Obukhov length 1/L

The suggested procedure to run the Python script to compute the inverse Obukhov length (1/L) is :

  1. Copy the script below to a text editor
  2. Save the script as 'calculate-Obukhov-length.py'
  3. Run the script 'calculate-Obukhov-length.py' with the correct arguments: python3 calculate-Obukhov-length.py data_2020010112.grib output.grib
calculate-Obukhov-length
# **************************** LICENSE START ***********************************
#
# Copyright 2022 ECMWF. This software is distributed under the terms
# of the Apache License version 2.0. In applying this license, ECMWF does not
# waive the privileges and immunities granted to it by virtue of its status as
# an Intergovernmental Organization or submit itself to any jurisdiction.
#
# ***************************** LICENSE END ************************************

import sys

import numpy as np

from eccodes import codes_count_in_file, codes_get, codes_get_values
from eccodes import codes_clone, codes_set, codes_set_values
from eccodes import codes_set_key_vals, codes_write
from eccodes import codes_grib_new_from_file, codes_release
from eccodes import CodesInternalError

def qswat(t, p):
    """
    Computes saturation q (with respect to water)
    In  t   : Temperature                   (K)
        p   : Pressure                      (Pa)
    Out qswt: Saturation specific humidity  (kg/kg)
    """
    rkbol = 1.380658e-23
    rnavo = 6.0221367e+23
    r = rnavo * rkbol
    rmd = 28.9644
    rmv = 18.0153
    rd = 1000 * r / rmd # Gas constant for air
    rv = 1000 * r / rmv # Gas constant for water vapour
    restt = 611.21
    r2es = restt * rd / rv
    r3les = 17.502
    r4les = 32.19
    retv = rv / rd - 1  # Tv = T(1+retv*q)
    rtt = 273.16        # Melting point (0 Celcius)
    foeew = r2es * np.exp((r3les * (t - rtt)) / (t - r4les))
    qs = foeew / p
    zcor = 1 / (1 - retv * qs)
    qs = qs * zcor

    return qs

def calculate_Obukhov_length(data_file, output_file):
    try:
        gid_list = None
        with open(data_file) as f_in:
            mcount = codes_count_in_file(f_in)
            gid_list = [codes_grib_new_from_file(f_in) for i in range(mcount)]

            # Input data: sp, 2d, 2t, ie, ishf, iews, inss
            data = {}
            for i in range(mcount):
                gid = gid_list[i]
                sn = codes_get(gid, 'shortName')
                vals = codes_get_values(gid) # Type: class 'numpy.ndarray'
                data[sn] = vals

            rd = 287.06     # Gas constant
            retv = 0.6078   # Tv = T*(1+retv*q)
            cp = 1004.7     # Air heat capacity at constant pressure
            lmelt = 2500800 # Latent heat of evaporation
            vk = 0.4        # VonKarman constant
            g = 9.81        # Gravity acceleration

            q2 = qswat(data['2d'], data['sp'])  # q2 is saturation value at 2d
            tv2 = data['2t'] * (1 + retv * q2)  # Virtual temperature
            rho = data['sp'] / (rd * tv2)       # Air density
            tau = np.sqrt(data['iews']**2 + data['inss']**2)
                                                # Turb. surface stress
            ust = np.maximum(np.sqrt(tau / rho), 0.001) # Friction velocity

            wt = -data['ishf'] / (rho * cp)     # Turbulent heat flux
            wq = -data['ie'] / rho              # Turbulent moisture flux
            wtv = wt + retv * data['2t'] * wq   # Virtual turbulent heat flux
            tvst = -wtv / ust                   # Turbulent temperature scale 
            qst = -wq / ust                     # Turbulent moisture scale
            Linv = vk * g * tvst / (tv2 * ust**2)   # Inverse Obukhov length
                                                # Stull eq. (5.7c)
            Linvrange = 100.    # Range of useful Linv values to control
                                # grib-coding
            Linv = np.maximum(np.minimum(Linv, Linvrange), -Linvrange)
                                                             
            print('Linv max/min/mean: %e %e %e' %
                (Linv.max(), Linv.min(), np.abs(Linv).mean()))
	    
            # Write to a file
            with open(output_file, 'wb') as f_out:
                gid = codes_clone(gid_list[0])
                codes_set(gid, 'paramId', 80) # Experimental parameter
                codes_set_values(gid, Linv)
                codes_write(gid, f_out)
                codes_release(gid)
    except CodesInternalError as ex:
        sys.stderr.write('Error: %s\n' % ex)
        sys.exit(-100)
    finally:
        if gid_list != None:
            for gid in gid_list:
                codes_release(gid)

def print_usage():
    print('python calculate-Obukhov-length.py data.grib output.grib')

def main(argv):
    if len(argv) != 2:
        print_usage()
        sys.exit(100)

    calculate_Obukhov_length(argv[0], argv[1])

if __name__ == '__main__':
    main(sys.argv[1:])

Running the script produces a file in the current working directory called 'output.grib' . This is a GRIB file containing the inverse of the Obukhov length. 

If you require the output data in netCDF, the GRIB file with the results can be converted to netCDF on your local system by using the ecCodes command grib_to_netcdf

Limitations 

In areas with significant orography, the inverse Obukhov length (1/L), as computed above, may be less usefull for two reasons: (a) Monin-Obkhov similarity becomes questionable in areas with orography, and (b) the ERA5 momentum fluxes, as needed for the computation of 1/L, contain the sum of turbulent flux and momentum flux due to unresolved small scale orography. The latter is parametrized by the so-called Turbulent Orographic Form Drag scheme (TOFD). To identify such areas, the retrieval job also includes parameter sdfor (standard deviation of filtered subgrid orography). It is recommended to use 1/L only in areas where sdfor < 50 m. In that case the contribution of TOFD to the total turbulent stress is small (Beljaars et al. 2004). 

References

Stull, Roland B., 1988: An Introduction to Boundary Layer Meteorology. Springer, 670 pp.

Beljaars, A., Brown, A. and Wood, N., 2004: A new parametrization of turbulent orographic form drag. QJRMS, 130, pp.1327-1347.

This document has been produced in the context of the Copernicus Climate Change Service (C3S).

The activities leading to these results have been contracted by the European Centre for Medium-Range Weather Forecasts, operator of C3S on behalf of the European Union (Delegation Agreement signed on 11/11/2014 and Contribution Agreement signed on 22/07/2021). All information in this document is provided "as is" and no guarantee or warranty is given that the information is fit for any particular purpose.

The users thereof use the information at their sole risk and liability. For the avoidance of all doubt , the European Commission and the European Centre for Medium - Range Weather Forecasts have no liability in respect of this document, which is merely representing the author's view.

Related articles