Skip to content

Commit

Permalink
Add spatial spreading rate option for mid oceanic ridges
Browse files Browse the repository at this point in the history
  • Loading branch information
danieldouglas92 committed Feb 15, 2024
1 parent d690a1d commit 2f1e963
Show file tree
Hide file tree
Showing 12 changed files with 803 additions and 516 deletions.
353 changes: 239 additions & 114 deletions doc/world_builder_declarations.schema.json

Large diffs are not rendered by default.

404 changes: 150 additions & 254 deletions doc/world_builder_declarations_closed.md

Large diffs are not rendered by default.

443 changes: 336 additions & 107 deletions doc/world_builder_declarations_open.md

Large diffs are not rendered by default.

Original file line number Diff line number Diff line change
Expand Up @@ -89,7 +89,7 @@ namespace WorldBuilder
Objects::Surface max_depth_surface;
double top_temperature;
double bottom_temperature;
std::vector<std::vector<double>> spreading_velocities;
std::pair<std::vector<double>,std::vector<double>> spreading_velocities;
std::vector<std::vector<Point<2> > > mid_oceanic_ridges;
std::vector<std::vector<double>> spreading_velocities_at_each_ridge_point;
Operations operation;
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -93,8 +93,9 @@ namespace WorldBuilder
Objects::Surface max_depth_surface;
double top_temperature;
double bottom_temperature;
double spreading_velocity;
std::pair<std::vector<double>,std::vector<double>> spreading_velocities;
std::vector<std::vector<Point<2> > > mid_oceanic_ridges;
std::vector<std::vector<double>> spreading_velocities_at_each_ridge_point;
Operations operation;

};
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -97,7 +97,8 @@ namespace WorldBuilder
double min_depth;
double max_depth;
double density;
double plate_velocity;
std::pair<std::vector<double>,std::vector<double>> plate_velocities;
std::vector<std::vector<double>> plate_velocities_at_each_ridge_point;
double mantle_coupling_depth;
double forearc_cooling_factor;
double thermal_conductivity;
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -83,7 +83,7 @@ namespace WorldBuilder
"in degree Kelvin for this feature. If the model has an adiabatic gradient"
"this should be the mantle potential temperature, and T = Tad + Thalf. ");

prm.declare_entry("spreading velocity", Types::OneOf(Types::Double(0.01),Types::Array(Types::Array(Types::Double(0.01)))),
prm.declare_entry("spreading velocity", Types::OneOf(Types::Double(0.01),Types::Array(Types::ValueAtPoints(0.01, std::numeric_limits<double>::max()))),
"The spreading velocity of the plate in meter per year. "
"This is the velocity with which one side moves away from the ridge.");

Expand All @@ -105,7 +105,7 @@ namespace WorldBuilder
operation = string_operations_to_enum(prm.get<std::string>("operation"));
top_temperature = prm.get<double>("top temperature");
bottom_temperature = prm.get<double>("bottom temperature");
spreading_velocities = prm.get_vector<std::vector<double>>("spreading velocity");
spreading_velocities = prm.get_value_at_array("spreading velocity");

mid_oceanic_ridges = prm.get_vector<std::vector<Point<2>>>("ridge coordinates");
const double dtr = prm.coordinate_system->natural_coordinate_system() == spherical ? Consts::PI / 180.0 : 1.0;
Expand All @@ -115,22 +115,22 @@ namespace WorldBuilder
ridge_coordinate *= dtr;
}

// unsigned int index_x = 0;
// unsigned int index_y = 0;
// unsigned int ridge_point_index = 0;
// for (index_x = 0; index_x < mid_oceanic_ridges.size(); index_x++)
// {
// std::vector<double> spreading_rates_for_ridge;
// for (index_y = 0; index_y < mid_oceanic_ridges[index_x].size(); index_y++)
// {
// if (spreading_velocities.second.size() <= 1)
// spreading_rates_for_ridge.push_back(spreading_velocities.first[0] / 31557600); // m/s
// else
// spreading_rates_for_ridge.push_back(spreading_velocities.second[ridge_point_index] / 31557600); // m/s
// ridge_point_index += 1;
// }
// spreading_velocities_at_each_ridge_point.push_back(spreading_rates_for_ridge);
// }
unsigned int index_x = 0;
unsigned int index_y = 0;
unsigned int ridge_point_index = 0;
for (index_x = 0; index_x < mid_oceanic_ridges.size(); index_x++)
{
std::vector<double> spreading_rates_for_ridge;
for (index_y = 0; index_y < mid_oceanic_ridges[index_x].size(); index_y++)
{
if (spreading_velocities.second.size() <= 1)
spreading_rates_for_ridge.push_back(spreading_velocities.first[0]);
else
spreading_rates_for_ridge.push_back(spreading_velocities.second[ridge_point_index]);
ridge_point_index += 1;
}
spreading_velocities_at_each_ridge_point.push_back(spreading_rates_for_ridge);
}
}

double
Expand Down Expand Up @@ -162,7 +162,7 @@ namespace WorldBuilder
}

std::pair<double, double> ridge_parameters = Utilities::calculate_ridge_distance_and_spreading(mid_oceanic_ridges,
spreading_velocities,
spreading_velocities_at_each_ridge_point,
world->parameters.coordinate_system,
position_in_natural_coordinates_at_min_depth);

Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -48,7 +48,6 @@ namespace WorldBuilder
max_depth(NaN::DSNAN),
top_temperature(NaN::DSNAN),
bottom_temperature(NaN::DSNAN),
spreading_velocity(NaN::DSNAN),
operation(Operations::REPLACE)
{
this->world = world_;
Expand Down Expand Up @@ -102,7 +101,7 @@ namespace WorldBuilder
operation = string_operations_to_enum(prm.get<std::string>("operation"));
top_temperature = prm.get<double>("top temperature");
bottom_temperature = prm.get<double>("bottom temperature");
spreading_velocity = prm.get<double>("spreading velocity")/31557600;
spreading_velocities = prm.get_value_at_array("spreading velocity");

mid_oceanic_ridges = prm.get_vector<std::vector<Point<2>>>("ridge coordinates");
const double dtr = prm.coordinate_system->natural_coordinate_system() == spherical ? Consts::PI / 180.0 : 1.0;
Expand All @@ -111,6 +110,23 @@ namespace WorldBuilder
{
ridge_coordinate *= dtr;
}

unsigned int index_x = 0;
unsigned int index_y = 0;
unsigned int ridge_point_index = 0;
for (index_x = 0; index_x < mid_oceanic_ridges.size(); index_x++)
{
std::vector<double> spreading_rates_for_ridge;
for (index_y = 0; index_y < mid_oceanic_ridges[index_x].size(); index_y++)
{
if (spreading_velocities.second.size() <= 1)
spreading_rates_for_ridge.push_back(spreading_velocities.first[0]);
else
spreading_rates_for_ridge.push_back(spreading_velocities.second[ridge_point_index]);
ridge_point_index += 1;
}
spreading_velocities_at_each_ridge_point.push_back(spreading_rates_for_ridge);
}
}


Expand Down Expand Up @@ -145,7 +161,7 @@ namespace WorldBuilder
const int sommation_number = 100;

std::pair<double, double> ridge_parameters = Utilities::calculate_ridge_distance_and_spreading(mid_oceanic_ridges,
spreading_velocity,
spreading_velocities_at_each_ridge_point,
world->parameters.coordinate_system,
position_in_natural_coordinates_at_min_depth);

Expand All @@ -159,23 +175,23 @@ namespace WorldBuilder
{
temperature = temperature + (bottom_temperature_local - top_temperature) *
((2 / (double(i) * Consts::PI)) * std::sin((double(i) * Consts::PI * depth) / max_depth) *
std::exp((((spreading_velocity * max_depth)/(2 * thermal_diffusivity)) -
std::sqrt(((spreading_velocity*spreading_velocity*max_depth*max_depth) /
std::exp((((ridge_parameters.first * max_depth)/(2 * thermal_diffusivity)) -
std::sqrt(((ridge_parameters.first*ridge_parameters.first*max_depth*max_depth) /
(4*thermal_diffusivity*thermal_diffusivity)) + double(i) * double(i) * Consts::PI * Consts::PI)) *
((spreading_velocity * age) / max_depth)));
((ridge_parameters.first * age) / max_depth)));

}

WBAssert(!std::isnan(temperature), "Temparture inside plate model is not a number: " << temperature
<< ". Relevant variables: bottom_temperature_local = " << bottom_temperature_local
<< ", top_temperature = " << top_temperature
<< ", max_depth = " << max_depth
<< ", spreading_velocity = " << spreading_velocity
<< ", spreading_velocity = " << ridge_parameters.first
<< ", thermal_diffusivity = " << thermal_diffusivity
<< ", age = " << age << '.');
WBAssert(std::isfinite(temperature), "Temparture inside plate model is not a finite: " << temperature << ". Relevant variables: bottom_temperature_local = " << bottom_temperature_local
<< ", top_temperature = " << top_temperature
<< ", spreading_velocity = " << spreading_velocity
<< ", spreading_velocity = " << ridge_parameters.first
<< ", thermal_diffusivity = " << thermal_diffusivity
<< ", age = " << age << '.');

Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -28,8 +28,10 @@
#include "world_builder/types/array.h"
#include "world_builder/types/bool.h"
#include "world_builder/types/double.h"
#include "world_builder/types/one_of.h"
#include "world_builder/types/object.h"
#include "world_builder/types/point.h"
#include "world_builder/types/value_at_points.h"
#include "world_builder/utilities.h"
#include "world_builder/world.h"

Expand All @@ -48,7 +50,6 @@ namespace WorldBuilder
min_depth(NaN::DSNAN),
max_depth(NaN::DSNAN),
density(NaN::DSNAN),
plate_velocity(NaN::DSNAN),
mantle_coupling_depth(NaN::DSNAN),
forearc_cooling_factor(NaN::DSNAN),
thermal_conductivity(NaN::DSNAN),
Expand Down Expand Up @@ -116,7 +117,7 @@ namespace WorldBuilder
prm.declare_entry("density", Types::Double(3300),
"The reference density of the subducting plate in $kg/m^3$");

prm.declare_entry("plate velocity", Types::Double(0.05),
prm.declare_entry("plate velocity", Types::OneOf(Types::Double(0.01),Types::Array(Types::Array(Types::Double(0.01)))),
"The velocity with which the plate subducts in meters per year. Default is 5 cm/yr");

prm.declare_entry("coupling depth", Types::Double(100e3),
Expand Down Expand Up @@ -180,7 +181,7 @@ namespace WorldBuilder

density = prm.get<double>("density");
thermal_conductivity = prm.get<double>("thermal conductivity");
plate_velocity = prm.get<double>("plate velocity");
plate_velocities = prm.get_value_at_array("plate velocity");

mantle_coupling_depth = prm.get<double>("coupling depth");
forearc_cooling_factor = prm.get<double>("forearc cooling factor");
Expand Down Expand Up @@ -216,6 +217,24 @@ namespace WorldBuilder
{
ridge_coordinate *= dtr;
}

unsigned int index_x = 0;
unsigned int index_y = 0;
unsigned int ridge_point_index = 0;
for (index_x = 0; index_x < mid_oceanic_ridges.size(); index_x++)
{
std::vector<double> plate_velocities_for_ridge;
for (index_y = 0; index_y < mid_oceanic_ridges[index_x].size(); index_y++)
{
if (plate_velocities.second.size() <= 1)
plate_velocities_for_ridge.push_back(plate_velocities.first[0]);
else
plate_velocities_for_ridge.push_back(plate_velocities.second[ridge_point_index]);
ridge_point_index += 1;
}
plate_velocities_at_each_ridge_point.push_back(plate_velocities_for_ridge);
}

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;
Expand Down Expand Up @@ -244,15 +263,17 @@ namespace WorldBuilder
const Point<3> trench_point = distance_from_planes.closest_trench_point;
const Objects::NaturalCoordinate trench_point_natural = Objects::NaturalCoordinate(trench_point,
*(world->parameters.coordinate_system));

std::pair<double, double> ridge_parameters = Utilities::calculate_ridge_distance_and_spreading(mid_oceanic_ridges,
plate_velocity,
plate_velocities_at_each_ridge_point,
world->parameters.coordinate_system,
trench_point_natural);

const double km2m = 1.0e3; // 1000 m/km
const double cm2m = 100; // 100 cm/m
const double my = 1.0e6; // 1e6 y/my
const double seconds_in_year = 60.0 * 60.0 * 24.0 * 365.25; // sec/y
const double plate_velocity = ridge_parameters.first * seconds_in_year;

const double age_at_trench = ridge_parameters.second / plate_velocity; // m/(m/y) = yr
const double plate_age_sec = age_at_trench * seconds_in_year; // y --> seconds
Expand Down
10 changes: 5 additions & 5 deletions source/world_builder/utilities.cc
Original file line number Diff line number Diff line change
Expand Up @@ -1342,17 +1342,17 @@ namespace WorldBuilder

distance_ridge = std::min(distance_ridge,
coordinate_system->distance_between_points_at_same_depth(Point<3>(position_in_natural_coordinates_at_min_depth.get_coordinates(),
position_in_natural_coordinates_at_min_depth.get_coordinate_system()),
compare_point1));
position_in_natural_coordinates_at_min_depth.get_coordinate_system()),
compare_point1));

distance_ridge = std::min(distance_ridge,
coordinate_system->distance_between_points_at_same_depth(Point<3>(position_in_natural_coordinates_at_min_depth.get_coordinates(),
position_in_natural_coordinates_at_min_depth.get_coordinate_system()),
compare_point2));
position_in_natural_coordinates_at_min_depth.get_coordinate_system()),
compare_point2));
}
}
std::pair<double, double> result;
result.first = spreading_velocity_at_ridge;
result.first = spreading_velocity_at_ridge / 31557600; // m/s;
result.second = distance_ridge;
return result;
}
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -7,4 +7,3 @@
250e3 30e3 0 50e3 489.645
250e3 60e3 0 50e3 574.848
250e3 90e3 0 50e3 639.408

Original file line number Diff line number Diff line change
Expand Up @@ -7,4 +7,3 @@
6321e3 22.5 25 50e3 647.504
6321e3 22.5 27 50e3 588.45
6321e3 22.5 29 50e3 510.556

0 comments on commit 2f1e963

Please sign in to comment.