Skip to content

Commit

Permalink
Include option for plate model in the mass conserving temperature
Browse files Browse the repository at this point in the history
  • Loading branch information
@lhy11009
lhy11009 committed Feb 2, 2024
1 parent 40655ff commit cfb5408
Show file tree
Hide file tree
Showing 11 changed files with 542 additions and 11 deletions.
1 change: 1 addition & 0 deletions CHANGELOG.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
7 changes: 6 additions & 1 deletion include/world_builder/features/subducting_plate_models/temperature/mass_conserving.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -110,7 +110,12 @@ namespace WorldBuilder
bool adiabatic_heating;
std::vector<std::vector<Point<2>>> mid_oceanic_ridges;
Operations operation;

enum ReferenceModelName
{
half_space_model,
plate_model
};
ReferenceModelName reference_model_name;
};
} // namespace Temperature
} // namespace SubductingPlateModels
Expand Down
105 changes: 95 additions & 10 deletions source/world_builder/features/subducting_plate_models/temperature/mass_conserving.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -212,6 +216,11 @@ namespace WorldBuilder
{
ridge_coordinate *= dtr;
}
std::string reference_model_name_str = prm.get<std::string>("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
Expand All @@ -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)
{
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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;

Check warning on line 568 in source/world_builder/features/subducting_plate_models/temperature/mass_conserving.cc

View check run for this annotation

Codecov / codecov/patch

source/world_builder/features/subducting_plate_models/temperature/mass_conserving.cc#L568 

Added line #L568 was not covered by tests
}
}
else
{
temperature = background_temperature + (min_temperature - background_temperature) *
std::erfc(adjusted_distance / (2 * std::sqrt(thermal_diffusivity * effective_plate_age)));

}
}
}
else
Expand Down
13 changes: 13 additions & 0 deletions tests/gwb-dat/mass_conserving_slab_temperature_plate_model_old.dat
View file
Original file line number Diff line number Diff line change
@@ -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
162 changes: 162 additions & 0 deletions tests/gwb-dat/mass_conserving_slab_temperature_plate_model_old.wb
View file
Original file line number Diff line number Diff line change
@@ -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"
}
]
}
]
}
9 changes: 9 additions & 0 deletions tests/gwb-dat/mass_conserving_slab_temperature_plate_model_old/screen-output.log
View file
Original file line number Diff line number Diff line change
@@ -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
15 changes: 15 additions & 0 deletions tests/gwb-dat/mass_conserving_slab_temperature_plate_model_young.dat
View file
Original file line number Diff line number Diff line change
@@ -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
Loading

0 comments on commit cfb5408

Please sign in to comment.