\n $\"No$ \n

Timeseries of FLEXPART and measured eBC\u00b6

In this example, FLEXPART modelled black carbon (BC) are compared with measured equivalent black carbon (eBC) at Birkenes. eBC might exist as a combined product in the ACTRIS data portal, if not eBC will be derieved from aerosol light absorption coefficients. The calculation of equivalent black carbon is shown in notebook \"Equivalent black carbon (eBC).ipynb\", and is the same for both combined products (https://vocabulary.actris.nilu.no/skosmos/actris_vocab/en/page/combinedproduct) (in this case datasets with pre-calculated equivalent black carbon) and observation products (https://vocabulary.actris.nilu.no/skosmos/actris_vocab/en/page/observation) (datasets with hourly values of the standard ACTRIS variable aerosol particle light absorption coefficient). In both cases eBC is calculated from the absorption coefficient at 637 nm, devided by the MAC conversion factor. If the data does not contain the reference wavelength, it is interpolated with the \u00c5ngstr\u00f8m exponent.

This example is organized as follows:

(1) \u00a0\u00a0 Define all python libraries needed

(2) \u00a0\u00a0 Specify input files, define variables, and set the month to be studied/plotted

(3) \u00a0\u00a0 Read FLEXPART data and make a simple test plot

(4) \u00a0\u00a0 Read ACTRIS data and check the content, e.g., absorption coefficients

(5) \u00a0 Check if eBC exists as a level 3 product. If not, calculate \u00c5ngst\u00f8m coefficients ($\\alpha$) and absorption coefficients at 637nm (if needed)

\u2003\u2003\u2003\u2003 $\\alpha=ln(\\frac{abs_{470nm}}{abs_{660nm}}) /ln(\\frac{660}{470}) $

\u2003\u2003\u2003\u2003 $abs_{637nm}={abs_{660nm}}(\\frac{637}{660})^{-\\alpha} $

\u2003\u2003\u2003\u2003 Plot the \u00c5ngstr\u00f8m coefficients as a test

\u2003\u2003 And calculate Equivalent Black carbon (eBC)

\u2003\u2003\u2003\u2003 $eBC_{637nm}=\\frac{abs_{637nm}}{MAC_{637nm}}$

(6) \u00a0\u00a0 Plot eBC from FLEXPART and measurements for one month

1. Import libraries\u00b6

In\u00a0[1]:

import numpy as np\nimport pandas as pd\nimport xarray as xr\nimport plotly.express as px\nimport re\nimport requests\nimport matplotlib.pyplot as plt\nfrom io import StringIO\n

2. Define file names and variables\u00b6

In\u00a0[2]:

# FLEXPART file (from ACTRIS data portal, Services: ATMO-ACCESS Flexpart)\n\n#Flex_file = 'https://atmo-access.nilu.no/WEBPLOTS/NO0002R/2022/OUT_202201_0/ECMWF/sourcespectrum_TOTAL_BC_a_OUT_202201_0.dat' # file that contains the BC sources for the month of January 2022\nFlex_file = 'https://atmo-access.nilu.no/WEBPLOTS/NO0002R/2023/OUT_202302_0/ECMWF/sourcespectrum_TOTAL_BC_a_OUT_202302_0.dat' # file that contains the BC sources for the month of February 2023\n\n# Observation file (from ACTRIS data portal) \n#Obs_file = 'https://thredds.nilu.no/thredds/dodsC/ebas_doi/EJ/DP/Q2/EJDP-Q2TQ.nc' # file that contains the data for the month of January 2022, this is a *observation* file, i.e. contains aerosol absorption coefficient measurements\nObs_file = 'https://thredds.nilu.no/thredds/dodsC/ebas_doi/KC/Q9/7U/KCQ9-7U9Q.nc' # file that contains data for the month of February 2023, this is a *combined* file, i.e. contains pre-calculated eBC concentrations\n\n# Constants for calculating eBC\nlambda_ref=637.\nMAC=10\n\n# Month to be analyzed, change this to the month you want to analyze\n#month='2022-01'\nmonth='2023-02' \n

3. Read FLEXPART data and make a test plot\u00b6

In\u00a0[3]:

# FLEXPART-file (Flex_file) is defined above. Send a GET request to the URL\nresponse = requests.get(Flex_file)\n\n# Check if the request was successful\nif response.status_code == 200:\n    # Read the content into a pandas DataFrame\n    data = StringIO(response.text)\n    df = pd.read_csv(data, sep=r'\\s+')\n    \n    # Display the first few rows to understand its structure\n    print(df.head())\n    \n    # Select date/time column\n    time_column = df.columns[1]  \n    \n    # Convert the time column to datetime\n    df[time_column] = pd.to_datetime(df[time_column], format='%Y%m%d%H%M%S')\n    \n    # Select various eBC sources\n    ebc_flexpart = 'eBC_Flexpart'\n    value_column1 = df.columns[5]  # DOM\n    value_column2 = df.columns[6]  # ENE\n    value_column3 = df.columns[7]  # FLR\n    value_column4 = df.columns[8]  # IND\n    value_column5 = df.columns[9]  # SHP\n    value_column6 = df.columns[10]  # WST\n    value_column7 = df.columns[11]  # TRA\n    value_column8 = df.columns[12]  # BB\n    # Total eBC from Flexpart\n    df[ebc_flexpart] = df[[value_column1, value_column2, value_column3, value_column4, value_column5, \n                           value_column6, value_column7, value_column8]].sum(axis=1)\n\n    # Plot the time series to check\n    plt.figure(figsize=(10, 4))\n    plt.plot(df[time_column], df[value_column1])\n    plt.plot(df[time_column], df[value_column2])\n    plt.plot(df[time_column], df[value_column3])\n    plt.plot(df[time_column], df[value_column4])\n    plt.plot(df[time_column], df[value_column5])\n    plt.plot(df[time_column], df[value_column6])\n    plt.plot(df[time_column], df[value_column7])\n    plt.plot(df[time_column], df[value_column8])\n    #plt.plot(df[time_column], df[ebc_flexpart])\n    \n    plt.xlabel('Date')\n    plt.ylabel('ng/m3')\n    month_label = pd.to_datetime(month).strftime('%B %Y')\n    plt.title(f'FLEXPART eBC sources for Birkenes in {month_label}')\n    plt.xticks(rotation=40)\n    plt.legend(['DOM', 'ENE', 'FLR', 'IND', 'SHP', 'WST', 'TRA', 'BB'])\n    plt.grid(True)\n    plt.show()\n\nelse:\n    print(f\"Failed to download the file. Status code: {response.status_code}\")\n

    UNIT           START             END    LON      LAT      DOM     ENE  \\\n0  ng/m3  20230201000000  20230201030000  8.252  58.3885  17.2900  0.2538   \n1  ng/m3  20230201030000  20230201060000  8.252  58.3885  23.4763  0.2789   \n2  ng/m3  20230201060000  20230201090000  8.252  58.3885  20.5765  0.3193   \n3  ng/m3  20230201090000  20230201120000  8.252  58.3885  13.8812  0.1859   \n4  ng/m3  20230201120000  20230201150000  8.252  58.3885  13.6053  0.1451   \n\n      FLR     IND     SHP     WST     TRA      BB  \n0  1.9440  1.7208  0.3720  0.4953  6.6918  1.6699  \n1  2.1129  2.2984  0.3422  0.6202  9.4350  1.3629  \n2  3.1455  2.0917  0.3651  0.5499  8.1086  1.0050  \n3  3.5289  1.4182  0.3136  0.3279  5.6985  0.4056  \n4  4.3181  1.3316  0.3571  0.2441  5.4814  0.1620  \n

\n $\"No$ \n

4. Read ACTRIS data\u00b6

In\u00a0[4]:

# Open the dataset \ninst_data = xr.open_dataset(Obs_file) \n\nall_wavelengths = inst_data.Wavelength.values\ntime=inst_data.time.values \n\nprint(\"Wavelengths:\",all_wavelengths)\n\n# Find the wavelength closest to the reference wavelength and the shortest wavelength\nlambda_long = min(all_wavelengths, key=lambda x:abs(x-lambda_ref))\nlambda_short = min(all_wavelengths)\n\nprint(\"Wavelength closest to ref:\",lambda_long)\nprint(\"Shortest wavelength:\",lambda_short)\n

Wavelengths: [637.]\nWavelength closest to ref: 637.0\nShortest wavelength: 637.0\n

5. Calculate \u00c5ngstr\u00f8m exponent from observations, calculate eBC\u00b6

In\u00a0[5]:

# Check if eBC exists as a level 3 product. If not, calculate eBC from absorption coefficients\n\nif 'equivalent_black_carbon_amean' in inst_data.variables:\n    # Read eBC level 3 product and convert to dataframe\n    df_eBC_ref = inst_data.equivalent_black_carbon_amean.rename('eBC_measured').to_dataframe().reset_index().set_index('time')\n    print(\"eBC exists as a level 3 product\")   \nelse:\n    # Select absorption coefficients at two wavelengths defined above\n    abs_short = inst_data.sel(Wavelength=lambda_short).aerosol_absorption_coefficient_amean\n    abs_long = inst_data.sel(Wavelength=lambda_long).aerosol_absorption_coefficient_amean\n\n    # Create empty arrays with length of dataset\n    A = len(abs_long)\n    abs_coeff = np.zeros(A)     #Abs. coeff. at reference wavelength\n    eBC_ref = np.zeros(A)       #eBC at reference wavelength\n    alpha = np.zeros(A)         #\u00c5ngstr\u00f6m exponent\n\n    # Calculate \u00c5ngstr\u00f6m exponent (alpha)\n    if lambda_long == lambda_ref:\n        alpha = 0.0\n    else:\n        alpha = (np.log(abs_short/abs_long)) / (np.log(lambda_long/lambda_short))\n        alpha[alpha < 0.1] = 0.7\n        alpha[alpha > 3.0] = 3.0\n        \n        #Plot the \u00c5ngstr\u00f6m exponent\n        plt.figure(figsize=(10, 4))\n        plt.plot(time,alpha)\n\n\n    # Calculate absorption coefficient and eBC at reference wavelength\n    abs_coeff = ((lambda_ref/lambda_long)**(-alpha))*abs_long\n    eBC_ref = (abs_coeff /MAC).rename('eBC_measured')\n    \n    # Convert to Dataframe\n    df_eBC_ref = eBC_ref.to_dataframe()\n

eBC exists as a level 3 product\n

6. Plot eBC from observations and Flexpart\u00b6

In\u00a0[6]:

# Extract the time period specified above\ndf_hour_ref_month = df_eBC_ref.eBC_measured[month]\nx=df_hour_ref_month.index[:]\n\n# Plot option 1 (not visible in pdf):\n#####################################\n# Create station name: a)selct the first part of the name, b)remove non-alphanumeric characters \nstation_name_full = inst_data.ebas_station_name\nstation_string = inst_data.ebas_station_name.split(\" \",1)[0]\nstation_name = re.sub(r'\\W+','', station_string)\n\n# Create plot with specified size\nfig=px.line(df_hour_ref_month)\nfig.add_scatter(x=df[time_column], y=df[ebc_flexpart]/1000, name=\"eBC_Flexpart\", xaxis='x', yaxis='y')\n\nfig.update_layout(\n    title=\"Equivalent Black Carbon (eBC) at Birkenes\",\n    xaxis_title=\"Time\",\n    yaxis_title=\"\u00b5g/m3\",\n    font=dict(\n        family=\"Courier New, monospace\",\n        size=16,\n        color=\"RebeccaPurple\"\n    ),\n    width=1500,  # Specify the width of the plot\n    height=500   # Specify the height of the plot\n)\nfig.show() \n

To see a more detailed calculation of equvalent black carbon from measurement aerosol absorption coefficient data, please have a look at \"Equivalent black carbon (eBC).ipynb\".

To have a look at footprint emission sensitivity analysis from the FLEXPART model, please have a look at \"FLEXPART footprint emission for black carbon.ipynb\".