SE_model.py

#load modules
import matplotlib.pyplot as plt
import matplotlib
import numpy as np
import pandas as pd
import sys
import numba
from scipy.optimize import minimize

sys.path.append('../')


@numba.njit
def semi_empirical_sea_level_model(T, T0_init, a, c0, tau1, tau2, dt=1.0):
    """
    Simulates the sea-level rise using the semi-empirical model.

    Parameters:
    T (array): Array of global mean temperatures (in degrees Celsius).
    T0_init (float): Initial equilibrium temperature (in degrees Celsius).
    a (float): Sensitivity of sea level rate to temperature deviation (in mm/yr/K).
    c0 (float): Initial temperature-independent rate term (in mm/yr).
    tau1 (float): Timescale for temperature relaxation (in years).
    tau2 (float): E-folding time for the temperature-independent rate term (in years).
    dt (float): Time step (in years). Default is 1.0 year.

    Returns:
    h (array): Simulated global sea level (in mm).
    T0 (array): Simulated equilibrium temperatures (in degrees Celsius).
    c (array): Simulated temperature-independent rate term (in mm/yr).
    dh_dt (array): Simulated rate of sea level rise (in mm/yr).
    """

    
    n = len(T)
    T0 = np.zeros(n)
    c = np.zeros(n)
    dh_dt = np.zeros(n)
    h = np.zeros(n)

    T0[0] = T0_init
    c[0] = c0
    dh_dt[0] = c0
    
    
    for t in range(1, n):
        if isinstance(dt, (int, float)):
            T0[t] = T0[t-1] + (T[t-1] - T0[t-1]) / tau1 * dt
            c[t] = c[t-1] - c[t-1] / tau2 * dt
        else:
            T0[t] = T0[t-1] + (T[t-1] - T0[t-1]) / tau1 * dt[t-1]
            c[t] = c[t-1] - c[t-1] / tau2 * dt[t-1]
    
    dh_dt[1:] = a * (T[1:] - T0[1:]) + c[1:]
    h[1:] = np.cumsum(dh_dt[1:] * dt) + h[0]
    return h, T0, c, dh_dt


def find_optimum_y_offset(ts1, ts2, ts_1_std):
    """
    Find the optimum Y value offset between two time series considering the uncertainty.

    Parameters:
    ts1 (np.ndarray): Time series with uncertainty (normal distribution).
    ts2 (np.ndarray): Time series without uncertainty.
    uncertainty (np.ndarray): Standard deviation of the uncertainties for ts1.

    Returns:
    float: The optimum Y value offset.
    """
    assert len(ts1) == len(ts2) == len(ts_1_std), "All input arrays must have the same length."

    # Define the objective function to minimize (weighted sum of squared residuals)
    def objective(offset):
        residuals = (ts1 + offset - ts2) / ts_1_std
        return np.sum(residuals**2)

    # Initial guess for the offset
    initial_guess = 0
    # Use scipy.optimize.minimize to find the optimal offset
    result = minimize(objective, initial_guess)
    
    return result.x[0]

@numba.njit
def calculate_C_at_year(C_500, tau2, target_year, reference_year=500):
    """
    Calculate C at a specific year given C at a reference year and tau2.

    Parameters:
    C_500 (float): Value of C at the reference year (500 CE).
    tau2 (float): E-folding time for the temperature-independent rate term (in years).
    target_year (float): The year for which to calculate C.
    reference_year (float): The reference year (default is 500 CE).

    Returns:
    float: Value of C at the target year.
    """
    return C_500 * np.exp(-(target_year - reference_year) / tau2)

@numba.njit
def linear_extrapolation(mean_tem_age, sea_level, x):
    # Extract the last two points
    x1, x2 = mean_tem_age[-2], mean_tem_age[-1]
    y1, y2 = sea_level[-2], sea_level[-1]
    
    # Calculate the slope (rate of change)
    slope = (y2 - y1) / (x2 - x1)
    
    # Extrapolate to the desired x value
    y_extrapolated = y2 + slope * (x - x2)
    
    return y_extrapolated