Finalyse isostasy calibration notebooks. Cleaning ETOPO1 datasets.

tth030 · Jan 9, 2022 · 8c48978 · 8c48978
1 parent f9ce50c
commit 8c48978
Show file tree

Hide file tree

Showing 22 changed files with 3,439 additions and 142 deletions.
diff --git a/...ntralvalley_GL-Tp1280_speed3-perplex.rhoz → ...alvalley_GL-Tp1280_speed3-perplex_M1.rhoz b/...ntralvalley_GL-Tp1280_speed3-perplex.rhoz → ...alvalley_GL-Tp1280_speed3-perplex_M1.rhoz
diff --git a/...ley_GL-Tp1280_speed3-perplex_corrected.tz → ...tralvalley_GL-Tp1280_speed3-perplex_M1.tz b/...ley_GL-Tp1280_speed3-perplex_corrected.tz → ...tralvalley_GL-Tp1280_speed3-perplex_M1.tz
diff --git a/...ridge_centralvalley_GL-Tp1280_speed3.rhoz → ...ge_centralvalley_GL-Tp1280_speed3_M2.rhoz b/...ridge_centralvalley_GL-Tp1280_speed3.rhoz → ...ge_centralvalley_GL-Tp1280_speed3_M2.rhoz
diff --git a/...ntralvalley_GL-Tp1280_speed3_corrected.tz → ...idge_centralvalley_GL-Tp1280_speed3_M2.tz b/...ntralvalley_GL-Tp1280_speed3_corrected.tz → ...idge_centralvalley_GL-Tp1280_speed3_M2.tz
diff --git a/data/isostasy_calibration/rhoc_lithdepletion_3300_bl400_waterload.txt b/data/isostasy_calibration/rhoc_lithdepletion_3300_bl400_waterload.txt
diff --git a/data/isostasy_calibration/rhoc_lithdepletion_3300_waterdepth2750_waterload.txt b/data/isostasy_calibration/rhoc_lithdepletion_3300_waterdepth2750_waterload.txt
diff --git a/data/isostasy_calibration/rhoc_lithdepletion_3311_bl400_waterload.txt b/data/isostasy_calibration/rhoc_lithdepletion_3311_bl400_waterload.txt
diff --git a/data/isostasy_calibration/rhoc_lithdepletion_3311_waterdepth2750_waterload.txt b/data/isostasy_calibration/rhoc_lithdepletion_3311_waterdepth2750_waterload.txt
diff --git a/data/isostasy_calibration/rhoc_lithdepletion_SP2008_speed3_0.55_bl400_waterload.txt b/data/isostasy_calibration/rhoc_lithdepletion_SP2008_speed3_0.55_bl400_waterload.txt
diff --git a/data/isostasy_calibration/rhoc_lithdepletion_SP2008_speed3_0.55_waterdepth2750_waterload.txt b/data/isostasy_calibration/rhoc_lithdepletion_SP2008_speed3_0.55_waterdepth2750_waterload.txt
diff --git a/data/topography/.gitignore b/data/topography/.gitignore
@@ -0,0 +1,6 @@
+ETOPO1_Bed_g_gmt4_1.grd
+ETOPO1_Bed_g_gmt4_2.grd
+ETOPO1_Bed_g_gmt4_3.grd
+ETOPO1_Bed_g_gmt4_4.grd
+ETOPO1_Bed_g_gmt4_5.grd
+ETOPO1_Bed_g_gmt4_6.grd
diff --git a/data/topography/rate.3.2.nc b/data/topography/rate.3.2.nc
diff --git a/isostasy_calib.ipynb b/isostasy_calib.ipynb
diff --git a/isostasy_calib2.ipynb b/isostasy_calib2.ipynb
diff --git a/isostasy_ref.ipynb b/isostasy_ref.ipynb
@@ -9,7 +9,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 1,
+   "execution_count": 5,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -19,7 +19,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 2,
+   "execution_count": 22,
    "metadata": {},
    "outputs": [
     {
@@ -39,42 +39,39 @@
       "Layer 1 T top =  295.3 Q top = 37.68 (matMC)\n",
       "Layer 0 T top =    0.0 Q top = 51.28 (matUC)\n",
       " Qsurface = 51.28 mW/m2\n",
-      "Pref = 3.79540e+09 Pa with compensation depth of 125.0 km\n",
+      "Pref = 3.91001e+09 Pa with compensation depth of 125.0 km\n",
       "====> Average crustal density lith125 is 2833.0813751716532 kg/m3\n",
       "\n",
       "# Column 2: MOR lithosphere ########################################################################\n",
       "\n",
       "# Subsidence at the MOR ###############################################################\n",
       "========= SUBSIDENCE AT THE RIDGE WITH FIXED WATER DEPTH ===========\n",
       "WATER DEPTH        = 2750.0\n",
-      "SEA LEVEL          = -412.0041304694064\n",
-      "SEDIMENT THICKNESS = 1e-09\n",
-      "dh                 = -3162.0041304704064\n",
-      "========= SUBSIDENCE AT THE RIDGE WITH FIXED SEA LEVEL ===========\n",
-      "WATER DEPTH        = 2767.7238304136886 m\n",
-      "SEA LEVEL          = -400.0 m\n",
-      "SEDIMENT THICKNESS = 1e-09 m\n",
-      "dh                 = -3167.7238304146886 m\n",
-      "========= SUBSIDENCE AT THE RIDGE WITH FIXED SEA LEVEL ===========\n",
+      "SEA LEVEL          = -395.7671958263859\n",
+      "dh                 = -3145.767195827386\n",
+      "========= SUBSIDENCE AT THE RIDGE WITH FIXED SEA LEVEL =========== (semi-analytical solution of the 2-D thermo-mechanical model without water-load)\n",
       "WATER DEPTH        = 0.0 m\n",
       "SEA LEVEL          = -10000.0 m (=> no water load)\n",
-      "SEDIMENT THICKNESS = 1e-09 m\n",
-      "dh                 = -2274.436383075561 m\n"
+      "dh                 = -2301.589291048629 m\n",
+      "========= SUBSIDENCE AT THE RIDGE WITH FIXED SEA LEVEL =========== (semi-analytical solution of the 2-D thermo-mechanical model with water-load)\n",
+      "WATER DEPTH        = 2743.890079699708 m\n",
+      "SEA LEVEL          = -400.0 m\n",
+      "dh                 = -3143.890079700708 m\n"
      ]
     }
    ],
    "source": [
-    "# 1-D isostatic comparison of reference columns\n",
+    "# 1-D isostatic comparison of reference columns, provide semi-analytical solution for the water depth at the ridge for models M1 and M2 (Figures 7 and 8)\n",
     "#\n",
     "# 1) definition of the reference 125 km thick continental lithosphere\n",
     "# 2) definition of the MOR column (model M2)\n",
+    "# 3) definition of reference 160, 180, 200, 240 and 280 km thick lithospheres\n",
     "#\n",
     "# Remark:\n",
     "# sea level and dh (elevation difference between continent and MOR)\n",
     "# are given relative to the top of the continental column, i.e. negative = down\n",
     "#\n",
     "water_depth_at_ridge     = 2750.0e0  # meters = reference water depth at the ridge based on average MOR elevation\n",
-    "sediment_thickness_fixed = 0.0e0     # meters\n",
     "sea_level_ref            = -400e0    # meters = reference sea level based on average elevation of continents\n",
     "dz                       = 200e0     # meters\n",
     "\n",
@@ -83,12 +80,12 @@
     "book_reflith             = book_lith125_GL # updated parameters (./scripts/book_lith.py)\n",
     "\n",
     "# Lithosphere 2 = Compared column\n",
-    "# The MOR from 2-D thermechanical model is used to calibrate the elevation of the reference continental column\n",
-    "# When you compare two continental column, it should be no elevation difference\n",
+    "# The MOR from 2-D thermechanical model is used to calibrate density structure of the reference continental column \n",
+    "# When you compare two continental columns, it should be no elevation differences!!\n",
     "#\n",
     "complith                 = 'MOR'  # 'MOR' 'lith160' 'lith180' 'lith200' 'lith240' 'lith280'\n",
     "\n",
-    "use_perplex_table_mantle = False\n",
+    "use_perplex_table_mantle = True\n",
     "perplex_name             = 'slm_sp2008'   # composition of the continental lithospheric mantle\n",
     "\n",
     "# Select right set of parameters, updated parameters (./scripts/book_lith.py)\n",
@@ -118,11 +115,11 @@
     "    # - define the correct perplex grid\n",
     "    # - use the right degree of depletion in the continental lithospheric mantle to fit MOR/continets elevation\n",
     "    #\n",
-    "    book_reflith  = update_mantle_density(reflith, book_reflith, perplex_name=perplex_name)  \n",
-    "    book_complith = update_mantle_density(complith,book_complith,perplex_name=perplex_name)\n",
+    "    book_reflith  = update_mantle_density_to_reference(reflith, book_reflith, perplex_name=perplex_name)  \n",
+    "    book_complith = update_mantle_density_to_reference(complith,book_complith,perplex_name=perplex_name)\n",
     "else:\n",
-    "    book_reflith  = update_mantle_density(reflith, book_reflith, perplex_name=None)  \n",
-    "    book_complith = update_mantle_density(complith,book_complith,perplex_name=None)\n",
+    "    book_reflith  = update_mantle_density_to_reference(reflith, book_reflith, perplex_name=None)  \n",
+    "    book_complith = update_mantle_density_to_reference(complith,book_complith,perplex_name=None)\n",
     "\n",
     "    \n",
     "print(\"# Column 1: Reference Continental lithosphere ###############################################################\")\n",
@@ -140,16 +137,16 @@
     "print(\"\\n# Column 2: \"+complith+\" lithosphere ########################################################################\")\n",
     "if complith == 'MOR':\n",
     "    if use_perplex_table_mantle:\n",
-    "        deltarho_ridge_correction        = -0.2 # kg/m3 # correction to reproduce the elevation difference \n",
-    "                                                # in the 2-D thermo-mechanical model\n",
-    "        geotherm_filename                = 'profile_temperature_ridge_centralvalley_GL-Tp1280_speed3-perplex_corrected.tz'\n",
-    "        density_profile_filename_ridge   = 'profile_density_ridge_centralvalley_GL-Tp1280_speed3-perplex.rhoz'\n",
+    "        deltarho_ridge_correction        = -0.42 # kg/m3 # correction to reproduce the elevation difference \n",
+    "                                                 # in the 2-D thermo-mechanical model M1\n",
+    "        geotherm_filename                = 'profile_temperature_ridge_centralvalley_GL-Tp1280_speed3-perplex_M1.tz'\n",
+    "        density_profile_filename_ridge   = 'profile_density_ridge_centralvalley_GL-Tp1280_speed3-perplex_M1.rhoz'\n",
     "        path                             = './data/geodyn2d/M1/'\n",
     "    else:\n",
-    "        deltarho_ridge_correction        = 4.5  # kg/m3 # correction to reproduce the elevation difference \n",
-    "                                                # in the 2-D thermo-mechanical model\n",
-    "        geotherm_filename                = 'profile_temperature_ridge_centralvalley_GL-Tp1280_speed3_corrected.tz'\n",
-    "        density_profile_filename_ridge   = 'profile_density_ridge_centralvalley_GL-Tp1280_speed3.rhoz'\n",
+    "        deltarho_ridge_correction        = 4.05 # kg/m3 # correction to reproduce the elevation difference \n",
+    "                                                        # in the 2-D thermo-mechanical model M2\n",
+    "        geotherm_filename                = 'profile_temperature_ridge_centralvalley_GL-Tp1280_speed3_M2.tz'\n",
+    "        density_profile_filename_ridge   = 'profile_density_ridge_centralvalley_GL-Tp1280_speed3_M2.rhoz'\n",
     "        path                             = './data/geodyn2d/M2/'\n",
     "\n",
     "    water_depth_used = water_depth_at_ridge\n",
@@ -173,7 +170,7 @@
     "    \n",
     "print(\"\\n# Subsidence at the \"+complith+\" ###############################################################\")\n",
     "#-------------------------------------------------------------------\n",
-    "lith2.set_sediment_thickness(sediment_thickness_fixed)\n",
+    "lith2.set_sediment_thickness(0)\n",
     "lith2.set_water_depth(water_depth_used)\n",
     "\n",
     "dh,water_depth,sea_level,sediment_thickness = lith2.get_subsidence_with_sediments(pref,\n",
@@ -184,53 +181,49 @@
     "print('========= SUBSIDENCE AT THE RIDGE WITH FIXED WATER DEPTH ===========')\n",
     "print('WATER DEPTH        = {}'.format(water_depth))\n",
     "print('SEA LEVEL          = {}'.format(sea_level))\n",
-    "print('SEDIMENT THICKNESS = {}'.format(sediment_thickness))\n",
     "print('dh                 = {}'.format(dh))\n",
     "# Stored for the plot in the next cell\n",
     "dh_fixed_water_depth = dh\n",
     "water_depth_final    = water_depth\n",
-    "\n",
     "#-------------------------------------------------------------------\n",
-    "lith2.set_sediment_thickness(sediment_thickness_fixed)\n",
-    "lith2.set_sea_level(sea_level_ref)   #-1.e-9)\n",
+    "lith2.set_sediment_thickness(0)\n",
+    "lith2.set_sea_level(-10000.)\n",
     "\n",
     "dh,water_depth,sea_level,sediment_thickness = lith2.get_subsidence_with_sediments(pref,\n",
     "                                                                                  ref_depth,\n",
     "                                                                                  SeaLevelFixed=True,\n",
     "                                                                                  printOut=False)\n",
-    "    \n",
     "\n",
-    "print('========= SUBSIDENCE AT THE RIDGE WITH FIXED SEA LEVEL ===========')\n",
+    "print('========= SUBSIDENCE AT THE RIDGE WITH FIXED SEA LEVEL =========== (semi-analytical solution of the 2-D thermo-mechanical model without water-load)')\n",
     "print('WATER DEPTH        = {} m'.format(water_depth))\n",
-    "print('SEA LEVEL          = {} m'.format(sea_level))\n",
-    "print('SEDIMENT THICKNESS = {} m'.format(sediment_thickness))\n",
+    "print('SEA LEVEL          = {} m (=> no water load)'.format(sea_level))\n",
     "print('dh                 = {} m'.format(dh))\n",
-    "\n",
     "#-------------------------------------------------------------------\n",
-    "lith2.set_sediment_thickness(sediment_thickness_fixed)\n",
-    "lith2.set_sea_level(-10000.)\n",
+    "lith2.set_sediment_thickness(0)\n",
+    "lith2.set_sea_level(sea_level_ref)\n",
     "\n",
     "dh,water_depth,sea_level,sediment_thickness = lith2.get_subsidence_with_sediments(pref,\n",
     "                                                                                  ref_depth,\n",
     "                                                                                  SeaLevelFixed=True,\n",
     "                                                                                  printOut=False)\n",
+    "    \n",
     "\n",
-    "print('========= SUBSIDENCE AT THE RIDGE WITH FIXED SEA LEVEL ===========')\n",
+    "print('========= SUBSIDENCE AT THE RIDGE WITH FIXED SEA LEVEL =========== (semi-analytical solution of the 2-D thermo-mechanical model with water-load)')\n",
     "print('WATER DEPTH        = {} m'.format(water_depth))\n",
-    "print('SEA LEVEL          = {} m (=> no water load)'.format(sea_level))\n",
-    "print('SEDIMENT THICKNESS = {} m'.format(sediment_thickness))\n",
-    "print('dh                 = {} m'.format(dh))\n"
+    "print('SEA LEVEL          = {} m'.format(sea_level))\n",
+    "print('dh                 = {} m'.format(dh))\n",
+    "\n"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 3,
+   "execution_count": 23,
    "metadata": {},
    "outputs": [
     {
      "data": {
       "application/vnd.jupyter.widget-view+json": {
-       "model_id": "fb6d3af1e1d7443f95b941a414bd87f2",
+       "model_id": "8a619d40ff944a08b54cd260928d1fab",
        "version_major": 2,
        "version_minor": 0
       },
@@ -244,7 +237,7 @@
     {
      "data": {
       "application/vnd.jupyter.widget-view+json": {
-       "model_id": "a96f54eb95dd4258b54921eef6a7f6ec",
+       "model_id": "16f4421e434c48ac8a4766d51acc5747",
        "version_major": 2,
        "version_minor": 0
       },
@@ -258,7 +251,7 @@
     {
      "data": {
       "application/vnd.jupyter.widget-view+json": {
-       "model_id": "61532d9975d2486abddcf493fba52d9a",
+       "model_id": "9587e4fe57eb450e8fe88d142f675591",
        "version_major": 2,
        "version_minor": 0
       },

diff --git a/scripts/book_deltarho.py b/scripts/book_deltarho.py
@@ -48,6 +48,7 @@
         },
     'onlytempdep': {
         'MOR':-12,
+        'MOR2':-12,
         'lith125':-20,
         'lith160':-28.5,
         'lith180':-31.5,

diff --git a/scripts/external.py b/scripts/external.py
@@ -425,9 +425,14 @@ def load_etopo(path='./data/topography/',filtered=True,resolution=2,corrected_fr
         else:
             filename = 'ETOPO1_Bed_g_gmt4_2m.grd'
         if ( not os.path.isfile(path+filename)):
-            print("Unfiltered 2 arc-min ETOPO1 dataset not available (check if you have downloaded FigShare Dataset) - Not available on Binder ({})".format(path+filename))
-            x=[] ; y=[] ; elev=[]
-            return x,y,elev
+            if filename == 'ETOPO1_Bed_g_gmt4.grd':
+                print("Unfiltered raw 1 arc-min ETOPO1 dataset not available (check if you have downloaded FigShare Dataset) - Not available on Binder ({})".format(path+filename))
+                x=[] ; y=[] ; elev=[]
+                return x,y,elev
+            name      = filename.split('.grd')[0]
+            nc_etopo1 = xr.open_mfdataset([path+name+'_1.grd',path+name+'_2.grd',path+name+'_3.grd',path+name+'_4.grd',path+name+'_5.grd',path+name+'_6.grd'],
+                                          concat_dim=['lon'], combine='nested',engine='netcdf4')
+            collection = True
         else:
             nc_etopo1 = netCDF4.Dataset(path+filename)
     #print(nc_etopo1.variables.keys())
@@ -438,7 +443,7 @@ def load_etopo(path='./data/topography/',filtered=True,resolution=2,corrected_fr
         x=x.to_numpy()
         y=y.to_numpy()
         elev=elev.to_numpy()
-    print("ETOPO 1 m min/max {} {}".format(np.nanmin(elev),np.nanmax(elev)))
+    print("ETOPO 1 m ({}) min/max {} {}".format(filename,np.nanmin(elev),np.nanmax(elev)))
     return x,y,elev
 
 def get_shape_etopo(path='./data/topography/'):

diff --git a/scripts/geodyn1d/isostasy/isostasy.py b/scripts/geodyn1d/isostasy/isostasy.py
@@ -105,12 +105,12 @@ def compute_pressure_density_profile(layers,geotherm,density_type=1,printOut=Tru
             p        = 0.0e0
             if filename:
                 rho_used = rho_read[0]
-                if (z>zdrho[0] and z<zdrho[1]):
-                    rho_used = rho_used + drho
                 melt_fraction_used = 0 #TODO melt_fraction_read[0] need to be added inside the file
             else:
                 rho_used = density().get_value(ilayer=ilayer,depth=z,pressure=p,temp=t,layers=layers,path=path)
                 melt_fraction_used = melt_fraction_obj().get_value(ilayer=ilayer,depth=z,pressure=p,temp=t,layers=layers,path=path)
+            if (z>zdrho[0] and z<zdrho[1]):
+                rho_used = rho_used + drho
             rho[nl-ni-1] = rho_used
             melt_fraction[nl-ni-1] = melt_fraction_used
             P[nl-ni-1]   = p
@@ -122,12 +122,12 @@ def compute_pressure_density_profile(layers,geotherm,density_type=1,printOut=Tru
                     rho_used = rho_read[len(z_read)-1]
                 else:
                     rho_used = pos2value_1d(z_read,rho_read,z)
-                if (z>zdrho[0] and z<zdrho[1]):
-                    rho_used = rho_used + drho
                 melt_fraction_used = 0 #TODO melt_fraction_read need to be added inside the file
             else:
                 rho_used = density().get_value(ilayer=ilayer,depth=z,pressure=p,temp=t,layers=layers,path=path)
                 melt_fraction_used = melt_fraction_obj().get_value(ilayer=ilayer,depth=z,pressure=p,temp=t,layers=layers,path=path)
+            if (z>zdrho[0] and z<zdrho[1]):
+                rho_used = rho_used + drho
             rho[nl-ni-1]       = rho_used
             melt_fraction[nl-ni-1] = melt_fraction_used
             p          = rho[nl-ni-1] * g * dz

diff --git a/scripts/geodyn1d/library.py b/scripts/geodyn1d/library.py
@@ -909,6 +909,43 @@
 }
 
 
+# =============================================================================
+ridge_Oc6_5_NoDiffLayer = {
+"numlayers": 6,
+"nature_layers": ['matUC','matMC','matLC','matLM1','matLM2','matSLM'],
+"thicknesses": [6.5e3,0.0e0,0.0e0,0e0,0e0,593.5e3],
+"thermalBc": ['thermBcUC','thermBcMC','thermBcLC','thermBcLM1','thermBcLM2','thermBcSLM'],
+"matUC": {
+        "temp_top": 0.0e0,
+        "H"  : 0.0e0,
+        "rho": 2900.e0
+        },
+"matMC": {
+        "H": 0.0e0,
+        "rho": 2900.e0
+        },
+"matLC": {
+        "H": 0.0e0,
+        "rho": 2900.e0
+        },
+"matLM1": {
+        "rho": 3300,
+        "H": 0.0e0
+        },
+"matLM2": {
+        "rho": 3300,
+        "H": 0.0e0
+        },
+"matSLM": {
+        "rho": 3300
+        }, 
+"thermBcSLM": {
+        "temp_bottom": 1793.15e0,
+        "temp_potential": 1553.15e0,
+        "q_bottom": 111.6718750e-3
+        }
+}
+
 
 # =============================================================================
 ridge_Oc6_5_DiffLayer26 = {