From cfb54082475ca3e3ee2509f1754767560d37bbe2 Mon Sep 17 00:00:00 2001 From: Haoyuan Li Date: Mon, 24 Apr 2023 20:39:19 -0700 Subject: [PATCH] Include option for plate model in the mass conserving temperature --- CHANGELOG.md | 1 + .../temperature/mass_conserving.h | 7 +- .../temperature/mass_conserving.cc | 105 ++++++++++-- ...rving_slab_temperature_plate_model_old.dat | 13 ++ ...erving_slab_temperature_plate_model_old.wb | 162 ++++++++++++++++++ .../screen-output.log | 9 + ...ing_slab_temperature_plate_model_young.dat | 15 ++ ...ving_slab_temperature_plate_model_young.wb | 162 ++++++++++++++++++ .../screen-output.log | 11 ++ .../world_builder_declarations_closed.md | 32 ++++ .../world_builder_declarations_open.md | 36 ++++ 11 files changed, 542 insertions(+), 11 deletions(-) create mode 100644 tests/gwb-dat/mass_conserving_slab_temperature_plate_model_old.dat create mode 100644 tests/gwb-dat/mass_conserving_slab_temperature_plate_model_old.wb create mode 100644 tests/gwb-dat/mass_conserving_slab_temperature_plate_model_old/screen-output.log create mode 100644 tests/gwb-dat/mass_conserving_slab_temperature_plate_model_young.dat create mode 100644 tests/gwb-dat/mass_conserving_slab_temperature_plate_model_young.wb create mode 100644 tests/gwb-dat/mass_conserving_slab_temperature_plate_model_young/screen-output.log diff --git a/CHANGELOG.md b/CHANGELOG.md index b3e2d978c..0ec3e61ed 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -11,6 +11,7 @@ This project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.htm ## [Unreleased] ### Added +- Added an option to use the plate model as the reference model for the mass conserving temperature of the slab. \[Haoyuan Li; 2024-02-02; [#471](https://github.com/GeodynamicWorldBuilder/WorldBuilder/pull/471) ### Changed 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 index bd7025a16..11980a10c 100644 --- 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,7 +110,12 @@ namespace WorldBuilder bool adiabatic_heating; std::vector>> mid_oceanic_ridges; Operations operation; - + enum ReferenceModelName + { + half_space_model, + plate_model + }; + ReferenceModelName reference_model_name; }; } // namespace Temperature } // namespace SubductingPlateModels 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 index 940c3543f..afbdb639a 100644 --- a/source/world_builder/features/subducting_plate_models/temperature/mass_conserving.cc +++ b/source/world_builder/features/subducting_plate_models/temperature/mass_conserving.cc @@ -164,6 +164,10 @@ namespace WorldBuilder "define the location of the ridge. You need to define at least one ridge." "So the an example with two ridges is " "[[[10,20],[20,30],[10,40]],[[50,10],[60,10]]]."); + + prm.declare_entry("reference model name", Types::String("half space model"), + "The type of thermal model to use in the mass conserving model of slab temperature. " + "Options are half space model and plate model"); } void @@ -212,6 +216,11 @@ namespace WorldBuilder { ridge_coordinate *= dtr; } + std::string reference_model_name_str = prm.get("reference model name"); + if (reference_model_name_str=="plate model") + reference_model_name = plate_model; + else if (reference_model_name_str=="half space model") + reference_model_name = half_space_model; } double @@ -226,6 +235,8 @@ namespace WorldBuilder { const double distance_from_plane = distance_from_planes.distance_from_plane; + const int plate_model_summation_number = 100; // for the plate model + if (distance_from_plane <= max_depth && distance_from_plane >= min_depth) { @@ -320,9 +331,35 @@ namespace WorldBuilder const double average_angle = distance_from_planes.average_angle; // 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) * - std::sqrt(plate_age_sec / (thermal_diffusivity * Consts::PI)); + // This uses the integral of the half-space temperature profile + // 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 (reference_model_name == plate_model) + { + initial_heat_content = thermal_conductivity / thermal_diffusivity * + (surface_temperature - potential_mantle_temperature) * max_depth / 2.0; + for (int i = 0; i< std::floor(plate_model_summation_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; + const double temp_heat_content = thermal_conductivity / thermal_diffusivity * + (surface_temperature - potential_mantle_temperature) * + 4 * max_depth / 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); + initial_heat_content -= temp_heat_content; + } + } + else + { + initial_heat_content = 2 * thermal_conductivity * (surface_temperature - potential_mantle_temperature) * + std::sqrt(plate_age_sec / (thermal_diffusivity * Consts::PI)); + } // Plate age increases with distance along the slab in the mantle double effective_plate_age = plate_age_sec + (distance_along_plane / plate_velocity) * seconds_in_year; // m/(m/y) = y(seconds_in_year) @@ -445,10 +482,33 @@ namespace WorldBuilder { // 3. Determine the heat content for side 1 (bottom) of the slab - // Comes from integrating the half-space cooling model temperature - - const double bottom_heat_content = 2 * thermal_conductivity * (min_temperature - potential_mantle_temperature) * - std::sqrt(effective_plate_age /(thermal_diffusivity * Consts::PI)); + // Comes from integrating the half-space cooling model or the plate model temperature + // 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 (reference_model_name == plate_model) + { + bottom_heat_content = thermal_conductivity / thermal_diffusivity * + (min_temperature - potential_mantle_temperature) * max_depth / 2.0; + for (int i = 0; i< std::floor(plate_model_summation_number/2.0); ++i) + { + // because n < sommation_number + 1 and n = 2k + 1 + const double plate_velocity_UI = plate_velocity / seconds_in_year; + const double temp_heat_content = thermal_conductivity / thermal_diffusivity * + (min_temperature - potential_mantle_temperature) * + 4 * max_depth / double(2*i + 1) / double(2*i + 1) / Consts::PI / Consts::PI * + exp((plate_velocity_UI * max_depth / 2.0 / 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); + bottom_heat_content -= temp_heat_content; + } + } + 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. // Should not be a positive value @@ -486,9 +546,34 @@ namespace WorldBuilder } else { - // use half-space cooling model for the bottom (side 1) of the slab - temperature = background_temperature + (min_temperature - background_temperature) * - std::erfc(adjusted_distance / (2 * std::sqrt(thermal_diffusivity * effective_plate_age))); + // use half-space cooling or plate model for the bottom (side 1) of the slab + if (reference_model_name == plate_model) + { + if (adjusted_distance < max_depth) + { + const double plate_velocity_UI = plate_velocity / seconds_in_year; + temperature = background_temperature + (min_temperature - background_temperature) * (1 - adjusted_distance / max_depth); + for (int i = 1; i< std::floor(plate_model_summation_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; + } + } + 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/mass_conserving_slab_temperature_plate_model_old.dat b/tests/gwb-dat/mass_conserving_slab_temperature_plate_model_old.dat new file mode 100644 index 000000000..1cd727f1d --- /dev/null +++ b/tests/gwb-dat/mass_conserving_slab_temperature_plate_model_old.dat @@ -0,0 +1,13 @@ +# Profiles across a slab dipping at 30 degrees +# Now define parameters: +# dim = 2 +# compositions = 0 +# x z d T +2490e3 0 5e3 +2495e3 0 5e3 +2460e3 0 25e3 +2475e3 0 35e3 +7553e3 0 145e3 +7714e3 0 44e3 +7804e3 0 126e3 +7875e3 0 150e3 \ No newline at end of file diff --git a/tests/gwb-dat/mass_conserving_slab_temperature_plate_model_old.wb b/tests/gwb-dat/mass_conserving_slab_temperature_plate_model_old.wb new file mode 100644 index 000000000..38e826083 --- /dev/null +++ b/tests/gwb-dat/mass_conserving_slab_temperature_plate_model_old.wb @@ -0,0 +1,162 @@ +{ + "version":"0.6", + "coordinate system":{"model":"cartesian"}, + "gravity model":{"model":"uniform", "magnitude":10}, + "cross section":[[0,0],[10000e3,0]],"surface temperature":273, "force surface temperature":true, + "potential mantle temperature":1673, "thermal expansion coefficient":3.1e-5, + "specific heat":1000, "thermal diffusivity":1.0e-6, + "features": + [ + { + "model": "oceanic plate", + "name": "sp plate", + "max depth": 150000.0, + "min depth": -100000.0, + "coordinates": [ + [ + 0.0, + -1000.0 + ], + [ + 0.0, + 1000.0 + ], + [ + 7500000.0, + 1000.0 + ], + [ + 7500000.0, + -1000.0 + ] + ], + "temperature models": [ + { + "model": "plate model", + "min depth": -10000.0, + "max depth": 150000.0, + "spreading velocity": 0.05, + "ridge coordinates": [ + [ + [ + 0, + -1000.0 + ], + [ + 0, + 1000.0 + ] + ] + ] + } + ], + "composition models": [ + { + "model": "uniform", + "min depth": -10000.0, + "max depth": 7500.0, + "compositions": [ + 0 + ] + }, + { + "model": "uniform", + "min depth": 7500.0, + "max depth": 35200.0, + "compositions": [ + 1 + ] + } + ] + }, + { + "model": "subducting plate", + "name": "initial slab", + "coordinates": [ + [ + 7500000.0, + -1000.0 + ], + [ + 7500000.0, + 1000.0 + ] + ], + "dip point": [ + 40000000.0, + 0.0 + ], + "segments": [ + { + "length": 418880.0, + "thickness": [ + 300000.0 + ], + "top truncation": [ + -100000.0 + ], + "angle": [ + 0, + 60 + ], + "composition models": [ + { + "model": "uniform", + "compositions": [ + 0 + ], + "max distance slab top": 7500.0 + }, + { + "model": "uniform", + "compositions": [ + 1 + ], + "min distance slab top": 7500.0, + "max distance slab top": 35200.0 + } + ] + }, + { + "length": 100000.0, + "thickness": [ + 300000.0 + ], + "top truncation": [ + -100000.0 + ], + "angle": [ + 60, + 60 + ] + } + ], + "temperature models": [ + { + "model": "mass conserving", + "density": 3300, + "thermal conductivity": 3.3, + "adiabatic heating": true, + "plate velocity": 0.05, + "ridge coordinates": [ + [ + [ + 0, + -1000.0 + ], + [ + 0, + 1000.0 + ] + ] + ], + "coupling depth": 50000.0, + "taper distance": 100000.0, + "min distance slab top": -100000.0, + "max distance slab top": 150000.0, + "reference model name": "plate model" + } + ] + } + ] +} diff --git a/tests/gwb-dat/mass_conserving_slab_temperature_plate_model_old/screen-output.log b/tests/gwb-dat/mass_conserving_slab_temperature_plate_model_old/screen-output.log new file mode 100644 index 000000000..05ad5a889 --- /dev/null +++ b/tests/gwb-dat/mass_conserving_slab_temperature_plate_model_old/screen-output.log @@ -0,0 +1,9 @@ +# x z d T +2490e3 0 5e3 391.412 +2495e3 0 5e3 391.313 +2460e3 0 25e3 775.586 +2475e3 0 35e3 948.682 +7553e3 0 145e3 1670.58 +7714e3 0 44e3 1695.98 +7804e3 0 126e3 1689.45 +7875e3 0 150e3 1752.63 diff --git a/tests/gwb-dat/mass_conserving_slab_temperature_plate_model_young.dat b/tests/gwb-dat/mass_conserving_slab_temperature_plate_model_young.dat new file mode 100644 index 000000000..3bbd2a065 --- /dev/null +++ b/tests/gwb-dat/mass_conserving_slab_temperature_plate_model_young.dat @@ -0,0 +1,15 @@ +# Profiles across a slab dipping at 30 degrees +# Now define parameters: +# dim = 2 +# compositions = 0 +# x z d T +590e3 0 5e3 +595e3 0 5e3 +660e3 0 25e3 +675e3 0 35e3 +2071e3 0 88e3 +2160e3 0 87e3 +2232e3 0 89e3 +2268e3 0 90e3 +2339e3 0 192e3 +2357e3 0 263e3 \ No newline at end of file diff --git a/tests/gwb-dat/mass_conserving_slab_temperature_plate_model_young.wb b/tests/gwb-dat/mass_conserving_slab_temperature_plate_model_young.wb new file mode 100644 index 000000000..cc3b53e47 --- /dev/null +++ b/tests/gwb-dat/mass_conserving_slab_temperature_plate_model_young.wb @@ -0,0 +1,162 @@ +{ + "version":"0.6", + "coordinate system":{"model":"cartesian"}, + "gravity model":{"model":"uniform", "magnitude":10}, + "cross section":[[0,0],[10000e3,0]],"surface temperature":273, "force surface temperature":true, + "potential mantle temperature":1673, "thermal expansion coefficient":3.1e-5, + "specific heat":1000, "thermal diffusivity":1.0e-6, + "features": + [ + { + "model": "oceanic plate", + "name": "sp plate", + "max depth": 150000.0, + "min depth": -100000.0, + "coordinates": [ + [ + 0.0, + -1000.0 + ], + [ + 0.0, + 1000.0 + ], + [ + 2000000.0, + 1000.0 + ], + [ + 2000000.0, + -1000.0 + ] + ], + "temperature models": [ + { + "model": "plate model", + "min depth": -10000.0, + "max depth": 150000.0, + "spreading velocity": 0.05, + "ridge coordinates": [ + [ + [ + 0, + -1000.0 + ], + [ + 0, + 1000.0 + ] + ] + ] + } + ], + "composition models": [ + { + "model": "uniform", + "min depth": -10000.0, + "max depth": 7500.0, + "compositions": [ + 0 + ] + }, + { + "model": "uniform", + "min depth": 7500.0, + "max depth": 35200.0, + "compositions": [ + 1 + ] + } + ] + }, + { + "model": "subducting plate", + "name": "initial slab", + "coordinates": [ + [ + 2000000.0, + -1000.0 + ], + [ + 2000000.0, + 1000.0 + ] + ], + "dip point": [ + 40000000.0, + 0.0 + ], + "segments": [ + { + "length": 418880.0, + "thickness": [ + 300000.0 + ], + "top truncation": [ + -100000.0 + ], + "angle": [ + 0, + 60 + ], + "composition models": [ + { + "model": "uniform", + "compositions": [ + 0 + ], + "max distance slab top": 7500.0 + }, + { + "model": "uniform", + "compositions": [ + 1 + ], + "min distance slab top": 7500.0, + "max distance slab top": 35200.0 + } + ] + }, + { + "length": 100000.0, + "thickness": [ + 300000.0 + ], + "top truncation": [ + -100000.0 + ], + "angle": [ + 60, + 60 + ] + } + ], + "temperature models": [ + { + "model": "mass conserving", + "density": 3300, + "thermal conductivity": 3.3, + "adiabatic heating": true, + "plate velocity": 0.05, + "ridge coordinates": [ + [ + [ + 0, + -1000.0 + ], + [ + 0, + 1000.0 + ] + ] + ], + "coupling depth": 50000.0, + "taper distance": 100000.0, + "min distance slab top": -100000.0, + "max distance slab top": 150000.0, + "reference model name": "plate model" + } + ] + } + ] +} diff --git a/tests/gwb-dat/mass_conserving_slab_temperature_plate_model_young/screen-output.log b/tests/gwb-dat/mass_conserving_slab_temperature_plate_model_young/screen-output.log new file mode 100644 index 000000000..30af49af0 --- /dev/null +++ b/tests/gwb-dat/mass_conserving_slab_temperature_plate_model_young/screen-output.log @@ -0,0 +1,11 @@ +# x z d T +590e3 0 5e3 494.275 +595e3 0 5e3 493.436 +660e3 0 25e3 1148.15 +675e3 0 35e3 1369.21 +2071e3 0 88e3 1544.7 +2160e3 0 87e3 1210.14 +2232e3 0 89e3 654.321 +2268e3 0 90e3 1650.18 +2339e3 0 192e3 1044.98 +2357e3 0 263e3 1553.22 diff --git a/tests/gwb-dat/world_builder_declarations_closed.md b/tests/gwb-dat/world_builder_declarations_closed.md index 4664ff566..a14563f27 100644 --- a/tests/gwb-dat/world_builder_declarations_closed.md +++ b/tests/gwb-dat/world_builder_declarations_closed.md @@ -9130,6 +9130,14 @@ :::::::::::::: +::::::::::::::{dropdown} /features/items/oneOf/5/segments/items/temperature models/items/oneOf/3/reference model name +:name: closed_features_items_oneOf_5_segments_items_temperature-models_items_oneOf_3_reference-model-name + +- **default value**:half space model +- **type**:string +- **documentation**:The type of thermal model to use in the mass conserving model of slab temperature. Options are half space model and plate model +:::::::::::::: + ::::::::::::::: @@ -10074,6 +10082,14 @@ :::::::::::::::: +::::::::::::::::{dropdown} /features/items/oneOf/5/temperature models/items/oneOf/3/reference model name +:name: closed_features_items_oneOf_5_temperature-models_items_oneOf_3_reference-model-name + +- **default value**:half space model +- **type**:string +- **documentation**:The type of thermal model to use in the mass conserving model of slab temperature. Options are half space model and plate model +:::::::::::::::: + ::::::::::::::::: @@ -11122,6 +11138,14 @@ :::::::::::: +::::::::::::{dropdown} /features/items/oneOf/5/sections/items/segments/items/temperature models/items/oneOf/3/reference model name +:name: closed_features_items_oneOf_5_sections_items_segments_items_temperature-models_items_oneOf_3_reference-model-name + +- **default value**:half space model +- **type**:string +- **documentation**:The type of thermal model to use in the mass conserving model of slab temperature. Options are half space model and plate model +:::::::::::: + ::::::::::::: @@ -12066,6 +12090,14 @@ :::::::::::::: +::::::::::::::{dropdown} /features/items/oneOf/5/sections/items/temperature models/items/oneOf/3/reference model name +:name: closed_features_items_oneOf_5_sections_items_temperature-models_items_oneOf_3_reference-model-name + +- **default value**:half space model +- **type**:string +- **documentation**:The type of thermal model to use in the mass conserving model of slab temperature. Options are half space model and plate model +:::::::::::::: + ::::::::::::::: diff --git a/tests/gwb-dat/world_builder_declarations_open.md b/tests/gwb-dat/world_builder_declarations_open.md index 214908111..2637c952f 100644 --- a/tests/gwb-dat/world_builder_declarations_open.md +++ b/tests/gwb-dat/world_builder_declarations_open.md @@ -10297,6 +10297,15 @@ :::::::::::::: +::::::::::::::{dropdown} /features/items/oneOf/5/segments/items/temperature models/items/oneOf/3/reference model name +:open: +:name: open_features_items_oneOf_5_segments_items_temperature-models_items_oneOf_3_reference-model-name + +- **default value**:half space model +- **type**:string +- **documentation**:The type of thermal model to use in the mass conserving model of slab temperature. Options are half space model and plate model +:::::::::::::: + ::::::::::::::: @@ -11350,6 +11359,15 @@ :::::::::::::::: +::::::::::::::::{dropdown} /features/items/oneOf/5/temperature models/items/oneOf/3/reference model name +:open: +:name: open_features_items_oneOf_5_temperature-models_items_oneOf_3_reference-model-name + +- **default value**:half space model +- **type**:string +- **documentation**:The type of thermal model to use in the mass conserving model of slab temperature. Options are half space model and plate model +:::::::::::::::: + ::::::::::::::::: @@ -12522,6 +12540,15 @@ :::::::::::: +::::::::::::{dropdown} /features/items/oneOf/5/sections/items/segments/items/temperature models/items/oneOf/3/reference model name +:open: +:name: open_features_items_oneOf_5_sections_items_segments_items_temperature-models_items_oneOf_3_reference-model-name + +- **default value**:half space model +- **type**:string +- **documentation**:The type of thermal model to use in the mass conserving model of slab temperature. Options are half space model and plate model +:::::::::::: + ::::::::::::: @@ -13575,6 +13602,15 @@ :::::::::::::: +::::::::::::::{dropdown} /features/items/oneOf/5/sections/items/temperature models/items/oneOf/3/reference model name +:open: +:name: open_features_items_oneOf_5_sections_items_temperature-models_items_oneOf_3_reference-model-name + +- **default value**:half space model +- **type**:string +- **documentation**:The type of thermal model to use in the mass conserving model of slab temperature. Options are half space model and plate model +:::::::::::::: + :::::::::::::::