update before resubmission

src/01-model/01-initial_conditions.py

@@ -3,7 +3,7 @@
 import xarray as xr
 from xhistogram.xarray import histogram
 from scipy.optimize import curve_fit
-
+from cmcrameri import cm
 
 chla = xr.open_dataset("../../data/interim/bioargo_north_atlantic.nc")
 chla = chla.assign(biomass=chla.CPHL_ADJUSTED*(16/(0.06*106*14)))
@@ -56,24 +56,64 @@ def func_shortwave(t, a, b, c, d):
 titles["chla"] = f"{popt_chla[0]}(tanh({popt_chla[1]}(z+{-popt_chla[2]}))+1)/2 "
 
 
+
+
+x = xr.DataArray(np.arange(-50,50+1), dims="x")*1e3
+y = xr.DataArray(np.arange(-250,250+1), dims="y")*1e3
+z = xr.DataArray(np.arange(-1000,0+1), dims="z")
+
+L = (99)*1e3/10
+amp = 1e3
+g = 9.82
+ρₒ = 1026
+
+# background density profile based on Argo data
+bg = lambda z: -0.147 * np.tanh( 2.6 * ( z + 623 ) / 1000 ) - 1027.6
+
+# decay function for fronts
+decay = lambda z: ( np.tanh( (z + 500) / 300) + 1 ) / 2
+
+# front function
+front = lambda x, y, z, cy: np.tanh( ( y - ( cy + np.sin(np.pi * x / L) * amp ) ) / 12e3 )
+
+
+D = lambda x, y, z: bg(z) + 0.8*decay(z)*((front(x, y, z, -100e3)+front(x, y, z, 0)+front(x, y, z, 100e3))-3)/6
+
+B = lambda x, y, z: -(g/ρₒ)*D(x, y, z)
+
+bi = B(x,y,z).assign_coords(x=x, y=y, z=z).mean("x")
+
+
+
+
+
 def make_figure():
     colors = {
         "obs": "#1E88E5",
         "model": "#D81B60",
     }
 
-    fig = plt.figure(figsize=(6,6), constrained_layout=True)
+    fig = plt.figure(figsize=(8,6), constrained_layout=True)
     ax = fig.subplot_mosaic(
         [
-            ["pden", "chla"],
-            ["no3", "shortwave"],
+            ["pden", "b", "b"],
+            ["shortwave", "chla", "no3"],
         ],
     )
 
     _ = (argo.PDEN-1000)[::100].plot.line(ax=ax["pden"], y="z", ylim=[1000,0], lw=0, marker="o", color=colors["obs"])
     _ = (argo.PDENf-1000).plot.line(ax=ax["pden"], y="z", ylim=[1000,0], color=colors["model"])
     ax["pden"].text(0.5, 0.02, "Argo", ha="center", color="0.3", fontsize=12, transform=ax["pden"].transAxes)
 
+    bi.plot.contour(ax=ax["b"], levels=15, colors="0.2")
+    bi.differentiate("y").plot.contourf(
+                    ax=ax["b"],
+                    vmin=-1e-7,
+                    cmap=cm.acton,
+                    levels=(-1e-7)*np.arange(0,1.2,0.2),
+                    cbar_kwargs=dict(label="$\partial_y$b [s$^{-2}$]"),
+        )    
+
     cdf_chla.plot.contourf(ax=ax["chla"], levels=[0.2,0.8], colors=colors["obs"], alpha=0.5, extend="neither", add_colorbar=False)
     cdf_chla.plot.contour(ax=ax["chla"], levels=[0.5], colors=colors["obs"])
     ax["chla"].text(0.5, 0.02, "BioArgo", ha="center", color="0.3", fontsize=12, transform=ax["chla"].transAxes)
@@ -100,18 +140,28 @@ def make_figure():
         xlim=[27.4,27.75],
     )
 
+    ax["b"].set(
+        xlabel="y [km]",
+        ylabel="z [m]",
+        yticks=-np.arange(0,1000,200),
+        ylim=[-1000,-10],
+        xticks = np.arange(-200,200+1,100)*1e3,
+        xticklabels = np.arange(-200,200+1,100),
+    )    
+
     ax["chla"].set(
         ylabel="z [m]",
         yticks=-np.arange(0,1000,200),
         ylim=[-1000,-10],
-        xlabel="Chl-a [mmol N m$^{-3}$]",
+        xlabel="Phytoplankton ($\mathcal{P}\,$)\n[mmol N m$^{-3}$]",
     )
 
     ax["no3"].set(
         ylabel="z [m]",
         yticks=-np.arange(0,1000,200),
         ylim=[-800,-10],
-        xlabel="Nitrate [mmol N m$^{-3}$]"
+        xlim=[10,20],
+        xlabel="Nitrate ($\mathcal{N}_n\,$)\n[mmol N m$^{-3}$]"
     )
 
 
@@ -121,18 +171,26 @@ def make_figure():
         title="",
         ylim=[20,280],
         xlim=[0,370],
-        xticks=np.arange(0,365,50),
+        xticks=np.arange(0,365,100),
     )
     # title=f"{a} sin( 2$\pi$ ( t + t$_0$ ) / {b} + {c} ) + {d}"
 
     _ = [ax[k].grid(True, linestyle="--", alpha=0.5) for k in ax]
-    letters = "a b c d".split()
-    _ = [ax[k].set(title=f"{letter})"+50*" ") for k,letter in zip(ax, letters)]
+    letters = "a b c d e".split()
+    for k,letter in zip(ax, letters):
+        whitespace = 70 if k=="b" else 30
+        ax[k].set(title=f"{letter})"+whitespace*" ")
 
     # kw = dict(fontsize=9, color=colors["model"], rotation=90, va="center")
     # _ = [ax[k].text(0.93,0.5,titles[k], **kw, transform=ax[k].transAxes) for k in titles]
 
     return fig,ax
 
 fig,ax = make_figure()
-fig.savefig("../../reports/figures/data_driven_initial_conditions.png", facecolor="w", dpi=300, bbox_inches="tight")
+fig.savefig("../../reports/figures/data_driven_initial_conditions.png", facecolor="w", dpi=300, bbox_inches="tight")
+
+
+
+
+
+

src/01-model/02-submesoscale.jl

@@ -163,7 +163,7 @@ const ρₒ = 1026
 @inline B(x, y, z) = -(g/ρₒ)*D(x, y, z)
 
 # initial phytoplankton profile
-@inline P(x, y, z) = 0.1 * ( tanh( 0.01 * (z + 300)) + 1) / 2
+@inline P(x, y, z) = 0.2 * ( tanh( 0.01 * (z + 300)) + 1) / 2
 
 # initial nitrate profile
 @inline N(x, y, z) = z * (12 - 16) / (0 + 800) + 12

src/02-analysis/03-Ro_restratification.py

@@ -111,7 +111,7 @@
         C1 = (
             (ps[p]["factor"]*data[k]["D"][p].mean(["xC"])).T
         ).plot.contourf(ax=a, **ps[p]["kw"])
-        a.set(xlim=[1,80], ylabel="y [km]", xlabel="time [days]")
+        a.set(xlim=[1,40], ylabel="y [km]", xlabel="time [days]")
         text = a.text(0.97, 0.03, data[k]["label"], 
                transform=a.transAxes, **kw_text)
         text.set_path_effects([path_effects.Stroke(linewidth=2, foreground='w'),

src/02-analysis/04-biogeochemistry.py

@@ -28,7 +28,7 @@
 mus = [0.5,0.75,1.0,1.25]
 
 
-for k in data:
+for k in tqdm(data):
     data[k]["D"] = dict()
     for mu in mus:
         ds = xr.open_dataset(f"../../data/raw/output_{k}{sinking_text}_mu{mu}.nc").isel(time=slice(None,None,3))
@@ -109,6 +109,8 @@
     integ.append(integi.expand_dims("mu").assign_coords(mu=[mu]))
 integ = xr.concat(integ,"mu")
 
+print(integ.sel(time=slice(12,32)).integrate("time")-integ.sel(time=slice(12,32),exp="submesoscale").integrate("time"))
+print(1-integ.sel(time=slice(12,32)).integrate("time")/integ.sel(time=slice(12,32)).integrate("time").sel(exp="submesoscale"))
 
 colors = ["0.1","0.4","0.7","0.8"]
 titles = [
@@ -129,172 +131,3 @@
 _ = [a.set(ylabel="") for a in ax[1:]]
 fig.savefig(f"../../reports/figures/integrated_new_production{sinking_text}.png", facecolor="w", dpi=300, bbox_inches="tight")  
 
-print(1-integ.sel(time=slice(12,32)).integrate("time")/integ.sel(time=slice(12,32)).integrate("time").sel(exp="submesoscale"))
-
-
-
-
-
-# step = 5
-# bins = np.arange(-(100+step/2),100+step/2,step)
-
-# nflux = (submeso.w*submeso.N*86400).sel(zC=slice(-100,0))
-# nflux.name = "N_flux"
-
-# H = histogram(nflux, bins=bins, dim=["xC", "yC"]).T
-
-# fig,ax = plt.subplots()
-# H.sum("zC").plot(robust=True)
-# nflux.median(["xC", "yC", "zC"]).plot(x="time")
-# ax.set(ylim=[-50,50])
-
-# submeso_mld = submeso.interp(zC=-submeso.h).compute()
-# coarse_mld = coarse.interp(zC=-coarse.h).compute()
-
-
-# tis = [15,25,50]
-# vm = 0.04
-# kw = dict(
-#     nflux=dict(vmin=-vm, vmax=vm, add_colorbar=False, cmap="seismic")
-# )
-# fig, ax = plt.subplots(2,len(tis))
-# for col,ti in enumerate(tis):
-#     C = (submeso_mld.w*submeso_mld.N).sel(time=ti).plot(ax=ax[0,col], **kw["nflux"])
-#     (coarse_mld.w*coarse_mld.N).sel(time=ti).plot(ax=ax[1,col], **kw["nflux"])
-# fig.colorbar(C, ax=ax)
-
-
-# fig, ax = plt.subplots(2,len(tis))
-# for col,ti in enumerate(tis):
-#     C = (submeso.w*submeso.N).sel(zC=-100, time=ti, method="nearest").plot(ax=ax[0,col], **kw["nflux"])
-#     (coarse.w*coarse.N).sel(zC=-100, time=ti, method="nearest").plot(ax=ax[1,col], **kw["nflux"])
-# fig.colorbar(C, ax=ax)
-
-
-# submeso.new_production.mean(["xC", "yC"]).integrate("zC").plot()
-# coarse.new_production.mean(["xC", "yC"]).integrate("zC").plot()
-
-
-
-# fig, ax = plt.subplots(1,2)
-# submeso.h.mean("xC").plot(ax=ax[0])
-# coarse.h.mean("xC").plot(ax=ax[1])
-
-# submeso_coarsen = submeso.h.mean("xC").coarsen(yC=10).mean().interp(yC=coarse.yC)
-
-# fig,ax = plt.subplots()
-# (submeso_coarsen-coarse.h.mean("xC")).plot(vmin=-100,vmax=100,cmap="RdBu")
-# vm = 10
-# (submeso_coarsen-coarse.h.mean("xC")).plot.contour(levels=[-vm,vm], colors="0.2")
-# ax.axvline(-50)
-# ax.axvline(170)
-
-# data = dict(
-#     submeso=dict(D=submeso, label="1 km"),
-#     coarse=dict(D=coarse, label="10 km"),
-# )
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-# fig,ax = plt.subplots(3,1)
-
-# for p,a in zip(["P", "N", "new_production"], np.ravel(ax)):
-
-#     dims = ["xC", "yC", "zC"]
-#     average = lambda A: A.mean(dims)
-
-#     for D in [submeso, coarse]:
-#         var = average((86400*D[p]))
-#         # print(var.sel(time=slice(12,80)).integrate("time"))
-#         var.plot(ax=a)
-#     a.set(xlim=[0,80])
-
-
-
-# # z_mld = (submeso.zC/submeso.h).chunk({"time":1})
-# # z_mld.name = "z_mld"
-
-# # nflux = (submeso.w*submeso.N).chunk({"time":1})*86400
-# # nflux.name = "nflux"
-
-# # bins = [
-# #     np.linspace(-2,0,20),
-# #     np.linspace(-1e3,1e3,700),
-# # ]
-# # H = histogram(z_mld, nflux, bins=bins).load()
-# # H = H/H.sum("nflux_bin")
-
-# # # H.sel(nflux_bin=slice(-60,60)).plot(vmax=0.01)
-# # ((H*H.nflux_bin).sum("nflux_bin")/H.sum("nflux_bin")).plot(y="z_mld_bin", color="r")
-# # # H.cumsum("nflux_bin").plot.contour(levels=[0.01,0.2,0.5,0.8,0.99], colors="0.5")
-
-
-# # z_mld = (coarse.zC/coarse.h).chunk({"time":1})
-# # z_mld.name = "z_mld"
-
-# # nflux = (coarse.w*coarse.N).chunk({"time":1})*86400
-# # nflux.name = "nflux"
-
-# # bins = [
-# #     np.linspace(-2,0,20),
-# #     np.linspace(-1e3,1e3,700),
-# # ]
-# # H = histogram(z_mld, nflux, bins=bins).load()
-# # H = H/H.sum("nflux_bin")
-
-# # # H.sel(nflux_bin=slice(-60,60)).plot(vmax=0.01)
-# # ((H*H.nflux_bin).sum("nflux_bin")/H.sum("nflux_bin")).plot(y="z_mld_bin", color="r")
-# # # H.cumsum("nflux_bin").plot.contour(levels=[0.01,0.2,0.5,0.8,0.99], colors="0.5")
-
-
-
-
-# # bins = [
-# #     np.linspace(0,1e-7,50),
-# #     np.linspace(0,0.1,50),
-# # ]
-
-# # zmin = -100
-# # H = histogram(submeso.sel(zC=slice(zmin,0))["∇b"].chunk({"time":1}), 
-# #               86400*submeso.sel(zC=slice(zmin,0)).new_production.chunk({"time":1}), bins=bins).compute()
-# # H.plot(robust=True)
-
-# # H = histogram(coarse.sel(zC=slice(zmin,0))["∇b"].chunk({"time":1}), 
-# #               86400*coarse.sel(zC=slice(zmin,0)).new_production.chunk({"time":1}), bins=bins).compute()
-# # H.plot(robust=True)
-
-# # submeso.new_production.mean(["xC", "yC"]).integrate("zC").plot()
-# # coarse.new_production.mean(["xC", "yC"]).integrate("zC").plot()
-
-# # tslice = slice(12,80)
-
-# # da = 86400*submeso.new_production.mean(["xC", "yC"]).integrate("zC").sel(time=tslice)
-# # da = (da*np.gradient(da.time)).cumsum("time")
-
-# # db = 86400*coarse.new_production.mean(["xC", "yC"]).integrate("zC").sel(time=tslice)
-# # db = (db*np.gradient(db.time)).cumsum("time")
-
-# # ((da-db)/db).plot()
-
-
-# submeso.new_production.mean(["xC", "yC"]).T.sel(zC=slice(-60,0), time=slice(12,80)).plot.contourf(levels=10,vmax=3e-5)
-# coarse.new_production.mean(["xC", "yC"]).T.sel(zC=slice(-60,0), time=slice(12,80)).plot.contourf(levels=10,vmax=3e-5)
-
-
-# submeso.new_production.mean(["xC"]).integrate("zC").T.sel(time=slice(12,80)).plot.contourf(levels=10,vmin=0, vmax=1e-3)
-# coarse.new_production.mean(["xC"]).integrate("zC").T.sel(time=slice(12,80)).plot.contourf(levels=10,vmin=0, vmax=1e-3)
-
-
-# submeso.new_production.mean(["xC", "yC"]).integrate("zC").T.sel(time=slice(12,80)).plot()
-# coarse.new_production.mean(["xC", "yC"]).integrate("zC").T.sel(time=slice(12,80)).plot()

src/02-analysis/05-nitrate_flux.py

@@ -148,18 +148,16 @@
 (
     (
         (86400*datamld.sel(exp="submesoscale").w*datamld.sel(exp="submesoscale").N)
-        .mean(["xC","yC"]).sel(time=slice(12,40),mld=1).integrate("time")
+        .mean(["xC","yC"]).sel(time=slice(12,32),mld=1).integrate("time")
     ),
     (
         (86400*(datamld.sel(exp="coarse_mle").w+datamld.sel(exp="coarse_mle").w_mle)*datamld.sel(exp="coarse_mle").N)
-        .mean(["xC","yC"]).sel(time=slice(12,40),mld=1).integrate("time")
+        .mean(["xC","yC"]).sel(time=slice(12,32),mld=1).integrate("time")
     )
 )
 
 
 
-
-
 kw = lambda vm: dict(vmin=-vm,vmax=vm,levels=np.arange(-vm,vm+20,20),cmap="bwr",add_colorbar=False)
 ti = 15
 fig, ax = plt.subplots(2,2,figsize=(10,7))