notebook.community

Edit and run 





In [ ]:



    


%matplotlib inline
import matplotlib as mpl
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
import numpy as np
from scipy.interpolate import interp1d
from scipy.optimize import brentq
from scipy.signal import tukey
from scipy.stats import gamma
from matplotlib.colors import Normalize, LogNorm
from sklearn.neighbors.kde import KernelDensity
import pandas as pd
import clemb
from clemb.forward_model import forward_model, fullness, esol, rk4, es
from clemb.syn_model import SynModel
import IPython.core.debugger
dbg = IPython.core.debugger.Pdb()



    











In [ ]:



    


def error(y1, y2, t1, t2):
    f1 = interp1d(t1, y1)
    err = y2 - f1(t2)
    return err



    











In [ ]:



    


if True:
    c = clemb.Clemb(clemb.LakeDataFITS(), clemb.WindDataCSV(), start='2018-02-01 ', end='2018-06-06')



    











In [ ]:



    


if True:
    rp = c.run_backward()



    











In [ ]:



    


pwr = rp['pwr']*0.0864
datetime = pd.to_datetime(rp['dates'].data)
melt = rp['fmelt']
orig_steam = rp['steam']
orig_mevap = rp['evfl']
ll = rp['llvl']
t = rp['t']

ndays = (datetime[-1] - datetime[0]).days
dt = 1.
nsteps = int(ndays/float(dt))
time = np.linspace(0,ndays,nsteps+1)
T0 = t.data[0]

y = np.zeros((nsteps+1, 3))
a,vol = fullness(ll.data)
density = 1.003 - 0.00033 * T0
M0 = vol[0]*density
y[0] = [T0, M0, 0]
steams = np.zeros(nsteps+1)
mevaps = np.zeros(nsteps+1)
esols = np.zeros(nsteps+1)
for i in range(nsteps):
    solar = esol(dt, a[i], datetime[i].month)
    esols[i] = solar
    Mout = 0.
    enthalpy = 6.0
    windspeed = 4.5
    y_new, steam, mevap = forward_model(y[i], dt, a[i], vol[i], pwr.values[i], 
                                 melt.values[i], Mout, solar, enthalpy, windspeed)
    steams[i] = steam
    mevaps[i] = mevap
    y[i+1] = y_new
    
mpl.rcParams['figure.subplot.hspace'] = 0.5
plt.figure(figsize=(10,8))
ax1 = plt.subplot(4,2,1)
ax1.plot(time, pwr.values[:]/0.0864, ls='--')
ax1.set_title('Heat input rate')

ax2 = plt.subplot(4,2,3)
ax2.plot(time, melt.values[:], ls='--')
ax2.set_title('Meltwater inflow')

ax3 = plt.subplot(4,2,5)
ax3.plot(time, mevaps)
ax3.plot(time, orig_mevap, ls='--')
ax3.set_title('Evaporation')

ax4 = plt.subplot(4,2,2)
ax4.plot(time, y[:,0])
ax4.plot(time, t.values[:], ls='--')
ax4.set_title('Lake temperature')

ax5 = plt.subplot(4,2,4)
ax5.plot(time, y[:,0] - t.values[:])
ax5.set_title('Temperature difference')

ax6 = plt.subplot(4,2,6)
ax6.plot(time, steams)
ax6.plot(time, orig_steam, ls='--')
ax6.set_title('Steam')

ax7 = plt.subplot(4,2,7)
ax7.plot(time, y[:,1])
_, vol = fullness(ll.values[:])
rho = 1.003 - 0.00033 * t.values[:]
mass = vol*rho
ax7.plot(time, mass, ls='--')
ax7.set_title('Lake mass')

ax8 = plt.subplot(4,2,8)
ax8.plot(time, y[:,1]- mass)
ax8.set_title('Mass difference')

Synthetic input





In [ ]:



    


# matplotlib.rcParams['axes.prop_cycle']
cl1 = '#1f77b4'
cl2 = '#ff7f0e'
cl3 = '#2ca02c'



    











In [ ]:



    


def plot_syn_model(df, tm):
    t = df.index
    
    mpl.rcParams['figure.subplot.hspace'] = 0.5
    fig, axs = plt.subplots(nrows=5, ncols=2, figsize=(14, 10))
    
    axs[0,0].plot(t, tm[:,5], ls='--')
    axs[0,0].set_title('Heat input rate')

    axs[0,1].plot(t, tm[:, 0])
    axs[0,1].set_title('Meltwater inflow')

    axs[1,0].plot(t, df['M'] + df['M_err'], 'k+')
    axs[1,0].plot(t, df['M'])
    axs[1,0].set_title('Lake mass')

    axs[1,1].plot(t, df['T'] + df['T_err'],'k+')
    axs[1,1].plot(t, df['T'])
    axs[1,1].set_title('Lake temperature')
    
    axs[2,0].plot(t, tm[:,1])
    axs[2,0].set_title('Outflow')
    
    axs[2,1].plot(t, df['z'] + df['z_err'],'k+')
    axs[2,1].plot(t, df['z'])
    axs[2,1].set_title('Lake level')
    
    axs[3,0].plot(t, df['X'] + df['X_err'], 'k+')
    axs[3,0].plot(t, df['X'])
    axs[3,0].set_title('Mg++ total amount')

    axs[3,1].plot(t, tm[:, 3])
    axs[3,1].set_title('Evaporated steam')

    axs[4,0].plot(t, df['v'] + df['v_err'], 'k+')
    axs[4,0].plot(t, df['v'])
    axs[4,0].set_title('Lake volume')
    
    axs[4,1].plot(t, df['Mg'] + df['Mg_err'], 'k+')
    axs[4,1].plot(t, df['Mg'])
    axs[4,1].set_title('Mg++ concentration')



    











In [ ]:



    


c = clemb.Clemb(None, None, None, None, pre_txt='syn1')



    











In [ ]:



    


df, prm = SynModel().run(100, 1000.)
c._df = df
c._dates = df.index
plot_syn_model(df, prm)



    











In [ ]:



    


c.drmg = True
c.fullness = SynModel(area=df['a'].values).synth_fullness
rsb = c.run_backward(new=True)



    











In [ ]:



    


rs = c.run_forward(nsamples=2000, nresample=1000, m_out_max=40., m_in_max=40., q_in_max=1500., new=True)



    











In [ ]:



    


def density_plot(ax, rs, prm, prm_lim=(0, 1800, 100),  mode='kde'):
    nsteps = rs.dims['dates']
    if mode == 'kde':
        bins = np.arange(*prm_lim)
        m = []
        for i in range(nsteps):
            y = rs[prm][i].data
            idx = np.isnan(y)
            _kde = KernelDensity(kernel='gaussian', bandwidth=60.).fit(y[~idx].reshape(-1,1))
            X_plot = np.linspace(0, 1600, 1000)
            log_dens = _kde.score_samples(X_plot[:, np.newaxis])
            bins=X_plot
            m.append(np.exp(log_dens))
        m = np.array(m)
        ax.contourf(np.arange(nsteps), bins,  m.T, 30, cmap=plt.cm.get_cmap('RdBu_r'))
    if mode == 'scatter':
        datet = rs['dates']
        nresample = rs[prm].shape[1]
        for k in range(nsteps):
            ax.scatter([datet[k].data]*nresample, rs[prm][k], s=2, c=rs['lh'][k],
                    cmap=plt.cm.get_cmap('RdBu_r'), alpha=0.3)
    return



    











In [ ]:



    


import matplotlib.colors as colors

def plot_inv_res(irs, df, tm, irsb=None):
    days = mdates.DayLocator()  # every day
    months = mdates.MonthLocator()
    monthFmt = mdates.DateFormatter('%Y-%m-%d')
    dayFmt = mdates.DateFormatter('%d')

    datet = irs['dates']
    data = df
    t_data = data.index
    nsteps = datet.size
    nresample = irs['q_in'].shape[1]

    show_exp = True

    mpl.rcParams['figure.subplot.hspace'] = 0.5
    fig, axs = plt.subplots(nrows=5, ncols=2, figsize=(14, 10))

    axs[0,0].plot(t_data, np.ones(t_data.size)*4.5, ls='--')
    axs[0,0].set_title('Wind speed [m/s]')

    axs[0,1].plot(t_data, data['X'], ls='--')

    model_X = irs['model'].sel(obs='X')
    axs[0,1].fill_between(t_data, data['X']-3*data['X_err'],
                          data['X']+3*data['X_err'], alpha=0.5)
    if show_exp:
        axs[0,1].plot(datet, model_X.mean(axis=1), 'k-')
        axs[0,1].plot(datet, model_X.mean(axis=1)+3*model_X.std(axis=1), 'k--')
        axs[0,1].plot(datet, model_X.mean(axis=1)-3*model_X.std(axis=1), 'k--')
    axs[0,1].set_title('Mg++ amount [kt]')

    axs[1,0].plot(t_data, data['T'], ls='--')
    axs[1,0].fill_between(t_data, data['T']-3*data['T_err'],
                          data['T']+3*data['T_err'], alpha=0.5)

    model_T = irs['model'].sel(obs='T')
    if show_exp:
        axs[1,0].plot(datet, model_T.mean(axis=1), 'k-')
        axs[1,0].plot(datet, model_T.mean(axis=1)+3*model_T.std(axis=1), 'k--')
        axs[1,0].plot(datet, model_T.mean(axis=1)-3*model_T.std(axis=1), 'k--')
    axs[1,0].set_title('Lake temperature [$^{\circ}C$]')

    axs[1,1].plot(t_data, data['M'], ls='--')
    axs[1,1].fill_between(t_data, data['M']-3*data['M_err'],
                          data['M']+3*data['M_err'], alpha=0.5)
    model_M = irs['model'].sel(obs='M')
    if show_exp:
        axs[1,1].plot(datet, model_M.mean(axis=1), 'k-')
        axs[1,1].plot(datet, model_M.mean(axis=1)+3*model_M.std(axis=1), 'k--')
        axs[1,1].plot(datet, model_M.mean(axis=1)-3*model_M.std(axis=1), 'k--')
    axs[1,1].set_title('Lake mass [kt]')

    density_plot(axs[2,0], irs, 'q_in', mode='scatter')
    exp_q_in = irs['exp'].loc[:,'q_in']
    var_q_in = irs['var'].loc[:,'q_in']
    idx = np.argmax(irs['lh'], axis=1)
    map_q_in = irs['q_in'][:, idx]
    if show_exp:
        axs[2,0].plot(datet, exp_q_in, 'k')
        axs[2,0].plot(datet, exp_q_in - 3*np.sqrt(var_q_in), 'k--')
        axs[2,0].plot(datet, exp_q_in + 3*np.sqrt(var_q_in), 'k--')
    axs[2,0].plot(datet, tm[:-1,5], 'b')
    axs[2,0].plot(datet, map_q_in, 'g')
    if irsb is not None:
        axs[2,0].plot(irsb['dates'], irsb['pwr'], 'g-')
    axs[2,0].set_title('Heat input rate [MW]')

    density_plot(axs[2,1], irs, 'h', mode='scatter')
    exp_h = irs['exp'].loc[:,'h']
    axs[2,1].plot(datet, exp_h, 'k')
    axs[2,1].set_title('Enthalpy [TJ/kt]')

    density_plot(axs[3,0], irs, 'm_in', mode='scatter')
    exp_m_in = irs['exp'].loc[:,'m_in']
    var_m_in = irs['var'].loc[:,'m_in']
    if show_exp:
        axs[3,0].plot(datet, exp_m_in, 'k')
        axs[3,0].plot(datet, exp_m_in - 3*np.sqrt(var_m_in), 'k--')
        axs[3,0].plot(datet, exp_m_in + 3*np.sqrt(var_m_in), 'k--')
    axs[3,0].plot(datet, tm[:-1,0], 'b')
    axs[3,0].set_title('Inflow [kt/day]')

    density_plot(axs[3,1], irs, 'm_out', mode='scatter')
    exp_m_out = irs['exp'].loc[:,'m_out']
    var_m_out = irs['var'].loc[:,'m_out']
    if show_exp:
        axs[3,1].plot(datet, exp_m_out, 'k')
        axs[3,1].plot(datet, exp_m_out - 3*np.sqrt(var_m_out), 'k--')
        axs[3,1].plot(datet, exp_m_out + 3*np.sqrt(var_m_out), 'k--')
    axs[3,1].plot(datet, tm[:-1,1], 'b')
    axs[3,1].set_title('Outflow [kt/day]')

    h = 6.0
    density_plot(axs[4,0], irs, 'steam', mode='scatter')
    axs[4,0].plot(datet, tm[:-1,5]*0.0864/h, 'b')
    axs[4,0].set_title('Steam input [kt/day]')

    density_plot(axs[4,1], irs, 'mevap', mode='scatter')
    _ = axs[4,1].set_title('Evaporation mass loss [kt/day]')

    for row in range(5):
        for col in range(2):
            axs[row,col].xaxis.set_major_locator(months)
            axs[row,col].xaxis.set_major_formatter(monthFmt)
            axs[row,col].xaxis.set_minor_locator(days)
            axs[row,col].set_xlim(datet[0].data, datet[-1].data)



    











In [ ]:



    


plot_inv_res(rs, c._df, prm)



    











In [ ]:



    


exp_q_in = rs['exp'].loc[:,'q_in']
fig = plt.figure(figsize=(16,6))
ax = fig.add_subplot(121)
density_plot(ax, rs, 'q_in', mode='kde')
ax.plot(prm[:,5], 'k-')
ymin, ymax = ax.get_ylim()
ax.vlines(20, ymin, ymax, linestyle='--', color='white')
new_labels = []
new_ticks = []
dates = rs['dates'].data
for _xt in plt.xticks()[0]:
    try:
        dt = dates[int(_xt)].astype('datetime64[us]').min()
        new_labels.append((pd.to_datetime(str(dt))
                            .strftime("%Y-%m-%d")))
        new_ticks.append(_xt)
    except IndexError:
        continue
_, _ = plt.xticks(new_ticks, new_labels, rotation=30,
           horizontalalignment='right')

ax.set_ylabel('Energy input rate [MW]')
ax.set_xlim(0, nsteps-1)

ax1 = fig.add_subplot(122)
x = rs['q_in'][20].data
idx = np.isnan(x)
kde = KernelDensity(kernel='gaussian', bandwidth=60.).fit(x[~idx].reshape(-1,1))
X_plot = np.linspace(0, 1600, 1000)[:, np.newaxis]
log_dens = kde.score_samples(X_plot)
y = np.exp(log_dens)
ax1.plot(X_plot[:, 0], y)
#plt.plot(bins[:-1], m[20]/m[20].max())
ax1.set_xlabel('Energy input rate [MW]')
fout = '/media/win_home/Geosciences_2018/syn_example.png'
#plt.savefig(fout, dpi=300, bbox_inches='tight')



    











In [ ]:



    


x = rs['q_in'][10].data
idx = np.isnan(x)
kde = KernelDensity(kernel='gaussian', bandwidth=20.).fit(x[~idx].reshape(-1,1))
X_plot = np.linspace(0, 1600, 1000)[:, np.newaxis]
log_dens = kde.score_samples(X_plot)
y = np.exp(log_dens)
plt.plot(X_plot[:, 0], y)



    











In [ ]:



    


from pandas.plotting import scatter_matrix

sdf = pd.DataFrame({'q_in': rs['q_in'][20].data, 'm_in': rs['m_in'][20].data,
                    'm_out': rs['m_out'][20].data})
scatter_matrix(sdf, alpha=0.5, figsize=(6, 6), diagonal='kde')
fout = '/media/win_home/Geosciences_2018/scatter_matrix.png'
#plt.savefig(fout, dpi=300, bbox_inches='tight')



    











In [ ]:



    


fig, axs = plt.subplots(nrows=1, ncols=2, figsize=(16,6))
for q in [100, 500, 1000, 1500]:
    df, prm = SynModel().run(100, q)
    T = df['T']
    qi = prm[:, 5]
    d = df.index
    #T, qi, d = get_T(q, mode='gamma')
    axs[1].plot(d, T)
    axs[0].plot(d, qi)
idx = np.argmax(qi)
ymin, ymax = axs[0].get_ylim()
axs[0].vlines(d[idx], ymin, ymax, linestyle='--')
axs[0].set_ylim(0,ymax)
axs[0].set_ylabel('Heat input rate [MW]')
ymin, ymax = axs[1].get_ylim()
axs[1].vlines(d[idx], ymin, ymax, linestyle='--')
axs[1].set_ylim(ymin,ymax)
axs[1].set_ylabel('Temperature $^{\circ}C$')
fout = '/media/win_home/Geosciences_2018/heatpulse.png'
#fig.savefig(fout, dpi=300, bbox_inches='tight')



    











In [ ]:



    


c1 = clemb.Clemb(None, None, None, None, pre_txt='syn2')
df1, prm1 = SynModel().run(31, 1200.)
c1._df = df1
c1._dates = df1.index
rs1 = c1.run_forward(nsamples=1000, q_in_max=1400., new=False)



    











In [ ]:



    


plot_inv_res(rs1, c1._df, prm1)



    











In [ ]:



    


fig = plt.figure()
ax = fig.add_subplot(111)
datet = rs1['dates']
nresample=500
k = 4
T_model = rs1['model'].sel(obs='T')
ax.scatter([datet[k].data]*nresample, np.random.normal(df1['T'][k], df1['T_err'][k], nresample), s=10, c=rs1['lh'][k],
           cmap=plt.cm.get_cmap('RdBu_r'), alpha=0.3,
           norm=colors.Normalize(vmin=rs1['lh'].min(), vmax=rs1['lh'].max()))
k = 5
ax.scatter([datet[k].data]*nresample, T_model[k], s=10, c=rs1['lh'][k],
           cmap=plt.cm.get_cmap('RdBu_r'), alpha=0.3,
           norm=colors.Normalize(vmin=rs1['lh'].min(), vmax=rs1['lh'].max()))
ax.plot(rs1['dates'],df1['T'][1:])
ax.set_xlim(rs1['dates'][3].data, rs1['dates'][7].data)
ax.set_ylim(13.5,20)
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.get_yaxis().set_ticks([])
days = mdates.DayLocator()  # every day
months = mdates.MonthLocator()
monthFmt = mdates.DateFormatter('%Y-%m-%d')
dayFmt = mdates.DateFormatter('%d')
ax.xaxis.set_major_locator(months)
ax.xaxis.set_major_formatter(monthFmt)
ax.xaxis.set_minor_locator(days)
fout = '/media/win_home/Geosciences_2018/mcmc_example.png'
#fig.savefig(fout, dpi=300, bbox_inches='tight')

Equations from Hurst et al., 1991

Energy balance for the lake:

$$ \frac{d}{dt}Q = Q_i - Q_e - cTM_o + cT_sM_s $$

$Q$ = total energy of the lake

$T$ = temperature at any time $t$

$c$ = specific heat of lake water

$Q_i$ = heat input at the base of the lake

$Q_e$ = heat loss due to evaporation, radiation, solar heating (gain)

$M_o$ = total rate of outflow

$M_s$ = inflow rate from melt

$T_s$ = temperature of inflow

Assuming $T_s = 0 ^{\circ}C$:

$$ \begin{aligned} \frac{d}{dt}Q & = Q_i - Q_e - cTM_o \\ & = \frac{d}{dt}[cMT] \\ & = cM\frac{dT}{dt} + cT\frac{dM}{dt} \end{aligned} $$$$ \Rightarrow \qquad \frac{dT}{dt} = \frac{1}{cM}\left(Q_i - Q_e - cTM_o -cT\frac{dM}{dt}\right) $$

$M$ = mass of the water in the lake at time t

Mass balance:

$$ \frac{dM}{dt} = M_i + M_s - M_o - M_e $$

$M_i$ = rate at which water or steam is added through the volcanic vent

$M_e$ = rate of evaporation losses

Ion concentration balance:

$$ \begin{aligned} \frac{d}{dt}[MZ] & = Z_i M_i - Z M_o\\ \frac{dM}{dt} + M_o & = \frac{1}{Z} ( Z_i M_i - M \frac{dZ}{dt} ) \end{aligned} $$

$Z$ = ion concentration in the lake

$Z_iM_i$ = rate of addition of ions through the lake bottom

If $Z_iM_i = 0$:

$$ \frac{dM}{dt} = -M_o - \frac{M}{Z}\frac{dZ}{dt} $$

Equations from code

Lake mass:

$$ \begin{aligned} \rho & = 1.003 - 0.00033T \\ V & = f(h)\\ M & = V * \rho \end{aligned} $$

$\rho$ = density of lake water

$V$ = lake volume

$h$ = lake level

Outflow:

$$ M_o = M_n - M_{n-1} + M_{n-1} \left (0.98\frac{M_{n-1}Z_{n-1}}{M_nZ_n} -1.0 \right ) $$

Long wave radiation loss:

$$ \dot{E}_{long} = 5.67e^{-8}A[0.97(T - 1 + 273.15)^4 - 0.8(0.9 + 273.15)^4] $$

$A$ = lake surface area

Free and forced convection:

$$ \begin{aligned} \dot{E}_{evap} & = \\ & A \sqrt{\left [ 2.2 (T - 1 - 0.9)^{\frac{1}{3}} (6.112e^{\frac{17.62(T-1)}{243.12(T-1)}} - 6.5) \right ] ^2+ \left [ \left ( 4.07e^{-3} \frac{w^{0.8}}{500^{0.2}} - \frac{1.107e^{-2}}{500} \right ) \left ( \frac{1}{800} ( 6.112e^{\frac{17.62(T-1)}{243.12(T-1)}} - 6.5) \right ) 2400000 \right ]^2} \end{aligned} $$

$w$ = wind speed

Solar heating:

$$ \dot{E}_{solar} = \Delta t A 1.5e^{-5} \left [ 1 + 0.5 cos \left (\frac{((m-1)3.14)}{6.0} \right ) \right ] $$

$m$ = month of the year

Heat loss due to evaporation radiation and solar heating

$$ Q_e = \dot{E}_{long} + \dot{E}_{evap} \left ( 1+0.948 \frac{1005}{2.4e^6} \frac{T - 1 - 0.9}{6.335e^{-3} + 6.718e^{-4}(T-1) -2.0887e^{-5}(T-1)^2 + 7.3095e^{-7}(T-1)^3 - 2.2e^{-3}} \right ) - \dot{E}_{solar} $$$$ M_e = \dot{E}_{evap}2.4e^{-12} $$ 




In [ ]:

Content source: yannikbehr/clemb

Similar notebooks:

notebook.community | gallery | about