What is the formula used for atmospheric visibility?

Question

I am working on NCEP GFS data to plot various maps. One such plot I would like to have is visibility. I have seen some of the sites showing visibility charts and the source of the data was shown as NCEP GFS.

I could not find any variable related to visibility in the GFS file.

Is there any specific formula available for calculating visibility?

Answer 1

I calculate surface visibility from WRF output using a calculation that I adapted from DTC's Unified Post Processor, specifically from their Fortran routine found in UPPV2.2/src/unipost/CALVIS.f. The calculation is based on hydrometeor mixing ratios, and air temperature and pressure, all from the lowest model layer. If your GFS output has hydrometeor mixing ratios, you can use this formula to calculate visibility. The documentation in the original code reads:

This routine computes horizontal visibility at the
surface or lowest model layer, from qc, qr, qi, and qs.  
qv--water vapor mixing ratio (kg/kg)
qc--cloud water mixing ratio (kg/kg)
qr--rain water mixing ratio  (kg/kg)
qi--cloud ice mixing ratio   (kg/kg)
qs--snow mixing ratio        (kg/kg)
tt--temperature              (k)
pp--pressure                 (Pa)

If iice=0:
   qprc=qr     qrain=qr and qclw=qc if T>0C
   qcld=qc          =0          =0  if T<0C
               qsnow=qs and qclice=qc  if T<0C
                    =0            =0   if T>0C
If iice=1:
   qprc=qr+qs   qrain=qr and qclw=qc
   qcld=qc+qi   qsnow=qs and qclice=qc

Independent of the above definitions, the scheme can use different
assumptions of the state of hydrometeors:
     meth='d': qprc is all frozen if T<0, liquid if T>0
     meth='b': Bocchieri scheme used to determine whether qprc
        is rain or snow. A temperature assumption is used to
        determine whether qcld is liquid or frozen.
     meth='r': Uses the four mixing ratios qrain, qsnow, qclw,
        and qclice

The routine uses the following
expressions for extinction coefficient, beta (in km**-1),
with C being the mass concentration (in g/m**3):

   cloud water:  beta = 144.7 * C ** (0.8800)
   rain water:   beta =  2.24 * C ** (0.7500)
   cloud ice:    beta = 327.8 * C ** (1.0000)
   snow:         beta = 10.36 * C ** (0.7776)

These expressions were obtained from the following sources:

   for cloud water: from Kunkel (1984)
   for rainwater: from M-P dist'n, with No=8e6 m**-4 and
      rho_w=1000 kg/m**3
   for cloud ice: assume randomly oriented plates which follow
      mass-diameter relationship from Rutledge and Hobbs (1983)
   for snow: from Stallabrass (1985), assuming beta = -ln(.02)/vis

The extinction coefficient for each water species present is
calculated, and then all applicable betas are summed to yield
a single beta. Then the following relationship is used to
determine visibility (in km), where epsilon is the threshhold
of contrast, usually taken to be .02:

   vis = -ln(epsilon)/beta      [found in Kunkel (1984)]

I have adapted the code from this routine to a Python function below, which you can use for your purposes:

def calculate_visibility(qv,qc,qr,qi,qs,T,p):
    """
    Calculates visibility based on the UPP algorithm.

    See documentation in UPPV2.2/src/unipost/CALVIS.f for the description of
    input arguments and references.
    """
    Rd = 287.
    COEFLC = 144.7
    COEFLP = 2.24
    COEFFC = 327.8
    COEFFP = 10.36
    EXPLC  = 0.88
    EXPLP  = 0.75
    EXPFC  = 1.
    EXPFP  = 0.7776

    Tv   = T * (1+0.61*qv) # Virtual temperature

    rhoa = p/(Rd*Tv) # Air density [kg m^-3]
    rhow = 1e3       # Water density [kg m^-3]
    rhoi = 0.917e3   # Ice density [kg m^-3]

    vovmd = (1+qv)/rhoa + (qc+qr)/rhow + (qi+qs)/rhoi

    conc_lc = 1e3*qc/vovmd
    conc_lp = 1e3*qr/vovmd
    conc_fc = 1e3*qi/vovmd
    conc_fp = 1e3*qs/vovmd

    # Make sure all concentrations are positive
    conc_lc[conc_lc < 0] = 0
    conc_lp[conc_lp < 0] = 0
    conc_fc[conc_fc < 0] = 0
    conc_fp[conc_fp < 0] = 0

    betav = COEFFC*conc_fc**EXPFC\
          + COEFFP*conc_fp**EXPFP\
          + COEFLC*conc_lc**EXPLC\
          + COEFLP*conc_lp**EXPLP+1E-10

    vis = -np.log(0.02)/betav # Visibility [km]
    vis[vis > 24.135] = 24.135

    return vis