Convolution/weighted average of two arrays in Python

Question

I have an equation that I need to calculate numerically, but I am having doubts about my approach. I am cross-posting this question from Stack Exchange, because I am not getting any responses.

This is the equation I need to program:

where $\chi$ is a volume mixing ratio in the stratosphere, $\chi_0$ is a volume mixing ratio in the troposphere and G is an age spectrum (a distribution of times it takes an air parcel to travel from troposphere to stratosphere). G is an inverse Gaussian:

where t' is defined as the transit time $\Delta$t divided by $\Gamma$. The tropospheric volume mixing ratios as a function of time form an almost linear curve. I'll put the curves below (ignore the dashed lines). The time lags are with respect to 2009, i.e. a time lag of 1 corresponds to 2008 etc. Now I need to apply the first equation to these two curves numerically. I have tried numpy.convolve(chi_0[::-1], G, mode='same') and variations thereof, but I'm having doubts about the results, so any feedback is greatly appreciated. Do I have to multiply by the time step? Even then, however, the result seems wrong to me (too small). I have to "weigh" the tropospheric curve with the age spectrum. So to state it clearly, my question is: How do you convolve the tropospheric curve with the spectrum G, as prescribed by the first equation?

For those interested, the method is described in Stiller et al. 2012 under "non-linearity correction".

The result I get when I do a convolution of the age spectrum and the tropospheric reference curve looks something like this:

This is just for one Γ. You can do it for any Γ you like. According to the method, the modeled SF6 (i.e. the convolution) at measurement time (2009) should be larger than the original SF6 (red dashed line), but in my case it is always smaller, unless I drop the *dt in the code, then it is always way too large. If you calculated the new SF6 for several different values of Γ, you apply the lookup method and pick the Γ that causes a value of SF6 that is closest to your measurement, to get the actual mean age (and not just the lag time). So if I understand the method correctly, the blue curve should come above the red dashed line to get an SF6' > SF6.

Here is the code:

import numpy as np
import xarray as xr
import gdf_tools as GDF
import matplotlib.pyplot as plt
import pandas as pd
import sys
from datetime import datetime as dt
import time
from scipy import optimize
import scipy as sp
from scipy import interpolate
import datetime
from scipy.signal import savgol_filter
from scipy import stats
from scipy import signal

#---------------------------------------------------------------------------
#                             FUNCTIONS
#---------------------------------------------------------------------------

def smooth(y):
    return signal.savgol_filter(y, 21, 2)  
    
    
def func(x, a, b):
    return a*x + b


def spec(t, *args):
    gamma, d = args
    tp = t / gamma
    return ( 1 / ( 2 * d * np.sqrt( np.pi * tp**3 ) ) ) * \
        np.exp( ( -(gamma**2) * ((tp - 1)**2 )) /( 4 * (d**2) * tp ) )

def getYfromX(y, xval, x):
    return np.interp(xval, x, y)
    
def getXfromY(y, yval, x):
    return np.interp(yval, y, x)
    

#---------------------------------------------------------------------------
#                               MAIN
#---------------------------------------------------------------------------
# Read tropospheric reference curve
ref=pd.read_csv('~/pyscripts1/aoa/refcurve/refcurve.csv')
ref.index=ref["time"].astype("datetime64[D]")

t = ref['decimal_time'].values # Decimal time
t_trunc = t[t<2009] # Truncate at measurement time
z = ref["GML_Global_SF6"].values
ref_trunc =  z[:len(t_trunc)] # Truncate at measurement time

# Linear approximation of the reference curve
a,b = np.polyfit(t_trunc, ref_trunc, 1)  
fit = func(t_trunc, a, b)

# Smoothed reference curve
y=smooth(ref_trunc)

# Age of air (lag time) sequence to loop over
ages = np.linspace(3.5, 8.5, 200)

# Modeled (i.e. "propagatec") SF6 vmr
concentrations = []

# The measured vmr in 2009 is approx 5.25e-12
# Calculate a first guess age of air by projecting this measurement on
# the tropospheric reference and calculating the lag time
proj=getXfromY(y, 5.25, t_trunc)
print("First guess: ", 2009-proj)

# Loop over lag times/ages
for i in ages:      
    gamma = i   # Mean age
    delta = np.sqrt(0.7*gamma)  # Spread of the age spectrum
    deltat = 2009 - t_trunc   # Lag times
    dt = t[1] - t[0]    # Time increment
    
    # Calculate age spectrum
    kernel = spec(deltat, gamma, delta)
    kernel[np.isnan(kernel)] = 0
    
    # Normalise age spectrum
    I = np.trapz(kernel, t_trunc)
    kernel = kernel / I

    # Calculate "weighted average" of the tropo. ref.
    model=np.convolve(ref_trunc, kernel, mode="same")*dt
    
    concentrations.append(model[-1])


# Plot kernel
plt.figure()
plt.plot(t_trunc, kernel)

# Plot tropospheric reference curve and the model
plt.figure()
plt.plot(t_trunc, model, label="convolution")
#plt.plot(t_trunc, ref_trunc, label="trop")
#plt.plot(t_trunc, fit, label="linear fit")
plt.plot(t_trunc, y, label="smooth trop")
plt.axhline(y=5.25, ls='--', c='red')
plt.axvline(x=proj, ls='--', c='black')
plt.legend()
plt.grid()
plt.xlabel("Time")
plt.ylabel("SF6 [vmr]")
plt.tight_layout()

# Plot modeled SF6 as a function of age
# This is the "look-up" method described in Fritsch et al. 2020
plt.figure()
plt.plot(ages, concentrations)
plt.axhline(y=5.25, ls='--', c='red')
plt.xlabel("Mean age of air")
plt.ylabel("SF6 [vmr]")
plt.grid()
plt.tight_layout()


plt.show()