{ "cells": [ { "cell_type": "markdown", "id": "02b47774", "metadata": {}, "source": [ "# Fractal Aggregate Aerosols in a Hot Jupiter atmosphere\n", "\n", "In this tutorial, we will learn:\n", "\n", "1. How to run virga for non-spherical cloud particles. \n", "2. How to compare the particle sizes and optical properties for aggregates vs. spheres\n", "\n", "Aerosols can come in a variety of shapes and sizes, which we typically characterise by their fractal dimension (Df). Under this framework, aerosols are made of collections of small spheres called monomers, and Df=3 would represent the most compact possible arrangement, whereas Df=1 would be something perfectly linear. Real aerosols are usually somewhere in-between 1 and 3, and the ultimate shape of an aerosol is usually a result of their formation conditions. VIRGA v2.0 allows us to silmultaneously model the effects of a range of different particle shapes on spectra so we can find the best fit to observational data. For further details and example of typical shapes, see Lodge & Moran et al. (2025).\n", "\n", "Striving to maintain the simplicity and ethod of virga, aerosols can be added by providing two simple parameters:\n", "\n", " a) The fractal dimension (general shape) of the aerosols (which can be an array).\n", " b) EITHER the number of monomers (N_mon) OR the radius of the monomer (r_mon).\n", "\n", "Figure 2 of Moran & Lodge et al. (2025) demonstrates both of the choices for (b) and explains them in detail -- see this for more guidance. in this tutorial we will choose an array of fractal dimensions, and we choose to fix the number of monomers ($N_{mon}$=1000).\n", "\n", "Before beginning, make sure you have downloaded the new optical database for aggregates from here: https://zenodo.org/records/16581692 (which are all calculated with the assumption that $N_{mon}$=1000) and update the filepaths to the VIRGA reference data in the first cell:" ] }, { "cell_type": "code", "execution_count": null, "id": "9fb3885c", "metadata": {}, "outputs": [], "source": [ "# --------------------------------------------- UPDATE THESE FILEPATHS FOR YOUR OWN SETUP ------------------------------------------------------------------\n", "\n", "# first, set the locations of databases needed for virga\n", "mieff_dir = '/Users/nbatalh1/Documents/data/VIRGA_2_mieff_files' # set file paths for VIRGA v2 optical database (.mieff files) -- make sure you have downloaded the new database that includes aggregates (v2!)\n", "\n", "# --------------------------------------------------------------------------------------------------------------------------------------------" ] }, { "cell_type": "code", "execution_count": null, "id": "944f9007", "metadata": {}, "outputs": [], "source": [ "#standard import statements\n", "import numpy as np\n", "import pandas as pd\n", "import astropy.units as u\n", "from virga import justdoit as vdi\n", "\n", "#plotting tools\n", "from virga import justplotit as vpi\n", "from bokeh.plotting import show, figure\n", "from bokeh.io import output_notebook\n", "output_notebook()\n", "import matplotlib.pyplot as plt " ] }, { "cell_type": "markdown", "id": "aa4a7158", "metadata": {}, "source": [ "## Set planet parameters\n", "\n", "To start with, let's create choose a standard hot Jupiter profile and a condensible species for the clouds. Here we choose Mg2SiO4, a common species for these environments." ] }, { "cell_type": "code", "execution_count": null, "id": "33c11958", "metadata": {}, "outputs": [], "source": [ "# set planet parameters: \n", "\n", "planet_mass = 1 # M_Jup\n", "planet_radius = 1 # R_Jup\n", "\n", "# choose condensible species\n", "aggregate_species = ['Mg2SiO4'] # Aggregate species\n", "\n", "# set Kzz, metallicity and mean molecular weight \n", "metallicity = 1 #atmospheric metallicity relative to Solar\n", "mean_molecular_weight = 2.3 # atmospheric mean molecular weight\n", "constant_kzz = 10**10 # constant Kzz\n", "f_sed = 0.5 # sedimentation efficiency\n", "sigma_width = 2 # width of lognormal distribution\n", "\n", "# create hot Jupiter atmosphere\n", "jupiter = vdi.Atmosphere(aggregate_species, \n", " fsed=f_sed,\n", " mh=metallicity,\n", " mmw = mean_molecular_weight, \n", " sig = sigma_width,\n", " aggregates=False)\n", "\n", "# set gravity\n", "jupiter.gravity(mass=planet_mass, mass_unit=vdi.u.Unit('M_jup'), radius=planet_radius, radius_unit=vdi.u.Unit('R_jup')) # order: planet mass, unit, planet radius, unit\n", "\n", "# load preset hot jupiter p-T profile\n", "jupiter.ptk(df = vdi.hot_jupiter())\n", "\n", "# add constant Kzz (value set by user at top of this cell) to atmosphere profile\n", "jupiter.kz = np.zeros(len(jupiter.kz))+ constant_kzz # make a constant Kzz profile as an array" ] }, { "cell_type": "markdown", "id": "070aa3e5", "metadata": {}, "source": [ "## Adding Mg2SiO4 Clouds (Spherical Particles)" ] }, { "cell_type": "code", "execution_count": null, "id": "d9b200e0", "metadata": {}, "outputs": [], "source": [ "# run virga: this will determine where the clouds are, what particles sizes they are made from, and how opaque they are at each pressure level\n", "clouds_from_virga_spheres = vdi.compute(jupiter, directory=mieff_dir) # Add the Mg2SiO4 clouds!" ] }, { "cell_type": "markdown", "id": "1ffb7f72", "metadata": {}, "source": [ "## p-T profile\n", "\n", "To view the p-T profile and get an idea of where in the atmosphere our Mg2SiO4 condensates are forming into clouds, we can use virga's usual pt() function." ] }, { "cell_type": "code", "execution_count": null, "id": "40b7bc81", "metadata": {}, "outputs": [], "source": [ "# Plot pressure-Temperature profile\n", "show(vpi.pt(clouds_from_virga_spheres,plot_height=450))" ] }, { "cell_type": "markdown", "id": "a624c429", "metadata": {}, "source": [ "## Adding Fractal Aggregate (Mg2SiO4) Clouds.\n", "\n", "Everything up to this point has been the same method as for Virga v1.0. Now let's add fractal aggregate clouds to our clear model instead of spheres. The function is very similar, but we just provide two extra bits of information (fractal dimension and number of monomers). We have added a loop below to make clouds for five different fractal dimensions and save the output to a list." ] }, { "cell_type": "code", "execution_count": null, "id": "14d2195b", "metadata": {}, "outputs": [], "source": [ "# Set a range of fractal dimensions to explore (a different model will be run each time). This can be a single value or an array. \n", "# In this example, we will choose the values [1.2, 1.6, 2.0, 2.4, 2.8]:\n", "d_f_values = np.linspace(1.2,2.8,5)\n", "\n", "# create a list to store all of the different cloud models for each fractal dimension\n", "clouds_from_virga_fractals=[]\n", "\n", "# loop over all particle shapes\n", "for i in range(len(d_f_values)):\n", "\n", "\n", " d_f = d_f_values[i]\n", " N_monomers = 1000 # Here we provide the number of monomers. If you prefer to provide r_mon instead, set it here (but note that it needs to be in cm, not um), and submit r_mon instead of N_mon in the .virga() function below\n", " \n", " jupiter.aggregates=True\n", " jupiter.Df= d_f\n", " jupiter.N_mon = N_monomers\n", " # That's it -- Df and N_mon are all you need to prescribe for aggregates! More details on other typical aggregate parameters below for experts:\n", " # k0 will usually be calculated using Eq. 14+15 of Moran & Lodge et al. (2025), but this can be prescribed here if you wish to use a fixed value.\n", " # Radius of gyration will be calculated within the code: virga can calculate the mass of the each aggregate, meaning that r_mon can be calculated\n", " # from N_mon (or vice versa). k0 is calculated as stated above, so including Df we have everything we need to find R_gyr using Eq. 2 of Moran & Lodge \n", " # et al. (2025)\n", "\n", " # take the clear atmosphere (cloudy_jupiter_list[i]) and run virga to find the clouds for each fractal dimension (using aggregates=True and our aggregate parameters). \n", " # this method saves each cloud model as a new model in a list (clouds_from_virga_fractals).\n", " clouds_from_virga_fractals.append(vdi.compute(jupiter, directory=mieff_dir)) # Add the Mg2SiO4 clouds!" ] }, { "cell_type": "markdown", "id": "f3ff32d1", "metadata": {}, "source": [ "# Aggregate Analysis Tools\n", "\n", "We've made our clouds! Now let's see how the aggregates compared to spheres. There are a number of tools that have been added to compare between the different models, and each one can give you interesting diagnostic information. Let's explore them one-by-one. First, set-up the figure size and the min/max wavelength to plot." ] }, { "cell_type": "code", "execution_count": null, "id": "b0597dad", "metadata": {}, "outputs": [], "source": [ "# set figure size in notebook\n", "plt.rcParams['figure.figsize'] = [15, 10]\n", "\n", "# set min and max wavelength to show on plots\n", "min_wavelength = 0.3 # in um\n", "max_wavelength = 15 # in um\n", "\n", "# set custom colors for plots\n", "cloud_colors = ['#ce3fce', '#7cf2fb', '#33b9cc', '#1c7bb8', '#295091', '#212255'] # order: spheres, 1.2, 1.6, 2.0, 2.4, 2.8" ] }, { "cell_type": "markdown", "id": "80550169", "metadata": {}, "source": [ "## Visualising Particle Sizes in each Pressure Layer\n", "\n", "Virga will typically find different particle sizes for different aggregate shapes because of how they behave dynamically (see Moran & Lodge et al. 2025). Lower fractal dimensions tend to form into larger particles (because they have slower fall speed vs compact spheres. However, the upwards forces are the same in both versions of the model. Therefore, the aggregates need to form into more massive particles to stay stable against the same upwards velocity from Kzz/convection). The mean particle sizes that virga found for each shape can be seen and compared using the tool below. " ] }, { "cell_type": "code", "execution_count": null, "id": "48f1c63d", "metadata": {}, "outputs": [], "source": [ "vpi.aggregates_pressure_vs_radius(clouds_from_virga_spheres, clouds_from_virga_fractals, d_f_list=d_f_values, colors=cloud_colors)" ] }, { "cell_type": "markdown", "id": "1426452b", "metadata": {}, "source": [ "\n", "Interpretation: The Df=1.2 model has the largest particles, which are at the bottom of the cloud-deck. In all models, cloud particles get smaller as altitude increases." ] }, { "cell_type": "markdown", "id": "ff27cb13", "metadata": {}, "source": [ "## Probing the relationship between Pressure, Radius and Number Density\n", "\n", "Because the same amount of condensate is condensing out in each model in virga (regardless of particle shape), if the low fractal dimensions form into larger particles, there must be fewer of them. Plotting the number density at each pressure layer shows demonstrates this balance between particle size and number density visually with altitude. The marker size indicates particle size. " ] }, { "cell_type": "code", "execution_count": null, "id": "52ca3460", "metadata": {}, "outputs": [], "source": [ "vpi.aggregates_pressure_vs_number_density(clouds_from_virga_spheres, clouds_from_virga_fractals, d_f_list=d_f_values, colors=cloud_colors)" ] }, { "cell_type": "markdown", "id": "c8bf4e5f", "metadata": {}, "source": [ "Interpretation: The Df=1.2 clouds are shown to be made of fewer (but larger) particles than in the other models, just like in Figure 4 of Lodge & Moran et al. (2025)." ] }, { "cell_type": "markdown", "id": "5aed211d", "metadata": {}, "source": [ "## Visualising .mieff files and Optical Properties\n", "\n", "This tool allows you to see how the optical properties of the cloud particles vary with shape. $Q_{ext}$ represents how efficiently they will absorb and scatter radiation, and it is a function of particle size, wavelength and shape. This reads all of the data from the .mieff database and plots $Q_{ext}$ for all particle sizes (blue = smallest particles in the grid, green = largest particles in the grid).\n" ] }, { "cell_type": "code", "execution_count": null, "id": "34549e44", "metadata": {}, "outputs": [], "source": [ "# plot Q_ext as a function of wavelength for each fractal dimension\n", "vpi.aggregates_optical_properties(aggregate=aggregate_species[0], mieff_dir=mieff_dir, d_f_list=d_f_values, min_wavelength=min_wavelength, max_wavelength=max_wavelength)" ] }, { "cell_type": "markdown", "id": "2d9145b9", "metadata": {}, "source": [ "Interpretation: This is usually a busy plot, but it is useful because it can provide intuition into how linear and compact aggregates compare to their spherical counterparts. Here the linear fractals with Df = 1.2 are the most opaque at short wavelengths. We saw that these particles are the least opaque in the transmission spectrum, which is because their number density is so much lower as shown in the previous figure.\n", "\n", "These plots can highlight resonance features and aid in checking that your optics are behaving as expected for a particular chemical species." ] }, { "cell_type": "markdown", "id": "fc8990e6", "metadata": {}, "source": [ "## Checking that Particle Radii are within the .mieff Grid\n", "\n", "It's important that virga finds particles that are within the range of optical properties in the .mieff grid, otherwise virga may underpredict or overpredict the opacity of the particles. Virga will give warnings if it finds particles off-grid, but this plot is a quick visual check to see if and where that happens:" ] }, { "cell_type": "code", "execution_count": null, "id": "37636d6b", "metadata": {}, "outputs": [], "source": [ "vpi.aggregates_wavelength_radius_grid(clouds_from_virga_spheres, clouds_from_virga_fractals, aggregate=aggregate_species[0], d_f_list=d_f_values, mieff_dir=mieff_dir)" ] }, { "cell_type": "markdown", "id": "ec2239e0", "metadata": {}, "source": [ "Interpretation: \n", "- The red outline box marks the extremes of the .mieff grid, with faint red vertical lines representing the grid spacing in radius.\n", "- The vertical extent of the box is just to show the wavelength range in the .mieff files. The colored markers show the column density so that we have an indication of how many particles were in each radius bin for each model.\n", "- The \"radius\" of an aggregate is defined as the compact radius (the radius of a sphere made from the same volume of material) and this is how they are stored in the .mieff files. The particle sizes obtained by virga at each pressure layer and for each fractal dimension are plotted on this grid as markers, and these should all ideally be inside the horizontal bounds of the red box. If not, you may need to adjust your Kzz/fsed values (both of these can effect particle size), or make a new .mieff grid (see next tutorial!). It's ok for the particles to be near the vertical bounds of the red box -- the marker y-positions just scale with column density, so they will always sit on the top/bottom edge of the box." ] }, { "cell_type": "markdown", "id": "23ddc1a2-976b-48bc-aade-8aa6b3437e87", "metadata": {}, "source": [ "# Next steps" ] }, { "cell_type": "markdown", "id": "66d245a5-d59b-46ff-b125-79433495fc93", "metadata": {}, "source": [ "To see how different the spectra would be for these fractal aggergate clouds, check out the next tutorial: \"11 Fractal Aggregate Spectra in PICASO\"" ] } ], "metadata": { "kernelspec": { "display_name": "pic312", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.12.2" }, "toc": { "base_numbering": 1, "nav_menu": {}, "number_sections": true, "sideBar": true, "skip_h1_title": false, "title_cell": "Table of Contents", "title_sidebar": "Contents", "toc_cell": true, "toc_position": {}, "toc_section_display": true, "toc_window_display": false } }, "nbformat": 4, "nbformat_minor": 5 }