{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Volcanic Ensemble SDE Simulation - Figure 7\n",
    "\n",
    "This notebook generates the final figure showing the impact of volcanic forcing on AMOC transitions in a simple stochastic differential equation (SDE) model. The analysis includes:\n",
    "\n",
    "1. **Age Test**: Three ensembles at different interstadial ages (t=12.2, 16.2, 18.2) with constant noise (σ=0.15)\n",
    "2. **Noise Level Test**: Three ensembles at the same interstadial age (t=18.2) with varying noise levels (σ=0.0, 0.07, 0.10)\n",
    "\n",
    "Each test includes 50 ensemble members comparing volcanic vs. non-volcanic scenarios.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.lines as mlines\n",
    "\n",
    "from scipy.integrate import odeint\n",
    "\n",
    "# widgets stuff\n",
    "import ipywidgets as widgets\n",
    "from ipywidgets import interact,interactive\n",
    "from ipywidgets import Checkbox\n",
    "from scipy.signal import argrelextrema\n",
    "from matplotlib.animation import FuncAnimation\n",
    "import os\n",
    "import sdeint\n",
    "import seaborn as sns \n",
    "from matplotlib import gridspec\n",
    "from matplotlib.patches import Rectangle\n",
    "\n",
    "import xarray as xr\n",
    "import pandas as pd\n",
    "\n",
    "from numba import jit\n",
    "from xclim import ensembles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def setaxesfsize(axl, fontsize):\n",
    "    \"\"\"Set font size for axis tick labels\"\"\"\n",
    "    for label in axl.get_xticklabels():\n",
    "        label.set_fontsize(fontsize)\n",
    "    for label in axl.get_yticklabels():\n",
    "        label.set_fontsize(fontsize)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Model parameters\n",
    "# This cell defines the parameters for the simple AMOC-AABW box model\n",
    "# The model represents the Atlantic Meridional Overturning Circulation (AMOC)\n",
    "# and Antarctic Bottom Water (AABW) as a coupled system of stochastic differential equations\n",
    "\n",
    "b0 = 0.625\n",
    "q0 = -9.\n",
    "q1 = 12.\n",
    "tau = 0.902\n",
    "tc = 166.66  # Characteristic time constant\n",
    "Bbar = 3.4164 \n",
    "eta = 0.5\n",
    "gammabar = 1.22\n",
    "\n",
    "def B0(gamma):\n",
    "    \"\"\"Calculate B0 based on gamma\"\"\"\n",
    "    return Bbar + eta * (gammabar - gamma)\n",
    "\n",
    "def dA(A, B, gamma, sigma, dt):\n",
    "    \"\"\"Stochastic differential equation for AMOC (A)\"\"\"\n",
    "    drift = -(B - B0(gamma)) - np.abs(q0 + q1*(A - b0))*(A - b0)\n",
    "    diffusion = sigma * np.random.normal()\n",
    "    return drift * dt + diffusion * np.sqrt(dt)\n",
    "\n",
    "def dB(A, gamma, sigma, dt):\n",
    "    \"\"\"Stochastic differential equation for AABW (B)\"\"\"\n",
    "    drift = (A - gamma) / tau\n",
    "    diffusion = sigma * np.random.normal()\n",
    "    return drift * dt + diffusion * np.sqrt(dt)\n",
    "\n",
    "def get_ensmem(sigma, gamma, volc=False, volct=15.0, volcdt=0.1, volcamp=0.005, \n",
    "               tlen=40.01, dt=0.001, t0=0, aabw0=1.0, amoc0=1.0, pert=0.0, vduration=0.07):\n",
    "    \"\"\"Generate a single ensemble member simulation\"\"\"\n",
    "    simtime = np.arange(t0, tlen, dt)\n",
    "    AMOC_nd = np.zeros(len(simtime))\n",
    "    AABW_nd = np.zeros(len(simtime))\n",
    "    \n",
    "    # Set initial conditions with optional perturbation\n",
    "    AABW_nd[0] = aabw0 + pert * (0.5 - np.random.rand())\n",
    "    AMOC_nd[0] = amoc0 + pert * (0.5 - np.random.rand())\n",
    "    \n",
    "    # Set up volcanic forcing\n",
    "    vtime = np.linspace(0, vduration, int(0.1/dt))\n",
    "    freq = 1000 * np.pi\n",
    "    volcseries = volcamp * np.cos(freq * vtime * np.pi / 180.)\n",
    "    \n",
    "    # Main integration loop\n",
    "    j = 0\n",
    "    for i, t in enumerate(simtime[:-1]):\n",
    "        if volc and volct <= t < volct + volcdt:\n",
    "            # Apply volcanic forcing\n",
    "            AMOC_nd[i+1] = AMOC_nd[i] + dA(AMOC_nd[i], AABW_nd[i], gamma, sigma, dt) + volcseries[j]\n",
    "            AABW_nd[i+1] = AABW_nd[i] + dB(AMOC_nd[i], gamma, sigma, dt)\n",
    "            j += 1\n",
    "        else:\n",
    "            # Normal evolution without volcanic forcing\n",
    "            AMOC_nd[i+1] = AMOC_nd[i] + dA(AMOC_nd[i], AABW_nd[i], gamma, sigma, dt)\n",
    "            AABW_nd[i+1] = AABW_nd[i] + dB(AMOC_nd[i], gamma, sigma, dt)\n",
    "    \n",
    "    return simtime, AMOC_nd, AABW_nd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Generating control run...\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Normalized AMOC Strength')"
      ]
     },
     "execution_count": 87,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 52664 characters]",
      "text/plain": [
       "<Figure size 1000x400 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Generate control run with noise to show realistic D-O oscillations\n",
    "# This creates a baseline simulation showing natural variability in the AMOC\n",
    "# The control run will be used as the background for all ensemble tests\n",
    "# Initial conditions are set to produce Dansgaard-Oeschger (D-O) type transitions\n",
    "\n",
    "print(\"Generating control run...\")\n",
    "\n",
    "gamma = 1.385\n",
    "\n",
    "novolcsim = get_ensmem(sigma=0.05, gamma=gamma, volc=False, tlen=200, dt=0.001, aabw0=2.5, amoc0=0.8, pert=0.0)\n",
    "\n",
    "# Create xarray Dataset for control run\n",
    "xr_novolcsim = xr.Dataset(\n",
    "    {\n",
    "        \"AMOC\": ((\"time\"), novolcsim[1]),\n",
    "        \"AABW\": ((\"time\"), novolcsim[2])\n",
    "    },\n",
    "    coords={\"time\": novolcsim[0]}\n",
    ")\n",
    "\n",
    "# we want the first random transitions to occur around t = 27 to 30\n",
    "# to make things agree with the complex model\n",
    "# just rerun until it is approximately there\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(10,4))\n",
    "xr_novolcsim.AMOC.plot(ax=ax)\n",
    "plt.xlim(0,40)\n",
    "plt.title(\"Control Run: AMOC Time Series with Stochastic Variability\")\n",
    "plt.xlabel(\"Time (arbitrary units)\")\n",
    "plt.ylabel(\"Normalized AMOC Strength\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Generating age test ensembles...\n",
      "Age test ensembles complete.\n",
      "Age test ensembles complete.\n"
     ]
    }
   ],
   "source": [
    "# Generate Age Test Ensembles (sigma = 0.15)\n",
    "# This cell creates three ensembles at different interstadial ages to test\n",
    "# how the timing of volcanic forcing affects AMOC transitions\n",
    "# Each ensemble has 50 members: one with volcanic forcing, one without\n",
    "# Volcanic amplitude: 0.006, Noise level: σ=0.15 (constant for all three)\n",
    "# Times tested: t=10.2, t=16.2, and t=26.2 (model time units)\n",
    "\n",
    "print(\"Generating age test ensembles...\")\n",
    "nruns = 50\n",
    "sig = 0.15\n",
    "gamma = 1.385\n",
    "#vamp = 0.005\n",
    "vamp = 0.006\n",
    "vduration = 0.07\n",
    "ipert = 1.e-2\n",
    "\n",
    "# Ensemble 1: Initial time at t=10.0\n",
    "it1 = 10000\n",
    "volc_init1 = 10.2\n",
    "tdur = 40\n",
    "x5novolc_t = xr_novolcsim.AMOC.time[it1]\n",
    "x5novolc_amoc = xr_novolcsim.AMOC[it1]\n",
    "x5novolc_aabw = xr_novolcsim.AABW[it1]\n",
    "\n",
    "ens_novolc_sim_1 = [get_ensmem(sig, gamma, volc=False, volct=volc_init1, volcdt=0.1, volcamp=vamp, \n",
    "                                tlen=tdur, dt=0.001, t0=x5novolc_t, aabw0=x5novolc_aabw, \n",
    "                                amoc0=x5novolc_amoc, pert=ipert, vduration=vduration) \n",
    "                    for i in range(nruns)]\n",
    "ens_volc_sim_1 = [get_ensmem(sig, gamma, volc=True, volct=volc_init1, volcdt=0.1, volcamp=vamp, \n",
    "                              tlen=tdur, dt=0.001, t0=x5novolc_t, aabw0=x5novolc_aabw, \n",
    "                              amoc0=x5novolc_amoc, pert=ipert, vduration=vduration) \n",
    "                  for i in range(nruns)]\n",
    "\n",
    "# Ensemble 2: Initial time at t=16.0\n",
    "it2 = 16000\n",
    "volc_init2 = 16.2\n",
    "tdur = 42\n",
    "x5novolc_t = xr_novolcsim.AMOC.time[it2]\n",
    "x5novolc_amoc = xr_novolcsim.AMOC[it2]\n",
    "x5novolc_aabw = xr_novolcsim.AABW[it2]\n",
    "\n",
    "ens_novolc_sim_2 = [get_ensmem(sig, gamma, volc=False, volct=volc_init2, volcdt=0.1, volcamp=vamp, \n",
    "                                tlen=tdur, dt=0.001, t0=x5novolc_t, aabw0=x5novolc_aabw, \n",
    "                                amoc0=x5novolc_amoc, pert=ipert, vduration=vduration) \n",
    "                    for i in range(nruns)]\n",
    "ens_volc_sim_2 = [get_ensmem(sig, gamma, volc=True, volct=volc_init2, volcdt=0.1, volcamp=vamp, \n",
    "                              tlen=tdur, dt=0.001, t0=x5novolc_t, aabw0=x5novolc_aabw, \n",
    "                              amoc0=x5novolc_amoc, pert=ipert, vduration=vduration) \n",
    "                  for i in range(nruns)]\n",
    "\n",
    "# Ensemble 3: Initial time at t=22.0\n",
    "#it3 = 18000\n",
    "it3 = 22000\n",
    "volc_init3 = 22.2\n",
    "tdur = 42\n",
    "x5novolc_t = xr_novolcsim.AMOC.time[it3]\n",
    "x5novolc_amoc = xr_novolcsim.AMOC[it3]\n",
    "x5novolc_aabw = xr_novolcsim.AABW[it3]\n",
    "\n",
    "ens_novolc_sim_3 = [get_ensmem(sig, gamma, volc=False, volct=volc_init3, volcdt=0.1, volcamp=vamp, \n",
    "                                tlen=tdur, dt=0.001, t0=x5novolc_t, aabw0=x5novolc_aabw, \n",
    "                                amoc0=x5novolc_amoc, pert=ipert, vduration=vduration) \n",
    "                    for i in range(nruns)]\n",
    "ens_volc_sim_3 = [get_ensmem(sig, gamma, volc=True, volct=volc_init3, volcdt=0.1, volcamp=vamp, \n",
    "                              tlen=tdur, dt=0.001, t0=x5novolc_t, aabw0=x5novolc_aabw, \n",
    "                              amoc0=x5novolc_amoc, pert=ipert, vduration=vduration) \n",
    "                  for i in range(nruns)]\n",
    "\n",
    "print(\"Age test ensembles complete.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_ensemble_xarray(ensemble_data_list):\n",
    "    \"\"\"\n",
    "    Create xarray Dataset for ensemble simulation data.\n",
    "    \n",
    "    Parameters:\n",
    "    ensemble_data_list: List of simulation outputs, each containing (time, AMOC, AABW)\n",
    "    \n",
    "    Returns:\n",
    "    xr.Dataset with dimensions for ensemble members and time\n",
    "    \"\"\"\n",
    "    time = ensemble_data_list[0][0]\n",
    "    num_members = len(ensemble_data_list)\n",
    "    \n",
    "    # Extract AMOC data for all ensemble members\n",
    "    ensemble_data = np.array([ensemble_data_list[i][1] for i in range(num_members)])\n",
    "    \n",
    "    return xr.Dataset(\n",
    "        {\"AMOC\": ((\"realization\", \"time\"), ensemble_data)},\n",
    "        coords={\"realization\": np.arange(num_members), \"time\": time}\n",
    "    )\n",
    "\n",
    "# Create xarray datasets for age test ensembles\n",
    "xr_ens_novolc_1 = create_ensemble_xarray(ens_novolc_sim_1)\n",
    "xr_ens_novolc_1_mean = xr_ens_novolc_1.mean(dim=\"realization\")\n",
    "xr_ens_volc_1 = create_ensemble_xarray(ens_volc_sim_1)\n",
    "xr_ens_volc_1_mean = xr_ens_volc_1.mean(dim=\"realization\")\n",
    "\n",
    "xr_ens_novolc_2 = create_ensemble_xarray(ens_novolc_sim_2)\n",
    "xr_ens_novolc_2_mean = xr_ens_novolc_2.mean(dim=\"realization\")\n",
    "xr_ens_volc_2 = create_ensemble_xarray(ens_volc_sim_2)\n",
    "xr_ens_volc_2_mean = xr_ens_volc_2.mean(dim=\"realization\")\n",
    "\n",
    "xr_ens_novolc_3 = create_ensemble_xarray(ens_novolc_sim_3)\n",
    "xr_ens_novolc_3_mean = xr_ens_novolc_3.mean(dim=\"realization\")\n",
    "xr_ens_volc_3 = create_ensemble_xarray(ens_volc_sim_3)\n",
    "xr_ens_volc_3_mean = xr_ens_volc_3.mean(dim=\"realization\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Generating noise level test ensembles...\n",
      "Noise level test ensembles complete.\n",
      "Noise level test ensembles complete.\n"
     ]
    }
   ],
   "source": [
    "# Generate Noise Level Test Ensembles (all at t=16.2 with different sigma)\n",
    "# This cell creates three ensembles at the same interstadial age (t=16.2) but with\n",
    "# different noise levels to test how background variability affects volcanic sensitivity\n",
    "# Each ensemble has 50 members: one with volcanic forcing, one without\n",
    "# Volcanic amplitude: 0.006 (same as age test), Time: t=16.2 (fixed)\n",
    "# Noise levels tested: σ=0.0 (no noise), σ=0.07 (moderate), σ=0.15 (high)\n",
    "\n",
    "print(\"Generating noise level test ensembles...\")\n",
    "nruns = 50\n",
    "gamma = 1.385\n",
    "# use same volcanic amplitude as age test ensembles\n",
    "#vamp = 0.0048\n",
    "#vamp = 0.007\n",
    "ipert = 1.e-2\n",
    "\n",
    "# Use initial condition at t=16.0 for noise level test\n",
    "itn = 16000\n",
    "volc_initn = 16.2\n",
    "x5novolc_t = xr_novolcsim.AMOC.time[itn]\n",
    "x5novolc_amoc = xr_novolcsim.AMOC[itn]\n",
    "x5novolc_aabw = xr_novolcsim.AABW[itn]\n",
    "\n",
    "# Noise level 1: sigma = 0.0\n",
    "sig = 0.0\n",
    "ens_novolc_sim_1n = [get_ensmem(sig, gamma, volc=False, volct=volc_initn, volcdt=0.1, volcamp=vamp, \n",
    "                                 tlen=40, dt=0.001, t0=x5novolc_t, aabw0=x5novolc_aabw, \n",
    "                                 amoc0=x5novolc_amoc, pert=ipert) for i in range(nruns)]\n",
    "ens_volc_sim_1n = [get_ensmem(sig, gamma, volc=True, volct=volc_initn, volcdt=0.1, volcamp=vamp, \n",
    "                               tlen=40, dt=0.001, t0=x5novolc_t, aabw0=x5novolc_aabw, \n",
    "                               amoc0=x5novolc_amoc, pert=ipert) for i in range(nruns)]\n",
    "\n",
    "# Noise level 2: sigma = 0.03\n",
    "#sig = 0.03\n",
    "sig = 0.07\n",
    "ens_novolc_sim_2n = [get_ensmem(sig, gamma, volc=False, volct=volc_initn, volcdt=0.1, volcamp=vamp, \n",
    "                                 tlen=40, dt=0.001, t0=x5novolc_t, aabw0=x5novolc_aabw, \n",
    "                                 amoc0=x5novolc_amoc, pert=ipert) for i in range(nruns)]\n",
    "ens_volc_sim_2n = [get_ensmem(sig, gamma, volc=True, volct=volc_initn, volcdt=0.1, volcamp=vamp, \n",
    "                               tlen=40, dt=0.001, t0=x5novolc_t, aabw0=x5novolc_aabw, \n",
    "                               amoc0=x5novolc_amoc, pert=ipert) for i in range(nruns)]\n",
    "\n",
    "# Noise level 3: sigma = 0.07\n",
    "#sig = 0.07\n",
    "sig = 0.15\n",
    "ens_novolc_sim_3n = [get_ensmem(sig, gamma, volc=False, volct=volc_initn, volcdt=0.1, volcamp=vamp, \n",
    "                                 tlen=40, dt=0.001, t0=x5novolc_t, aabw0=x5novolc_aabw, \n",
    "                                 amoc0=x5novolc_amoc, pert=ipert) for i in range(nruns)]\n",
    "\n",
    "ens_volc_sim_3n = [get_ensmem(sig, gamma, volc=True, volct=volc_initn, volcdt=0.1, volcamp=vamp, \n",
    "                               tlen=40, dt=0.001, t0=x5novolc_t, aabw0=x5novolc_aabw, \n",
    "                               amoc0=x5novolc_amoc, pert=ipert) for i in range(nruns)]\n",
    "\n",
    "print(\"Noise level test ensembles complete.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create xarray datasets for noise level test ensembles\n",
    "xr_ens_novolc_1n = create_ensemble_xarray(ens_novolc_sim_1n)\n",
    "xr_ens_novolc_1n_mean = xr_ens_novolc_1n.mean(dim=\"realization\")\n",
    "xr_ens_volc_1n = create_ensemble_xarray(ens_volc_sim_1n)\n",
    "xr_ens_volc_1n_mean = xr_ens_volc_1n.mean(dim=\"realization\")\n",
    "\n",
    "xr_ens_novolc_2n = create_ensemble_xarray(ens_novolc_sim_2n)\n",
    "xr_ens_novolc_2n_mean = xr_ens_novolc_2n.mean(dim=\"realization\")\n",
    "xr_ens_volc_2n = create_ensemble_xarray(ens_volc_sim_2n)\n",
    "xr_ens_volc_2n_mean = xr_ens_volc_2n.mean(dim=\"realization\")\n",
    "\n",
    "xr_ens_novolc_3n = create_ensemble_xarray(ens_novolc_sim_3n)\n",
    "xr_ens_novolc_3n_mean = xr_ens_novolc_3n.mean(dim=\"realization\")\n",
    "xr_ens_volc_3n = create_ensemble_xarray(ens_volc_sim_3n)\n",
    "xr_ens_volc_3n_mean = xr_ens_volc_3n.mean(dim=\"realization\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Note the model is very sensitive the control D-O cycle (i.e. whether we get a tipping event half way throught the cycle)\n",
    "#### This is necessary to get the right behavior as observed in Figure 7 from the paper.\n",
    "#### Noise plays a significant role in reproducing the results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plotting Functions\n",
    "# This cell defines helper functions for creating the final multi-panel figure:\n",
    "# - add_rectangle: Draws zoom region boxes on the overview plots\n",
    "# - calc_tprob: Calculates transition probability (% of ensemble members transitioning)\n",
    "# - plot_zoom: Creates detailed zoom plots with ensemble members, means, and histograms\n",
    "#\n",
    "# The plot_zoom function handles:\n",
    "#   * Plotting all 50 ensemble members (volcanic vs. non-volcanic)\n",
    "#   * Adding ensemble means and basin boundary markers\n",
    "#   * Creating histogram insets showing final AMOC distribution\n",
    "#   * Calculating and displaying transition probabilities\n",
    "\n",
    "def add_rectangle(ax, x_start, x_end, y_start, y_end, color, border_color=\"black\", border_width=0.7):\n",
    "    \"\"\"Add a rectangle to the given axis\"\"\"\n",
    "    rect_width = x_end - x_start\n",
    "    rect_height = y_end - y_start\n",
    "    rect = Rectangle((x_start, y_start), width=rect_width, height=rect_height, color=color)\n",
    "    border = Rectangle((x_start, y_start), width=rect_width, height=rect_height,\n",
    "                       facecolor=\"none\", edgecolor=border_color, linewidth=border_width)\n",
    "    ax.add_patch(rect)\n",
    "    ax.add_patch(border)\n",
    "\n",
    "def calc_tprob(data, threshold_val):\n",
    "    \"\"\"Calculate transition probability below threshold\"\"\"\n",
    "    data_series = pd.Series(data)\n",
    "    count_below = data_series[data_series < threshold_val].count()\n",
    "    return 100.0 * count_below / 50  # 50 ensemble members\n",
    "\n",
    "def plot_zoom(axtop, axzoom, axins, xr_ens_novolc, xr_ens_volc, xr_ens_novolc_mean, xr_ens_volc_mean,\n",
    "              eq_v, xzoom_start, xzoom_end, endtime, threshold_val, nruns, labelx=True, labely=True):\n",
    "    \"\"\"Plot zoomed ensemble with histogram inset\"\"\"\n",
    "    sns.set_style(\"darkgrid\")\n",
    "    \n",
    "    # Plot parameters\n",
    "    esalp = 0.6  # ensemble line transparency\n",
    "    esmalp = 0.9  # ensemble mean line transparency\n",
    "    eslw = 0.5  # ensemble line width\n",
    "    esmlw = 2.0  # ensemble mean line width\n",
    "    axzoomfs = 12\n",
    "    axinsfs = 12\n",
    "    axlbls = 14\n",
    "    ymin, ymax = 0.4, 1.8\n",
    "    \n",
    "    # Plot control run\n",
    "    xr_novolcsim.AMOC.plot(ax=axzoom, color=\"blue\", linewidth=1.0, alpha=0.7)\n",
    "    \n",
    "    # Plot ensemble members with smoothing\n",
    "    window = 6\n",
    "    # Plot first member of each ensemble with label, rest without\n",
    "    xr_ens_novolc.AMOC[0].rolling(time=window, center=True).mean().plot(\n",
    "        ax=axzoom, color=\"blue\", linewidth=eslw, alpha=esalp, label=\"Noise Ensemble\")\n",
    "    [xr_ens_novolc.AMOC[i].rolling(time=window, center=True).mean().plot(\n",
    "        ax=axzoom, color=\"blue\", linewidth=eslw, alpha=esalp) for i in range(1, nruns)]\n",
    "    \n",
    "    xr_ens_volc.AMOC[0].rolling(time=window, center=True).mean().plot(\n",
    "        ax=axzoom, color=\"red\", linewidth=eslw, alpha=esalp, label=\"Volcanic Ensemble\")\n",
    "    [xr_ens_volc.AMOC[i].rolling(time=window, center=True).mean().plot(\n",
    "        ax=axzoom, color=\"red\", linewidth=eslw, alpha=esalp) for i in range(1, nruns)]\n",
    "    \n",
    "    # Plot ensemble means\n",
    "    #xr_ens_novolc_mean.AMOC.plot(ax=axzoom, color=\"cyan\", linewidth=esmlw, alpha=esmalp, label=\"Noise (Ens Mean)\")\n",
    "    #xr_ens_volc_mean.AMOC.plot(ax=axzoom, color=\"yellow\", linewidth=esmlw, alpha=esmalp, label=\"Volc (Ens Mean)\")\n",
    "    \n",
    "    # Add basin boundary as transparent green shaded region\n",
    "    axzoom.axhspan(ymin=1.198, ymax=1.302, color='green', alpha=0.15, zorder=0, label='Basin Boundary')\n",
    "    \n",
    "    # Format zoom plot - set limits AFTER all plotting and turn off autoscale\n",
    "    axzoom.autoscale(False)\n",
    "    axzoom.set_xlim(xzoom_start, xzoom_end)\n",
    "    axzoom.set_ylim(ymin, ymax)\n",
    "    axzoom.vlines(eq_v, ymin, ymax, color=\"red\", alpha=0.7, linewidths=2.0, zorder=1)\n",
    "    setaxesfsize(axzoom, axzoomfs)\n",
    "    axzoom.set_xlabel(\"Model Time (non-dim)\" if labelx else \"\", size=axlbls)\n",
    "    axzoom.set_ylabel(\"AMOC Max (non-dim)\" if labely else \"\", size=axlbls)\n",
    "    axzoom.grid()\n",
    "    axzoom.legend(loc=3, fontsize=10, facecolor='#ffffff', framealpha=1.0)\n",
    "    axzoom.set_title(\"\")\n",
    "    \n",
    "    # Add volcanic event annotation\n",
    "    vlabxy = (0.254, 0.54)\n",
    "    vlabxyt = (0.02, 0.54)\n",
    "    axzoom.annotate('Volcanic\\nEvent', xy=vlabxy, xycoords='axes fraction', xytext=vlabxyt,\n",
    "                    textcoords='axes fraction', arrowprops=dict(facecolor='red', shrink=0.05),\n",
    "                    horizontalalignment='left', verticalalignment='top', size=12,\n",
    "                    bbox=dict(boxstyle=\"round,pad=0.3\", edgecolor=\"black\", facecolor='#ffffff', alpha=1.0))\n",
    "    \n",
    "    # Clean up tick labels\n",
    "    xticks = axzoom.get_xticks()\n",
    "    labels = [item.get_text() for item in axzoom.get_xticklabels()]\n",
    "    if len(labels) > 1:\n",
    "        labels[-1] = \"\"\n",
    "    axzoom.set_xticks(xticks)\n",
    "    axzoom.set_xticklabels(labels)\n",
    "    axzoom.tick_params(bottom=False, right=False)\n",
    "    axtop.set_title(\"\")\n",
    "    \n",
    "    # Inset histogram\n",
    "    axins.set_xlim(0, 20)\n",
    "    \n",
    "    xr_ens_novolc_lastyr = xr_ens_novolc.AMOC.isel(time=endtime)\n",
    "    xr_ens_volc_lastyr = xr_ens_volc.AMOC.isel(time=endtime)\n",
    "    \n",
    "    sns.histplot(ax=axins, y=xr_ens_novolc_lastyr, kde=True, orientation=\"vertical\",\n",
    "                 binwidth=0.03, binrange=[ymin, ymax], color=\"blue\", alpha=0.3, label=\"NoVolc\")\n",
    "    sns.histplot(ax=axins, y=xr_ens_volc_lastyr, kde=True, orientation=\"vertical\",\n",
    "                 binwidth=0.03, binrange=[ymin, ymax], color=\"red\", alpha=0.3, label=\"Volc\")\n",
    "    \n",
    "    # Format inset - set ylim AFTER plotting and turn off autoscale\n",
    "    axins.autoscale(False)\n",
    "    axins.set_ylim(ymin, ymax)\n",
    "    axins.tick_params(axis=\"x\", direction=\"in\", pad=-15)\n",
    "    axins.set_ylabel(\"Final AMOC Max\", size=axlbls)\n",
    "    axins.set_xlabel(\"Member Count\", size=axlbls)\n",
    "    axins.xaxis.set_label_coords(0.5, -0.04)\n",
    "    axins.yaxis.set_label_coords(0.7, 0.5)\n",
    "    axins.yaxis.set_label_position(\"right\")\n",
    "    axins.yaxis.tick_right()\n",
    "    axins.tick_params(axis=\"y\", direction=\"in\", pad=-25)\n",
    "    setaxesfsize(axins, axinsfs)\n",
    "    axins.grid()\n",
    "    axins.axhspan(ymin=1.198, ymax=1.302, color='green', alpha=0.15, zorder=0)\n",
    "    \n",
    "    # Add probability labels - volcanic on timeseries, noise on histogram\n",
    "    tprob_volc = calc_tprob(xr_ens_volc_lastyr, threshold_val)\n",
    "    tprob_novolc = calc_tprob(xr_ens_novolc_lastyr, threshold_val)\n",
    "    axzoom.text(0.36, 0.94, f'Volcanic Transitions: {tprob_volc:.0f}%', transform=axzoom.transAxes,\n",
    "               size=axinsfs, color='red', bbox=dict(boxstyle=\"round,pad=0.3\", edgecolor=\"black\",\n",
    "               facecolor='#ffffff', alpha=1.0), zorder=10)\n",
    "    axins.text(0.35, 0.94, f'Noise Transitions: {tprob_novolc:.0f}%', transform=axins.transAxes,\n",
    "               size=axinsfs, color='blue', bbox=dict(boxstyle=\"round,pad=0.3\", edgecolor=\"black\",\n",
    "               facecolor='#ffffff', alpha=1.0), zorder=10, ha='center')\n",
    "    \n",
    "    # Clean up inset tick labels\n",
    "    xticks = axins.get_xticks()\n",
    "    labels = [item.get_text() for item in axins.get_xticklabels()]\n",
    "    if len(labels) > 1:\n",
    "        labels[0] = labels[-1] = \"\"\n",
    "    axins.set_xticks(xticks)\n",
    "    axins.set_xticklabels(labels)\n",
    "    \n",
    "    yticks = axins.get_yticks()\n",
    "    labels = [item.get_text() for item in axins.get_yticklabels()]\n",
    "    if len(labels) > 1:\n",
    "        labels[0] = labels[-1] = \"\"\n",
    "    axins.set_yticks(yticks)\n",
    "    axins.set_yticklabels(labels)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 971036 characters]",
      "text/plain": [
       "<Figure size 1600x1600 with 14 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#### Main Figure Generation\n",
    "# This cell creates the complete multi-panel figure showing:\n",
    "#\n",
    "# LEFT SIDE (Age Test):\n",
    "#   - Panel A: Overview showing control run with three zoom regions\n",
    "#   - Panels B, C, D: Zoomed views at t=10.2, 16.2, 26.2 with σ=0.15\n",
    "#\n",
    "# RIGHT SIDE (Noise Level Test):\n",
    "#   - Panel E: Overview showing control run with one zoom region (t=16.2)\n",
    "#   - Panels F, G, H: Zoomed views at t=16.2 with σ=0.0, 0.07, 0.15\n",
    "#\n",
    "# Each zoom panel includes:\n",
    "#   * 50 ensemble members (blue=noise only, red=volcanic)\n",
    "#   * Ensemble means (cyan=noise, yellow=volcanic)\n",
    "#   * Histogram showing final AMOC distribution\n",
    "#   * Transition probabilities for both ensembles\n",
    "#\n",
    "# Grey backgrounds on zoom panels match grey boxes in overview plots\n",
    "\n",
    "fig = plt.figure(figsize=(16, 16))\n",
    "\n",
    "# set the style\n",
    "sns.set_style(\"ticks\")\n",
    "sns.set_context(\"paper\")\n",
    "\n",
    "# Create a 1x6 grid (this can be adjusted based on your needs)\n",
    "gs = gridspec.GridSpec(140, 200)\n",
    "\n",
    "# Three horizontal axes spanning the entire width\n",
    "zoomcolor = (0.8, 0.8, 0.8, 0.3)  # lightgray with 0.3 transparency, we are going to zoom out this area\n",
    "tlmid = 97 #vertical center between subplots\n",
    "trmid = 103\n",
    "ax0 = fig.add_subplot(gs[0:20, :tlmid])\n",
    "ax1 = fig.add_subplot(gs[26:60, 0:70], facecolor=zoomcolor)\n",
    "axins1 = fig.add_subplot(gs[26:60, 70:tlmid], facecolor=zoomcolor)\n",
    "ax2 = fig.add_subplot(gs[66:100, 0:70], facecolor=zoomcolor)\n",
    "axins2 = fig.add_subplot(gs[66:100, 70:tlmid], facecolor=zoomcolor)\n",
    "ax3 = fig.add_subplot(gs[106:140, 0:70], facecolor=zoomcolor)\n",
    "axins3 = fig.add_subplot(gs[106:140, 70:tlmid], facecolor=zoomcolor)\n",
    "\n",
    "# Duplicate axes on the right\n",
    "ax0_dup = fig.add_subplot(gs[0:20, trmid:])\n",
    "ax1_dup = fig.add_subplot(gs[26:60, trmid:173], facecolor=zoomcolor)\n",
    "axins1_dup = fig.add_subplot(gs[26:60, 173:200], facecolor=zoomcolor)\n",
    "ax2_dup = fig.add_subplot(gs[66:100, trmid:173], facecolor=zoomcolor)\n",
    "axins2_dup = fig.add_subplot(gs[66:100, 173:200], facecolor=zoomcolor)\n",
    "ax3_dup = fig.add_subplot(gs[106:140, trmid:173], facecolor=zoomcolor)\n",
    "axins3_dup = fig.add_subplot(gs[106:140, 173:200], facecolor=zoomcolor)\n",
    "\n",
    "ax = [ax0, ax1, axins1, ax2, axins2, ax3, axins3, ax0_dup, ax1_dup, axins1_dup, ax2_dup, axins2_dup, ax3_dup, axins3_dup]\n",
    "\n",
    "################## Main Plot\n",
    "lw = 1.5\n",
    "lalp = 0.7\n",
    "ax0ymax = 2.0\n",
    "ax0ymin = 0.\n",
    "axlbls = 14  # axes label size\n",
    "axts = 12  # axes tick label size\n",
    "\n",
    "################ Define volcanic event times (CHANGE ONLY THESE!)\n",
    "# These should match the values used in ensemble generation\n",
    "volc_time_1 = 10.2  # corresponds to it1 = 10000\n",
    "volc_time_2 = 16.2  # corresponds to it2 = 16000\n",
    "volc_time_3 = 22.2  # corresponds to it3 = 22000\n",
    "volc_time_n = 16.2  # corresponds to itn = 16000\n",
    "\n",
    "################ Define zoom regions (automatically calculated from volcanic times)\n",
    "zoom_width = 1.6  # Width of zoom window (symmetric around volcanic event)\n",
    "xzoom1_start = volc_time_1 - 0.4\n",
    "xzoom1_end = volc_time_1 + 1.2\n",
    "xzoom2_start = volc_time_2 - 0.4\n",
    "xzoom2_end = volc_time_2 + 1.2\n",
    "xzoom3_start = volc_time_3 - 0.4\n",
    "xzoom3_end = volc_time_3 + 1.2\n",
    "xzoom3_dup_start = volc_time_n - 0.4\n",
    "xzoom3_dup_end = volc_time_n + 1.2\n",
    "\n",
    "# Add vertical lines for volcanic events (automatically from defined times)\n",
    "vgrid = [volc_time_1, volc_time_2, volc_time_3]\n",
    "rec_ymin = 0.4\n",
    "rec_ymax = 1.8\n",
    "[ax0.vlines(vgrid[i], ax0ymin, ax0ymax, color=\"red\", alpha=0.7, linewidths=2.0, zorder=3)\n",
    " for i in [0,1,2]]\n",
    "add_rectangle(ax0, xzoom1_start, xzoom1_end, rec_ymin, rec_ymax, zoomcolor)\n",
    "add_rectangle(ax0, xzoom2_start, xzoom2_end, rec_ymin, rec_ymax, zoomcolor)\n",
    "add_rectangle(ax0, xzoom3_start, xzoom3_end, rec_ymin, rec_ymax, zoomcolor)\n",
    "\n",
    "vgrid_dup = [volc_time_n]\n",
    "[ax0_dup.vlines(vgrid_dup[i], ax0ymin, ax0ymax, color=\"red\", alpha=0.7, linewidths=2.0, zorder=3)\n",
    " for i in [0]]\n",
    "add_rectangle(ax0_dup, xzoom3_dup_start, xzoom3_dup_end, rec_ymin, rec_ymax, zoomcolor)\n",
    "\n",
    "cntcol = \"black\"\n",
    "cntalp = 0.5\n",
    "# volcanic run colors\n",
    "vcols = [\"purple\", \"blue\", \"red\"]\n",
    "# eruption colors\n",
    "ecols = [\"magenta\", \"cyan\", \"orange\"]\n",
    "\n",
    "# Define threshold value before using it\n",
    "threshold_val = 1.18\n",
    "nruns = 50\n",
    "endtime = 1400  # this is in 0.001 dt units\n",
    "\n",
    "# Plot control run\n",
    "xr_novolcsim.AMOC.plot(\n",
    "    ax=ax0,\n",
    "    color=\"blue\",\n",
    "    linewidth=1.0,\n",
    "    alpha=0.7,\n",
    "    label=\"Control: No Volcanism\",\n",
    "    zorder=2\n",
    ")\n",
    "\n",
    "ax0.axhline(\n",
    "    1.388,\n",
    "    0,\n",
    "    1,\n",
    "    color=\"k\",\n",
    "    alpha=0.5,\n",
    "    linewidth=1.0,\n",
    "    linestyle=\"dashed\",\n",
    "    label=\"Upper Saddle Node\",\n",
    ")\n",
    "ax0.axhline(\n",
    "    0.97,\n",
    "    0,\n",
    "    1,\n",
    "    color=\"k\",\n",
    "    alpha=0.5,\n",
    "    linewidth=1.0,\n",
    "    linestyle=\"dotted\",\n",
    "    label=\"Lower Saddle Node\",\n",
    ")\n",
    "# Add basin boundary as transparent green shaded region\n",
    "ax0.axhspan(1.198, 1.302, color='green', alpha=0.15, zorder=0, label=\"Basin Boundary\")\n",
    "\n",
    "# Duplicate the main plot on the right\n",
    "ax0_dup.plot(xr_novolcsim.time, xr_novolcsim.AMOC, color=\"blue\", linewidth=1.0, alpha=0.7, label=\"Control: No Volcanism\", zorder=2)\n",
    "ax0_dup.axhline(1.388, 0, 1, color=\"k\", alpha=0.5, linewidth=1.0, linestyle=\"dashed\", label=\"Upper Saddle Node\")\n",
    "ax0_dup.axhline(0.97, 0, 1, color=\"k\", alpha=0.5, linewidth=1.0, linestyle=\"dotted\", label=\"Lower Saddle Node\")\n",
    "# Add basin boundary as transparent green shaded region\n",
    "ax0_dup.axhspan(1.198, 1.302, color='green', alpha=0.15, zorder=0, label=\"Basin Boundary\")\n",
    "\n",
    "# Create zoom plots (using automatically defined times)\n",
    "zoomplot_1 = plot_zoom(ax0, ax1, axins1, xr_ens_novolc_1, xr_ens_volc_1, xr_ens_novolc_1_mean, xr_ens_volc_1_mean,\n",
    "            volc_time_1, xzoom1_start, xzoom1_end, endtime, threshold_val, nruns, labelx=True, labely=True)\n",
    "zoomplot_1_dup = plot_zoom(ax0_dup, ax1_dup, axins1_dup, xr_ens_novolc_1n, xr_ens_volc_1n, xr_ens_novolc_1n_mean, xr_ens_volc_1n_mean,\n",
    "            volc_time_n, xzoom3_dup_start, xzoom3_dup_end, endtime, threshold_val, nruns, labelx=True, labely=False)\n",
    "\n",
    "zoomplot_2 = plot_zoom(ax0, ax2, axins2, xr_ens_novolc_2, xr_ens_volc_2, xr_ens_novolc_2_mean, xr_ens_volc_2_mean,\n",
    "            volc_time_2, xzoom2_start, xzoom2_end, endtime, threshold_val, nruns, labelx=True, labely=True)\n",
    "zoomplot_2_dup = plot_zoom(ax0_dup, ax2_dup, axins2_dup, xr_ens_novolc_2n, xr_ens_volc_2n, xr_ens_novolc_2n_mean, xr_ens_volc_2n_mean,\n",
    "            volc_time_n, xzoom3_dup_start, xzoom3_dup_end, endtime, threshold_val, nruns, labelx=True, labely=False)\n",
    "\n",
    "zoomplot_3 = plot_zoom(ax0, ax3, axins3, xr_ens_novolc_3, xr_ens_volc_3, xr_ens_novolc_3_mean, xr_ens_volc_3_mean,\n",
    "            volc_time_3, xzoom3_start, xzoom3_end, endtime, threshold_val, nruns, labelx=True, labely=True)\n",
    "zoomplot_3_dup = plot_zoom(ax0_dup, ax3_dup, axins3_dup, xr_ens_novolc_3n, xr_ens_volc_3n, xr_ens_novolc_3n_mean, xr_ens_volc_3n_mean,\n",
    "            volc_time_n, xzoom3_dup_start, xzoom3_dup_end, endtime, threshold_val, nruns, labelx=True, labely=False)\n",
    "\n",
    "# Format left overview plot (ax0)\n",
    "setaxesfsize(ax0, axts)\n",
    "ax0.set_ylabel(\"AMOC Max (non-dim)\", size=axlbls)\n",
    "ax0.set_xlabel(\"Model Time (non-dim)\", size=axlbls)\n",
    "ax0.set_title(\"Simple Model Interstadial Age Test\", size=14, loc=\"center\")\n",
    "ax0.set_xlim(0, 40)\n",
    "ax0.set_ylim(ax0ymin, ax0ymax)\n",
    "ax0.legend(loc=3, fontsize=10, ncol=2, facecolor='white', framealpha=1.0)\n",
    "\n",
    "# Format right overview plot (ax0_dup)\n",
    "setaxesfsize(ax0_dup, axts)\n",
    "ax0_dup.set_xlabel(\"Model Time (non-dim)\", size=axlbls)\n",
    "ax0_dup.set_title(\"Simple Model Noise Level Test\", size=14, loc=\"center\")\n",
    "ax0_dup.set_xlim(0, 40)\n",
    "ax0_dup.set_ylim(ax0ymin, ax0ymax)\n",
    "ax0_dup.legend(loc=3, fontsize=10, ncol=2, facecolor='white', framealpha=1.0)\n",
    "ax0_dup.set_ylabel(\"\")\n",
    "\n",
    "# title the zoom plots\n",
    "## add size and position variables\n",
    "zoomtitlesize = 12\n",
    "zoomtitlepos = \"left\"\n",
    "ax1.set_title(r\"$\\sigma$ = 0.15\", size=zoomtitlesize, loc=zoomtitlepos)\n",
    "ax2.set_title(r\"$\\sigma$ = 0.15\", size=zoomtitlesize, loc=zoomtitlepos)\n",
    "ax3.set_title(r\"$\\sigma$ = 0.15\", size=zoomtitlesize, loc=zoomtitlepos)\n",
    "ax1_dup.set_title(r\"$\\sigma$ = 0.00\", size=zoomtitlesize, loc=zoomtitlepos)\n",
    "ax2_dup.set_title(r\"$\\sigma$ = 0.07\", size=zoomtitlesize, loc=zoomtitlepos)\n",
    "ax3_dup.set_title(r\"$\\sigma$ = 0.15\", size=zoomtitlesize, loc=zoomtitlepos)\n",
    "\n",
    "# Annotating my plots:\n",
    "flsize=20\n",
    "xannleft = 0.01\n",
    "yannleft = 0.86\n",
    "ax0.annotate('A', \n",
    "             xy=(xannleft, yannleft), \n",
    "             xycoords='axes fraction', \n",
    "             fontsize=flsize, \n",
    "             weight='bold'\n",
    ")\n",
    "yannleft = 0.91\n",
    "ax1.annotate('B', \n",
    "             xy=(xannleft, yannleft), \n",
    "             xycoords='axes fraction', \n",
    "             fontsize=flsize, \n",
    "             weight='bold'\n",
    ")\n",
    "ax2.annotate('C', \n",
    "             xy=(xannleft, yannleft), \n",
    "             xycoords='axes fraction', \n",
    "             fontsize=flsize, \n",
    "             weight='bold' \n",
    ") \n",
    "ax3.annotate('D', \n",
    "             xy=(xannleft, yannleft), \n",
    "             xycoords='axes fraction', \n",
    "             fontsize=flsize, \n",
    "             weight='bold' \n",
    ")\n",
    "## annotate ax0 - position labels at volcanic event times (automatically calculated)\n",
    "flax0size=10\n",
    "# Convert volcanic event times to axes fraction coordinates\n",
    "# ax0 xlim is 0-40, add small offset to shift right of the vertical line\n",
    "offset = 0.3  # shift 0.3 time units to the right\n",
    "xb = (volc_time_1 + offset) / 40.0  # it1 timing\n",
    "yb = 0.8\n",
    "ax0.annotate('B', \n",
    "             xy=(xb, yb), \n",
    "             xycoords='axes fraction', \n",
    "             fontsize=flax0size, \n",
    "             weight='bold'\n",
    ")\n",
    "xc = (volc_time_2 + offset) / 40.0  # it2 timing\n",
    "yc = 0.8\n",
    "ax0.annotate('C', \n",
    "             xy=(xc, yc), \n",
    "             xycoords='axes fraction', \n",
    "             fontsize=flax0size, \n",
    "             weight='bold'\n",
    ")\n",
    "xd = (volc_time_3 + offset) / 40.0  # it3 timing\n",
    "yd = 0.8\n",
    "ax0.annotate('D', \n",
    "             xy=(xd, yd), \n",
    "             xycoords='axes fraction', \n",
    "             fontsize=flax0size, \n",
    "             weight='bold'\n",
    ")\n",
    "xannleft = 0.01\n",
    "yannleft = 0.86\n",
    "ax0_dup.annotate('E', \n",
    "             xy=(xannleft, yannleft), \n",
    "             xycoords='axes fraction', \n",
    "             fontsize=flsize, \n",
    "             weight='bold'\n",
    ")\n",
    "yannleft = 0.91\n",
    "ax1_dup.annotate('F', \n",
    "             xy=(xannleft, yannleft), \n",
    "             xycoords='axes fraction', \n",
    "             fontsize=flsize, \n",
    "             weight='bold'\n",
    ")\n",
    "ax2_dup.annotate('G', \n",
    "             xy=(xannleft, yannleft), \n",
    "             xycoords='axes fraction', \n",
    "             fontsize=flsize, \n",
    "             weight='bold' \n",
    ") \n",
    "ax3_dup.annotate('H', \n",
    "             xy=(xannleft, yannleft), \n",
    "             xycoords='axes fraction', \n",
    "             fontsize=flsize, \n",
    "             weight='bold' \n",
    ") \n",
    "## annotate ax0_dup - position labels at volcanic event time (automatically calculated)\n",
    "flax0size=10\n",
    "# ax0_dup xlim is also 0-40, all noise tests at same volcanic time (itn)\n",
    "offset_dup = 0.3  # shift 0.3 time units to the right\n",
    "xf = (volc_time_n + offset_dup) / 40.0  # itn timing for F\n",
    "yf = 0.6\n",
    "ax0_dup.annotate('F', \n",
    "             xy=(xf, yf), \n",
    "             xycoords='axes fraction', \n",
    "             fontsize=flax0size, \n",
    "             weight='bold'\n",
    ")\n",
    "xg = (volc_time_n + offset_dup) / 40.0  # itn timing for G\n",
    "yg = 0.45\n",
    "ax0_dup.annotate('G', \n",
    "             xy=(xg, yg), \n",
    "             xycoords='axes fraction', \n",
    "             fontsize=flax0size, \n",
    "             weight='bold'\n",
    ")\n",
    "xh = (volc_time_n + offset_dup) / 40.0  # itn timing for H\n",
    "yh = 0.32\n",
    "ax0_dup.annotate('H', \n",
    "             xy=(xh, yh), \n",
    "             xycoords='axes fraction', \n",
    "             fontsize=flax0size, \n",
    "             weight='bold'\n",
    ")\n",
    "\n",
    "# Save figure\n",
    "icepaperdir = \"/home/guido/Dropbox/plots/volcanic/\"\n",
    "saveplot = False\n",
    "if saveplot:\n",
    "    plt.savefig(\n",
    "        icepaperdir + \"Figure7_volc_perturb_simple_model_v2.png\",\n",
    "        dpi=300,\n",
    "        bbox_inches=\"tight\",\n",
    "        facecolor=fig.get_facecolor(),\n",
    "        transparent=False,\n",
    "    )\n",
    "    plt.savefig(\n",
    "        icepaperdir + \"Figure7_volc_perturb_simple_model_v2.pdf\",\n",
    "        dpi=300,\n",
    "        bbox_inches=\"tight\",\n",
    "        facecolor=fig.get_facecolor(),\n",
    "        transparent=False,\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (geosci)",
   "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
