Include option for plate model in the mass conserving temperature

GeodynamicWorldBuilder · Feb 2, 2023 · a4ab4ff · a4ab4ff
1 parent 1f9e88e
commit a4ab4ff
Show file tree

Hide file tree

Showing 4 changed files with 141 additions and 4 deletions.
diff --git a/include/world_builder/features/subducting_plate_models/temperature/mass_conserving.h b/include/world_builder/features/subducting_plate_models/temperature/mass_conserving.h
@@ -110,6 +110,8 @@ namespace WorldBuilder
             bool adiabatic_heating;
             std::vector<std::vector<Point<2>>> mid_oceanic_ridges;
             Utilities::Operations operation;
+            //todo: a variable to decide whether half space or plate model is used
+            bool todo_var;
 
         };
       } // namespace Temperature

diff --git a/source/world_builder/features/subducting_plate_models/temperature/mass_conserving.cc b/source/world_builder/features/subducting_plate_models/temperature/mass_conserving.cc
@@ -152,6 +152,12 @@ namespace WorldBuilder
                             "define the location of the ridge. You need to define at least one ridge."
                             "So the an exmple with two ridges is "
                             "[[[10,20],[20,30],[10,40]],[[50,10],[60,10]]].");
+
+          prm.declare_entry("todo variables",  Types::Bool(false),
+                            "The type of model to use in the mass conserving model of slab temperature"
+          );
+
+
         }
 
         void
@@ -200,6 +206,7 @@ namespace WorldBuilder
               {
                 ridge_coordinate *= dtr;
               }
+          todo_var = false;
         }
 
         double
@@ -216,6 +223,8 @@ namespace WorldBuilder
           const double distance_from_plane = distance_from_planes.distance_from_plane;
           const double distance_along_plane = distance_from_planes.distance_along_plane;
           const double total_segment_length = additional_parameters.total_local_segment_length;
+          const int plate_model_sommation_number = 100; // for the plate model
+
 
           if (distance_from_plane <= max_depth && distance_from_plane >= min_depth)
             {
@@ -302,8 +311,30 @@ namespace WorldBuilder
 
               // 1. Determine initial heat content of the slab based on age of plate at trench
               //    This uses the integral of the half-space temperature profile
-              double initial_heat_content = 2 * thermal_conductivity * (surface_temperature - potential_mantle_temperature) *
+              //    The initial heat content is also decided from the type of thermal model to use in the 
+              //    mass conserving model
+              double initial_heat_content;
+              if (todo_var)
+              {
+                initial_heat_content = (surface_temperature - potential_mantle_temperature) * max_depth / 2.0;
+                for (int i = 1; i< std::floor(plate_model_sommation_number/2.0); ++i){
+                  // because n < sommation_number + 1 and n = 2k + 1
+                  // The "plate_velocity" instead of "plate_velocity_UI" is used for the last instance as "age_at_trench" has
+                  // a unit of yr.
+                  const double plate_velocity_UI = plate_velocity / seconds_in_year;
+                  initial_heat_content -= (surface_temperature - potential_mantle_temperature) *
+                                        (4 / double(2*i + 1) / double(2*i + 1) / Consts::PI / Consts::PI * 
+                                        exp((plate_velocity_UI * max_depth / 2 / thermal_diffusivity - 
+                                          std::sqrt(plate_velocity_UI * plate_velocity_UI * max_depth * max_depth / 4.0 / thermal_diffusivity / thermal_diffusivity + 
+                                            double(2*i + 1) * double(2*i + 1) * Consts::PI * Consts::PI)) *
+                                          plate_velocity * age_at_trench / max_depth));
+                }
+              }
+              else
+              {
+                initial_heat_content = 2 * thermal_conductivity * (surface_temperature - potential_mantle_temperature) *
                                             std::sqrt(plate_age_sec / (thermal_diffusivity * Consts::PI));
+              }
 
               //  2. Get Tmin and distance_offset given distance along slab and depth of point on the slab.
               //  These equations are empirical based on fitting the temperature profiles from dynamic subduction models.
@@ -396,8 +427,26 @@ namespace WorldBuilder
                       effective_plate_age = effective_plate_age * (total_segment_length - distance_along_plane)/ (taper_distance);
                     }
 
-                  const double bottom_heat_content = 2 * thermal_conductivity * (min_temperature - potential_mantle_temperature) *
-                                                     std::sqrt(effective_plate_age /(thermal_diffusivity * Consts::PI));
+                  // the bottom heat content is also decided from the type of thermal model to use in the 
+                  // mass conserving model
+                  double bottom_heat_content;
+                  if (todo_var){
+                    bottom_heat_content = (min_temperature - potential_mantle_temperature) * max_depth / 2.0;
+                    for (int i = 1; i< std::floor(plate_model_sommation_number/2.0); ++i){
+                    // because n < sommation_number + 1 and n = 2k + 1
+                    const double plate_velocity_UI = plate_velocity / seconds_in_year;
+                    bottom_heat_content -= (min_temperature - potential_mantle_temperature) *
+                                        (4 / double(2*i + 1) / double(2*i + 1) / Consts::PI / Consts::PI * 
+                                        exp((plate_velocity_UI * max_depth / 2 / thermal_diffusivity - 
+                                          std::sqrt(plate_velocity_UI * plate_velocity_UI * max_depth * max_depth / 4.0 / thermal_diffusivity / thermal_diffusivity + 
+                                            double(2*i + 1) * double(2*i + 1) * Consts::PI * Consts::PI)) *
+                                          plate_velocity_UI * effective_plate_age / max_depth));
+                    }
+                  }
+                  else{
+                    bottom_heat_content = 2 * thermal_conductivity * (min_temperature - potential_mantle_temperature) *
+                                          std::sqrt(effective_plate_age /(thermal_diffusivity * Consts::PI));
+                  } 
 
                   // 4. The difference in heat content goes into the temperature above where Tmin occurs.
                   double top_heat_content = initial_heat_content - bottom_heat_content;
@@ -428,12 +477,30 @@ namespace WorldBuilder
                         }
                       // temperature = temperature_ + (2 * top_heat_content / (2 * density * specific_heat * std::sqrt(Consts::PI * thermal_diffusivity * time_top_slab))) *
                       //              std::exp(-(adjusted_distance * adjusted_distance) / (4 * thermal_diffusivity * time_top_slab));
+                      // maybe things need to be changed here as well
                     }
                   else
                     {
                       // use half-space cooling model for the bottom (side 1) of the slab
-                      temperature = background_temperature + (min_temperature - background_temperature) *
+                      if (todo_var){
+                        const double plate_velocity_UI = plate_velocity / seconds_in_year;
+                        temperature = background_temperature;
+                        // todo: the second operator is either "+" or "-" here
+                        for (int i = 1; i< std::floor(plate_model_sommation_number/2.0); ++i){
+                            temperature = temperature - (min_temperature - background_temperature) *
+                                      ((2 / (double(i) * Consts::PI)) * std::sin((double(i) * Consts::PI * adjusted_distance) / max_depth) *
+                                       std::exp((((plate_velocity_UI * max_depth)/(2 * thermal_diffusivity)) -
+                                                 std::sqrt(((plate_velocity_UI*plate_velocity_UI*max_depth*max_depth) /
+                                                            (4*thermal_diffusivity*thermal_diffusivity)) + double(i) * double(i) * Consts::PI * Consts::PI)) *
+                                                ((plate_velocity_UI * effective_plate_age) / max_depth)));
+                        }              
+                      }
+                      else
+                      {
+                        temperature = background_temperature + (min_temperature - background_temperature) *
                                     std::erfc(adjusted_distance / (2 * std::sqrt(thermal_diffusivity * effective_plate_age)));
+
+                      }
                     }
                 }
               else

diff --git a/tests/gwb-dat/world_buider_declarations_closed.md b/tests/gwb-dat/world_buider_declarations_closed.md
@@ -9130,6 +9130,14 @@
 
 ::::::::::::::
 
+::::::::::::::{dropdown} /features/items/oneOf/5/segments/items/temperature models/items/oneOf/3/todo variables
+:name: closed_features_items_oneOf_5_segments_items_temperature-models_items_oneOf_3_todo-variables
+
+- **default value**:false
+- **type**:boolean
+- **documentation**:The type of model to use in the mass conserving model of slab temperature
+::::::::::::::
+
 
 
 :::::::::::::::
@@ -9974,6 +9982,14 @@
 
 ::::::::::::::::
 
+::::::::::::::::{dropdown} /features/items/oneOf/5/temperature models/items/oneOf/3/todo variables
+:name: closed_features_items_oneOf_5_temperature-models_items_oneOf_3_todo-variables
+
+- **default value**:false
+- **type**:boolean
+- **documentation**:The type of model to use in the mass conserving model of slab temperature
+::::::::::::::::
+
 
 
 :::::::::::::::::
@@ -10922,6 +10938,14 @@
 
 ::::::::::::
 
+::::::::::::{dropdown} /features/items/oneOf/5/sections/items/segments/items/temperature models/items/oneOf/3/todo variables
+:name: closed_features_items_oneOf_5_sections_items_segments_items_temperature-models_items_oneOf_3_todo-variables
+
+- **default value**:false
+- **type**:boolean
+- **documentation**:The type of model to use in the mass conserving model of slab temperature
+::::::::::::
+
 
 
 :::::::::::::
@@ -11766,6 +11790,14 @@
 
 ::::::::::::::
 
+::::::::::::::{dropdown} /features/items/oneOf/5/sections/items/temperature models/items/oneOf/3/todo variables
+:name: closed_features_items_oneOf_5_sections_items_temperature-models_items_oneOf_3_todo-variables
+
+- **default value**:false
+- **type**:boolean
+- **documentation**:The type of model to use in the mass conserving model of slab temperature
+::::::::::::::
+
 
 
 :::::::::::::::

diff --git a/tests/gwb-dat/world_buider_declarations_open.md b/tests/gwb-dat/world_buider_declarations_open.md
@@ -10297,6 +10297,15 @@
 
 ::::::::::::::
 
+::::::::::::::{dropdown} /features/items/oneOf/5/segments/items/temperature models/items/oneOf/3/todo variables
+:open:
+:name: open_features_items_oneOf_5_segments_items_temperature-models_items_oneOf_3_todo-variables
+
+- **default value**:false
+- **type**:boolean
+- **documentation**:The type of model to use in the mass conserving model of slab temperature
+::::::::::::::
+
 
 
 :::::::::::::::
@@ -11239,6 +11248,15 @@
 
 ::::::::::::::::
 
+::::::::::::::::{dropdown} /features/items/oneOf/5/temperature models/items/oneOf/3/todo variables
+:open:
+:name: open_features_items_oneOf_5_temperature-models_items_oneOf_3_todo-variables
+
+- **default value**:false
+- **type**:boolean
+- **documentation**:The type of model to use in the mass conserving model of slab temperature
+::::::::::::::::
+
 
 
 :::::::::::::::::
@@ -12300,6 +12318,15 @@
 
 ::::::::::::
 
+::::::::::::{dropdown} /features/items/oneOf/5/sections/items/segments/items/temperature models/items/oneOf/3/todo variables
+:open:
+:name: open_features_items_oneOf_5_sections_items_segments_items_temperature-models_items_oneOf_3_todo-variables
+
+- **default value**:false
+- **type**:boolean
+- **documentation**:The type of model to use in the mass conserving model of slab temperature
+::::::::::::
+
 
 
 :::::::::::::
@@ -13242,6 +13269,15 @@
 
 ::::::::::::::
 
+::::::::::::::{dropdown} /features/items/oneOf/5/sections/items/temperature models/items/oneOf/3/todo variables
+:open:
+:name: open_features_items_oneOf_5_sections_items_temperature-models_items_oneOf_3_todo-variables
+
+- **default value**:false
+- **type**:boolean
+- **documentation**:The type of model to use in the mass conserving model of slab temperature
+::::::::::::::
+
 
 
 :::::::::::::::