New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Climate Model Research #9

Open
astrographer opened this Issue Apr 8, 2017 · 23 comments

Comments

Projects
None yet
2 participants
@astrographer

astrographer commented Apr 8, 2017

Found this while trying to find support for my k1 + k2cos(6lat) + k3(1-lat/90) - k4dist_to_cont_edge*sin(lat) model.

Evaporation and precipitation vs. latitude.
http://www.roperld.com/science/Evaporation_PrecipitationVSLatitude.htm

This might prove useful…

@astrographer

This comment has been minimized.

Show comment
Hide comment
@astrographer

astrographer Apr 8, 2017

Along the same vein, P-E might be a better input to the NPPp element of the Miami model, although that would require different coefficients on the curve fit to existing net primary production data.

astrographer commented Apr 8, 2017

Along the same vein, P-E might be a better input to the NPPp element of the Miami model, although that would require different coefficients on the curve fit to existing net primary production data.

@davidson16807

This comment has been minimized.

Show comment
Hide comment
@davidson16807

davidson16807 Apr 15, 2017

Owner

We should be able to easily work that into the shader code - relevant code is in "scripts/view/fragment/template.glsl.c" and "scripts/view/fragment/FragmentShaders.js"

Something relevant to mention: for the last year or so I've been working on a experimental branch to overhaul the crust model to address problems while working on issue #8 and #1 (it's the "tests/field" branch if you want to try it out on your machine, if so let me know what you think). If I can get it working, it throws the door wide open for proper climate modeling. So now is a good time to think about pie-in-the-sky climate models, because they just might be feasible soon.

There's maybe just one issue I need to work out before I'm comfortable merging that branch with master. It needs to handle plate collision a bit more sensibly: continental docking, accretion, etc. It already generates sensible looking planets, but it doesn't have code in place to accrete continental crust, and I'm unsatisfied with that

Owner

davidson16807 commented Apr 15, 2017

We should be able to easily work that into the shader code - relevant code is in "scripts/view/fragment/template.glsl.c" and "scripts/view/fragment/FragmentShaders.js"

Something relevant to mention: for the last year or so I've been working on a experimental branch to overhaul the crust model to address problems while working on issue #8 and #1 (it's the "tests/field" branch if you want to try it out on your machine, if so let me know what you think). If I can get it working, it throws the door wide open for proper climate modeling. So now is a good time to think about pie-in-the-sky climate models, because they just might be feasible soon.

There's maybe just one issue I need to work out before I'm comfortable merging that branch with master. It needs to handle plate collision a bit more sensibly: continental docking, accretion, etc. It already generates sensible looking planets, but it doesn't have code in place to accrete continental crust, and I'm unsatisfied with that

@davidson16807

This comment has been minimized.

Show comment
Hide comment
@davidson16807

davidson16807 Sep 20, 2017

Owner

I can finally see how 2d fluid simulation could be done with our model. The only problem: most fluids on a real planet require 3d fluid simulation.

Possible implementation:

var FluidModeling = {};

FluidModeling.get_advected_ids = function (velocity, timestep, advected_ids) {
	var past_pos = scatch || Float32Raster(velocity.grid);
	ScalarField.add_scalar_term	(velocity.grid.pos, velocity, -timestep, 	past_pos);
	crust.grid.getNearestIds	(past_pos, 									advected_ids);
}

// get pressure of an advected fluid
// An advected fluid has regions where fluid is compressed, 
// but fluid simulation usually assumes fluid is incompressible.
// This method takes an advected fluid velocity, and
// isolates the "compressing" component from the "swirling" component.
// The "swirling" component is declared to be the actual fluid velocity.
// The "compressing" component is declared to be the pressure gradient.
// We use something called "Jacobi iteration" to isolate the two components.
// For more info: http://jamie-wong.com/2016/08/05/webgl-fluid-simulation/
FluidModeling.refine_pressure_estimate = function (advected_velocity, pressure, result) {
  result = result || Float32Raster(grid);
  VectorField.divergence(advected_velocity, result);

  //TODO: need density and timestep somewhere here
  var neighbor_count = grid.neighbor_count;
  var average_differential = grid.average_differential * grid.average_differential;
  for (var i = 0; i < result.length; i++) {
  	result[i] *= -density;
  	result[i] /= timestep;
  	result[i] /= 4;
  	result[i] *= neighbor_count[i] * average_differential;
  }

  var arrows = field.grid.arrows;
  var arrow;
  for (var i=0, li=arrows.length; i<li; ++i) {
      arrow = arrows[i];
      result[arrow[0]] += pressure[arrow[1]] ;
  }
  for (var i = 0; i < result.length; i++) {
  	result[i] /= 4;
  }
  return result;
}

FluidModeling.get_advected_ids = function (velocity, timestep, advected_ids, scratch) {
	var past_pos = scatch || VectorRaster(velocity.grid);
	ScalarField.add_scalar_term	(velocity.grid.pos, velocity, -timestep, 	past_pos);
	crust.grid.getNearestIds	(past_pos, 									advected_ids);
	return advected_ids;
}

FluidModeling.simulate = function(velocity, pressure, timestep, iterations, scratch_f32, scratch_ui8, scratch_vec) {
	// Step 1, advect velocity:
    var advected_velocity = scratch_f32;
	FluidModeling.advect 	(velocity, velocity, timestep, scratch_ui8, 	advected_velocity)
	// copy so we can free up scratch_f32
	Float32Raster.copy 		(advected_velocity, velocity);

	// Step 2, find pressure:
	var temp = pressure;
	var refined_pressure = scratch_f32;
	for (var i = 0; i < iterations; i++) {
		FluidModeling.refine_pressure_estimate(advected_velocity, pressure, refined_pressure);
		temp = pressure;
		pressure = refined_pressure;
		refined_pressure = temp;
	}
	// copy so we can free up scratch_f32
	Float32Raster.copy 		(pressure, refined_pressure);

	// Step 3, find velocity
    var pressure_gradient = scratch_vec;
    ScalarField.gradient 	(pressure, 										pressure_gradient);
	VectorRaster.sub_field	(velocity, pressure_gradient, 					velocity);

	return velocity;
}
Owner

davidson16807 commented Sep 20, 2017

I can finally see how 2d fluid simulation could be done with our model. The only problem: most fluids on a real planet require 3d fluid simulation.

Possible implementation:

var FluidModeling = {};

FluidModeling.get_advected_ids = function (velocity, timestep, advected_ids) {
	var past_pos = scatch || Float32Raster(velocity.grid);
	ScalarField.add_scalar_term	(velocity.grid.pos, velocity, -timestep, 	past_pos);
	crust.grid.getNearestIds	(past_pos, 									advected_ids);
}

// get pressure of an advected fluid
// An advected fluid has regions where fluid is compressed, 
// but fluid simulation usually assumes fluid is incompressible.
// This method takes an advected fluid velocity, and
// isolates the "compressing" component from the "swirling" component.
// The "swirling" component is declared to be the actual fluid velocity.
// The "compressing" component is declared to be the pressure gradient.
// We use something called "Jacobi iteration" to isolate the two components.
// For more info: http://jamie-wong.com/2016/08/05/webgl-fluid-simulation/
FluidModeling.refine_pressure_estimate = function (advected_velocity, pressure, result) {
  result = result || Float32Raster(grid);
  VectorField.divergence(advected_velocity, result);

  //TODO: need density and timestep somewhere here
  var neighbor_count = grid.neighbor_count;
  var average_differential = grid.average_differential * grid.average_differential;
  for (var i = 0; i < result.length; i++) {
  	result[i] *= -density;
  	result[i] /= timestep;
  	result[i] /= 4;
  	result[i] *= neighbor_count[i] * average_differential;
  }

  var arrows = field.grid.arrows;
  var arrow;
  for (var i=0, li=arrows.length; i<li; ++i) {
      arrow = arrows[i];
      result[arrow[0]] += pressure[arrow[1]] ;
  }
  for (var i = 0; i < result.length; i++) {
  	result[i] /= 4;
  }
  return result;
}

FluidModeling.get_advected_ids = function (velocity, timestep, advected_ids, scratch) {
	var past_pos = scatch || VectorRaster(velocity.grid);
	ScalarField.add_scalar_term	(velocity.grid.pos, velocity, -timestep, 	past_pos);
	crust.grid.getNearestIds	(past_pos, 									advected_ids);
	return advected_ids;
}

FluidModeling.simulate = function(velocity, pressure, timestep, iterations, scratch_f32, scratch_ui8, scratch_vec) {
	// Step 1, advect velocity:
    var advected_velocity = scratch_f32;
	FluidModeling.advect 	(velocity, velocity, timestep, scratch_ui8, 	advected_velocity)
	// copy so we can free up scratch_f32
	Float32Raster.copy 		(advected_velocity, velocity);

	// Step 2, find pressure:
	var temp = pressure;
	var refined_pressure = scratch_f32;
	for (var i = 0; i < iterations; i++) {
		FluidModeling.refine_pressure_estimate(advected_velocity, pressure, refined_pressure);
		temp = pressure;
		pressure = refined_pressure;
		refined_pressure = temp;
	}
	// copy so we can free up scratch_f32
	Float32Raster.copy 		(pressure, refined_pressure);

	// Step 3, find velocity
    var pressure_gradient = scratch_vec;
    ScalarField.gradient 	(pressure, 										pressure_gradient);
	VectorRaster.sub_field	(velocity, pressure_gradient, 					velocity);

	return velocity;
}
@davidson16807

This comment has been minimized.

Show comment
Hide comment
@davidson16807

davidson16807 Oct 1, 2017

Owner

OK, really excited now.

I'm prototyping a climate model based on what I see here: https://imgur.com/a/UNvLF#vrrcARn. Nowhere near scientifically valid (at least not yet), and probably has lots of bugs.

We start with a simple model for air pressure based on latitude. It's just a few bands of alternating high and low pressure, represented by cosine of latitude:

AtmosphericModeling.surface_air_pressure_lat_effect = function (lat, pressure) {
	pressure = pressure || Float32Raster(lat.grid);
	var cos = Math.cos;
	for (var i=0, li=lat.length; i<li; ++i) {
	    pressure[i] = -cos(5*(lat[i]));
	}
	return pressure
}
view.setScalarDisplay(new ScalarHeatDisplay( { min: '-1.', max: '1.', 
		getField: function (world, effect, scratch) {
			var axial_tilt = Math.PI*23.5/180;
			var lat = Float32SphereRaster.latitude(world.grid.pos.y);
			AtmosphericModeling.surface_air_pressure_lat_effect(lat, effect);
			return effect;
		} 
	} ));

1506828375337-nanyr 1

Offset the latitude to simulate winter (I think this is the right thing to do):

view.setScalarDisplay(new ScalarHeatDisplay( { min: '-1.', max: '1.', 
		getField: function (world, effect, scratch) {
			var axial_tilt = Math.PI*23.5/180;
			var lat = Float32SphereRaster.latitude(world.grid.pos.y);
			var effective_lat = ScalarField.add_scalar(lat, axial_tilt);
			Float32RasterInterpolation.clamp(effective_lat, -Math.PI/2, Math.PI/2, effective_lat);
			AtmosphericModeling.surface_air_pressure_lat_effect(effective_lat, effect);
			return effect;
		} 
	} ));

1506828375337-nanyr 4

Result looks boring because there's no land. So let's find out where the land is:

view.setScalarDisplay(new ScalarHeatDisplay( { min: '-1.', max: '1.', 
		getField: function (world, effect, scratch) {
			var lat = Float32SphereRaster.latitude(world.grid.pos.y);
			var is_land = ScalarField.gt_scalar(world.displacement, world.SEALEVEL);
			return is_land;
		}
	} ));

1506828375337-nanyr 2

Smooth it out with binary morphology and the diffusion equation. This part is costly, performance-wise.

view.setScalarDisplay(new ScalarHeatDisplay( { min: '-1.', max: '1.', 
		getField: function (world, effect, scratch) {
			var lat = Float32SphereRaster.latitude(world.grid.pos.y);
			var is_land = ScalarField.gt_scalar(world.displacement, world.SEALEVEL);
			BinaryMorphology.erosion(is_land, 2, is_land);
			Float32Raster.FromUint8Raster(is_land, effect);
			var diffuse = ScalarField.diffusion_by_constant;
			var smoothing_iterations =  2;
			for (var i=0; i<smoothing_iterations; ++i) {
				diffuse(effect, 1, effect, scratch);
			}
			return effect;
		}
	} ));

1506828375337-nanyr 3

Land has different pressure effects depending on hemisphere. IIRC land in winter gets higher pressure:

view.setScalarDisplay(new ScalarHeatDisplay( { min: '-1.', max: '1.', 
		getField: function (world, effect, scratch) {
			var lat = Float32SphereRaster.latitude(world.grid.pos.y);
			var is_land = ScalarField.gt_scalar(world.displacement, world.SEALEVEL);
			BinaryMorphology.erosion(is_land, 2, is_land);
			Float32Raster.FromUint8Raster(is_land, effect);
			var diffuse = ScalarField.diffusion_by_constant;
			var smoothing_iterations =  2;
			for (var i=0; i<smoothing_iterations; ++i) {
				diffuse(effect, 1, effect, scratch);
			}
			var pos = world.grid.pos;
			ScalarField.mult_field(effect, lat, effect);
			return effect;
		}
	} ));

1506828375337-nanyr

Turn this into a function in the AtmosphericModeling namespace, and add a "season" parameter (season==1 for winter in the northern hemisphere and season==-1 for summer in the northern hemisphere)

AtmosphericModeling.surface_air_pressure_land_effect = function(displacement, lat, sealevel, season, effect, scratch) {
	var is_land = ScalarField.gt_scalar(displacement, sealevel);
	BinaryMorphology.erosion(is_land, 2, is_land);
	Float32Raster.FromUint8Raster(is_land, effect);
	var diffuse = ScalarField.diffusion_by_constant;
	var smoothing_iterations =  2;
	for (var i=0; i<smoothing_iterations; ++i) {
		diffuse(effect, 1, effect, scratch);
	}
	var pos = world.grid.pos;
	ScalarField.mult_field(effect, lat, effect);
    ScalarField.mult_scalar(effect, season, effect);
	return effect;
}

OK, now we add the two effects together:

AtmosphericModeling.surface_air_pressure = function(displacement, lat, sealevel, season, axial_tilt, pressure, scratch) {
	scratch = scratch || Float32Raster(lat.grid);
	// "effective latitude" is what the latitude is like weather-wise, when taking axial tilt into account
	var effective_lat = scratch; 
	ScalarField.add_scalar(lat, season*axial_tilt, effective_lat);
	Float32RasterInterpolation.clamp(effective_lat, -Math.PI/2, Math.PI/2, effective_lat);

	AtmosphericModeling.surface_air_pressure_lat_effect(effective_lat, pressure);

	var land_effect = Float32Raster(lat.grid);
	var scratch2 = Float32Raster(lat.grid);
	AtmosphericModeling.surface_air_pressure_land_effect(displacement, effective_lat, sealevel, season, land_effect, scratch2);
	ScalarField.add_scalar_term(pressure, land_effect, 1, pressure);

	return pressure;
}

1506828375337-nanyr 5

Density times the gradient of pressure is acceleration, but we pretend it's velocity.

view.setVectorDisplay( new VectorFieldDisplay( {
		getField: function (world) {
			var lat = Float32SphereRaster.latitude(world.grid.pos.y);
			var scratch = Float32Raster(world.grid);
			var pressure = Float32Raster(world.grid);
			AtmosphericModeling.surface_air_pressure(world.displacement, lat, world.SEALEVEL, 1, Math.PI*23.5/180, pressure, scratch);
			var velocity = ScalarField.gradient(pressure);
			return velocity;
		} 
	} ));

1506828375337-nanyr 6

Now just add the coriolis effect:

AtmosphericModeling.surface_air_velocity_coriolis_effect = function(pos, velocity, angular_speed, effect) {
	effect = effect || VectorRaster(pos.grid);
	VectorField.cross_vector_field	(velocity, pos, 			effect);
	VectorField.mult_scalar 		(effect, 2 * angular_speed, effect);
	VectorField.mult_scalar_field	(effect, pos.y, 			effect);
	return effect;
}
view.setVectorDisplay( new VectorFieldDisplay( {
		getField: function (world) {
			var lat = Float32SphereRaster.latitude(world.grid.pos.y);
			var scratch = Float32Raster(world.grid);
			var pressure = Float32Raster(world.grid);
			AtmosphericModeling.surface_air_pressure(world.displacement, lat, world.SEALEVEL, 1, Math.PI*23.5/180, pressure, scratch);
			var velocity = ScalarField.gradient(pressure);
			var coriolis_effect = AtmosphericModeling.surface_air_velocity_coriolis_effect(world.grid.pos, velocity, ANGULAR_SPEED)
			VectorField.add_vector_field(velocity, coriolis_effect, velocity);
			return velocity;
		} 
	} ));

1506828375337-nanyr 7

You can see how air gets deflected around landmass:

1506828375337-nanyr 9

BTW, you should be able to run all this code in the Developer Tools console. I've also got a "atmo-model" branch out with this stuff added.

Owner

davidson16807 commented Oct 1, 2017

OK, really excited now.

I'm prototyping a climate model based on what I see here: https://imgur.com/a/UNvLF#vrrcARn. Nowhere near scientifically valid (at least not yet), and probably has lots of bugs.

We start with a simple model for air pressure based on latitude. It's just a few bands of alternating high and low pressure, represented by cosine of latitude:

AtmosphericModeling.surface_air_pressure_lat_effect = function (lat, pressure) {
	pressure = pressure || Float32Raster(lat.grid);
	var cos = Math.cos;
	for (var i=0, li=lat.length; i<li; ++i) {
	    pressure[i] = -cos(5*(lat[i]));
	}
	return pressure
}
view.setScalarDisplay(new ScalarHeatDisplay( { min: '-1.', max: '1.', 
		getField: function (world, effect, scratch) {
			var axial_tilt = Math.PI*23.5/180;
			var lat = Float32SphereRaster.latitude(world.grid.pos.y);
			AtmosphericModeling.surface_air_pressure_lat_effect(lat, effect);
			return effect;
		} 
	} ));

1506828375337-nanyr 1

Offset the latitude to simulate winter (I think this is the right thing to do):

view.setScalarDisplay(new ScalarHeatDisplay( { min: '-1.', max: '1.', 
		getField: function (world, effect, scratch) {
			var axial_tilt = Math.PI*23.5/180;
			var lat = Float32SphereRaster.latitude(world.grid.pos.y);
			var effective_lat = ScalarField.add_scalar(lat, axial_tilt);
			Float32RasterInterpolation.clamp(effective_lat, -Math.PI/2, Math.PI/2, effective_lat);
			AtmosphericModeling.surface_air_pressure_lat_effect(effective_lat, effect);
			return effect;
		} 
	} ));

1506828375337-nanyr 4

Result looks boring because there's no land. So let's find out where the land is:

view.setScalarDisplay(new ScalarHeatDisplay( { min: '-1.', max: '1.', 
		getField: function (world, effect, scratch) {
			var lat = Float32SphereRaster.latitude(world.grid.pos.y);
			var is_land = ScalarField.gt_scalar(world.displacement, world.SEALEVEL);
			return is_land;
		}
	} ));

1506828375337-nanyr 2

Smooth it out with binary morphology and the diffusion equation. This part is costly, performance-wise.

view.setScalarDisplay(new ScalarHeatDisplay( { min: '-1.', max: '1.', 
		getField: function (world, effect, scratch) {
			var lat = Float32SphereRaster.latitude(world.grid.pos.y);
			var is_land = ScalarField.gt_scalar(world.displacement, world.SEALEVEL);
			BinaryMorphology.erosion(is_land, 2, is_land);
			Float32Raster.FromUint8Raster(is_land, effect);
			var diffuse = ScalarField.diffusion_by_constant;
			var smoothing_iterations =  2;
			for (var i=0; i<smoothing_iterations; ++i) {
				diffuse(effect, 1, effect, scratch);
			}
			return effect;
		}
	} ));

1506828375337-nanyr 3

Land has different pressure effects depending on hemisphere. IIRC land in winter gets higher pressure:

view.setScalarDisplay(new ScalarHeatDisplay( { min: '-1.', max: '1.', 
		getField: function (world, effect, scratch) {
			var lat = Float32SphereRaster.latitude(world.grid.pos.y);
			var is_land = ScalarField.gt_scalar(world.displacement, world.SEALEVEL);
			BinaryMorphology.erosion(is_land, 2, is_land);
			Float32Raster.FromUint8Raster(is_land, effect);
			var diffuse = ScalarField.diffusion_by_constant;
			var smoothing_iterations =  2;
			for (var i=0; i<smoothing_iterations; ++i) {
				diffuse(effect, 1, effect, scratch);
			}
			var pos = world.grid.pos;
			ScalarField.mult_field(effect, lat, effect);
			return effect;
		}
	} ));

1506828375337-nanyr

Turn this into a function in the AtmosphericModeling namespace, and add a "season" parameter (season==1 for winter in the northern hemisphere and season==-1 for summer in the northern hemisphere)

AtmosphericModeling.surface_air_pressure_land_effect = function(displacement, lat, sealevel, season, effect, scratch) {
	var is_land = ScalarField.gt_scalar(displacement, sealevel);
	BinaryMorphology.erosion(is_land, 2, is_land);
	Float32Raster.FromUint8Raster(is_land, effect);
	var diffuse = ScalarField.diffusion_by_constant;
	var smoothing_iterations =  2;
	for (var i=0; i<smoothing_iterations; ++i) {
		diffuse(effect, 1, effect, scratch);
	}
	var pos = world.grid.pos;
	ScalarField.mult_field(effect, lat, effect);
    ScalarField.mult_scalar(effect, season, effect);
	return effect;
}

OK, now we add the two effects together:

AtmosphericModeling.surface_air_pressure = function(displacement, lat, sealevel, season, axial_tilt, pressure, scratch) {
	scratch = scratch || Float32Raster(lat.grid);
	// "effective latitude" is what the latitude is like weather-wise, when taking axial tilt into account
	var effective_lat = scratch; 
	ScalarField.add_scalar(lat, season*axial_tilt, effective_lat);
	Float32RasterInterpolation.clamp(effective_lat, -Math.PI/2, Math.PI/2, effective_lat);

	AtmosphericModeling.surface_air_pressure_lat_effect(effective_lat, pressure);

	var land_effect = Float32Raster(lat.grid);
	var scratch2 = Float32Raster(lat.grid);
	AtmosphericModeling.surface_air_pressure_land_effect(displacement, effective_lat, sealevel, season, land_effect, scratch2);
	ScalarField.add_scalar_term(pressure, land_effect, 1, pressure);

	return pressure;
}

1506828375337-nanyr 5

Density times the gradient of pressure is acceleration, but we pretend it's velocity.

view.setVectorDisplay( new VectorFieldDisplay( {
		getField: function (world) {
			var lat = Float32SphereRaster.latitude(world.grid.pos.y);
			var scratch = Float32Raster(world.grid);
			var pressure = Float32Raster(world.grid);
			AtmosphericModeling.surface_air_pressure(world.displacement, lat, world.SEALEVEL, 1, Math.PI*23.5/180, pressure, scratch);
			var velocity = ScalarField.gradient(pressure);
			return velocity;
		} 
	} ));

1506828375337-nanyr 6

Now just add the coriolis effect:

AtmosphericModeling.surface_air_velocity_coriolis_effect = function(pos, velocity, angular_speed, effect) {
	effect = effect || VectorRaster(pos.grid);
	VectorField.cross_vector_field	(velocity, pos, 			effect);
	VectorField.mult_scalar 		(effect, 2 * angular_speed, effect);
	VectorField.mult_scalar_field	(effect, pos.y, 			effect);
	return effect;
}
view.setVectorDisplay( new VectorFieldDisplay( {
		getField: function (world) {
			var lat = Float32SphereRaster.latitude(world.grid.pos.y);
			var scratch = Float32Raster(world.grid);
			var pressure = Float32Raster(world.grid);
			AtmosphericModeling.surface_air_pressure(world.displacement, lat, world.SEALEVEL, 1, Math.PI*23.5/180, pressure, scratch);
			var velocity = ScalarField.gradient(pressure);
			var coriolis_effect = AtmosphericModeling.surface_air_velocity_coriolis_effect(world.grid.pos, velocity, ANGULAR_SPEED)
			VectorField.add_vector_field(velocity, coriolis_effect, velocity);
			return velocity;
		} 
	} ));

1506828375337-nanyr 7

You can see how air gets deflected around landmass:

1506828375337-nanyr 9

BTW, you should be able to run all this code in the Developer Tools console. I've also got a "atmo-model" branch out with this stuff added.

@davidson16807

This comment has been minimized.

Show comment
Hide comment
@davidson16807

davidson16807 Nov 12, 2017

Owner

I have these features implemented in prod.

The models aren't strictly rigorous, not scientifically. Pressure and velocity are scaled to the right units, but the pressure model is hand wavey (kind of like existing models for precip and temp), and using gradients to find velocity is not strictly correct. Strictly, density times the gradient of pressure ought to be acceleration, not velocity, but it is suggestive of velocity and it's much more performant than actual fluid simulation.

I need a model that can be run in less than a single frame. It needs to return the average climate over millions of years, not years or decades. There's not a lot of GCMs that meet that criteria. I expect the correct model I'm looking for either relies on statistical mechanics or a steady state assumption. I need to do some more reading.

The next thing I'm going to work on is probably temperature. I have a branch out there to model orbital mechanics, specifically distance and orientation relative to the sun. I'd sample over several points in the day and several points in the year, find equilibrium temperature for a black body at each point, then average them out to get summer/winter temperature.

Owner

davidson16807 commented Nov 12, 2017

I have these features implemented in prod.

The models aren't strictly rigorous, not scientifically. Pressure and velocity are scaled to the right units, but the pressure model is hand wavey (kind of like existing models for precip and temp), and using gradients to find velocity is not strictly correct. Strictly, density times the gradient of pressure ought to be acceleration, not velocity, but it is suggestive of velocity and it's much more performant than actual fluid simulation.

I need a model that can be run in less than a single frame. It needs to return the average climate over millions of years, not years or decades. There's not a lot of GCMs that meet that criteria. I expect the correct model I'm looking for either relies on statistical mechanics or a steady state assumption. I need to do some more reading.

The next thing I'm going to work on is probably temperature. I have a branch out there to model orbital mechanics, specifically distance and orientation relative to the sun. I'd sample over several points in the day and several points in the year, find equilibrium temperature for a black body at each point, then average them out to get summer/winter temperature.

@astrographer

This comment has been minimized.

Show comment
Hide comment
@astrographer

astrographer Nov 12, 2017

astrographer commented Nov 12, 2017

@davidson16807

This comment has been minimized.

Show comment
Hide comment
@davidson16807

davidson16807 Nov 12, 2017

Owner

Pressure does take landmass into account, but the effect only exists during the solstices. Try moving the season slider.

Precip is still your model. I'm unsure how to continue on that front. I know high pressure regions like the "horse" latitudes should be drier. We could also advect wind velocity to determine where water evaporates from the ocean or where rain shadows exist. I don't have a good mechanistic model for evaporation or precipitation, but that's what I want as well.

Erosion is already flow based, in a way. It's sort of related to the gradient of elevation. I spent a lot of time thinking about ways to simulate river valleys this year though and it will take a different approach to do so. A problem I foresee will be grid resolution. It's already working close to the limits to what we can do in a single frame. We can speed up diffused pressure fields for the mantle and atmosphere by switching to a coarser grid, but we can't take that approach for river valleys.

EDIT: I misunderstood what you meant by flow based erosion. Yeah, any improvements we make to precipitation can be fed immediately into our erosion model.

Owner

davidson16807 commented Nov 12, 2017

Pressure does take landmass into account, but the effect only exists during the solstices. Try moving the season slider.

Precip is still your model. I'm unsure how to continue on that front. I know high pressure regions like the "horse" latitudes should be drier. We could also advect wind velocity to determine where water evaporates from the ocean or where rain shadows exist. I don't have a good mechanistic model for evaporation or precipitation, but that's what I want as well.

Erosion is already flow based, in a way. It's sort of related to the gradient of elevation. I spent a lot of time thinking about ways to simulate river valleys this year though and it will take a different approach to do so. A problem I foresee will be grid resolution. It's already working close to the limits to what we can do in a single frame. We can speed up diffused pressure fields for the mantle and atmosphere by switching to a coarser grid, but we can't take that approach for river valleys.

EDIT: I misunderstood what you meant by flow based erosion. Yeah, any improvements we make to precipitation can be fed immediately into our erosion model.

@astrographer

This comment has been minimized.

Show comment
Hide comment
@astrographer

astrographer Nov 12, 2017

astrographer commented Nov 12, 2017

@davidson16807

This comment has been minimized.

Show comment
Hide comment
@davidson16807

davidson16807 Nov 12, 2017

Owner

It doesn't look like you have an arctic/antarctic high pressure zone.

Yeah, how should I fix that? Right now it's pressure = -cos(5*lat). Should it be 6*lat as with precip?

What it really is is an insolation model pretending to be a temperature model, much as pressure gradient is acceleration masquerading as velocity.

I think there's two things that can help with that in the not-too-distant future: 1.) black body equilibrium temperature (mentioned earlier, this is what I'm working on today), and 2.) advection.

If you are familiar with octave/matlab code, I can send you a decent-ish 1d+time insolation model and try to work up from that. May or may not be helpful.

Go for it.

I might look into a "fetch"-based model. Trace winds back to the nearest patch of ocean and base precipitation on inverse distance traversed.

Sounds like advection in essence. For each parcel of air, back-calculate previous position at some timestep. Sample winds from that position and use it to back-calculate further. Repeat as many times as you can afford in a single frame (it might only be once). However many times you choose, it traces a path. Follow that path forward and calculate evaporation and precipitation along the way in some time sensitive manner that's consistent with your back-calculation timestep (you could possibly do this going backward, but it's probably easier to think about if it goes forward). The back-calculation timestep is completely decoupled from the tectonics timestep that occurs every frame. It's shorter by many orders of magnitude, and it should be proportionate to the time it takes for global average windspeed to traverse one grid cell.

Regretfully, I think that should sleep on the back burner.

yeah, totally agree

Owner

davidson16807 commented Nov 12, 2017

It doesn't look like you have an arctic/antarctic high pressure zone.

Yeah, how should I fix that? Right now it's pressure = -cos(5*lat). Should it be 6*lat as with precip?

What it really is is an insolation model pretending to be a temperature model, much as pressure gradient is acceleration masquerading as velocity.

I think there's two things that can help with that in the not-too-distant future: 1.) black body equilibrium temperature (mentioned earlier, this is what I'm working on today), and 2.) advection.

If you are familiar with octave/matlab code, I can send you a decent-ish 1d+time insolation model and try to work up from that. May or may not be helpful.

Go for it.

I might look into a "fetch"-based model. Trace winds back to the nearest patch of ocean and base precipitation on inverse distance traversed.

Sounds like advection in essence. For each parcel of air, back-calculate previous position at some timestep. Sample winds from that position and use it to back-calculate further. Repeat as many times as you can afford in a single frame (it might only be once). However many times you choose, it traces a path. Follow that path forward and calculate evaporation and precipitation along the way in some time sensitive manner that's consistent with your back-calculation timestep (you could possibly do this going backward, but it's probably easier to think about if it goes forward). The back-calculation timestep is completely decoupled from the tectonics timestep that occurs every frame. It's shorter by many orders of magnitude, and it should be proportionate to the time it takes for global average windspeed to traverse one grid cell.

Regretfully, I think that should sleep on the back burner.

yeah, totally agree

@davidson16807

This comment has been minimized.

Show comment
Hide comment
@davidson16807

davidson16807 Nov 15, 2017

Owner

I've been looking at potential evapotranspiration models.

Priestly-Taylor sounds like it's the simplest but requires an empirical factor that varies by location, so that's probably out for the long run.

Penman-Monteith is pretty well regarded and we may be able to supply it with all the parameters we need, but it's probably most complex and it involves things like stomatal conductance that would require assumptions about the planet's plant life.

Penman is simple and only seems to require parameters that we either already have or could potentially represent in the model.

Penman sounds preferable for the long run, but we might get away with Priestly-Taylor for a first pass.

Owner

davidson16807 commented Nov 15, 2017

I've been looking at potential evapotranspiration models.

Priestly-Taylor sounds like it's the simplest but requires an empirical factor that varies by location, so that's probably out for the long run.

Penman-Monteith is pretty well regarded and we may be able to supply it with all the parameters we need, but it's probably most complex and it involves things like stomatal conductance that would require assumptions about the planet's plant life.

Penman is simple and only seems to require parameters that we either already have or could potentially represent in the model.

Penman sounds preferable for the long run, but we might get away with Priestly-Taylor for a first pass.

@davidson16807

This comment has been minimized.

Show comment
Hide comment
@davidson16807

davidson16807 Nov 16, 2017

Owner

Ritchie1998.pdf <- very easy to read description of priestly-taylor model
Penman1948.pdf <- original description of penman model

Owner

davidson16807 commented Nov 16, 2017

Ritchie1998.pdf <- very easy to read description of priestly-taylor model
Penman1948.pdf <- original description of penman model

@davidson16807

This comment has been minimized.

Show comment
Hide comment
@davidson16807

davidson16807 Nov 29, 2017

Owner

My first pass at a mechanistic temperature model was to calculate black body equilibrium temp for each grid cell given that cell's insolation:

https://www.desmos.com/calculator/uhq3bzb5du

Obviously, this doesn't work because the poles reach absolute zero, but I figure it could be a starting point.

Next, I took black body equilibrium and tried to simulate some advection/diffusion.

https://www.desmos.com/calculator/goiuvmxstl

This resulted in an unrealistic graph with jagged peaks and valleys. This is because I started with black body equilibrium, which features a very sharp drop off where the poles reach absolute zero. The jagged drop off propogrates from advection/diffusion. Advection and diffusion need to operate on temperature fields where advection and diffusion have already been taken into consideration, so its sort of a a chicken-egg problem.

The best solution I found was to take a moving average of black body equilibrium.

https://www.desmos.com/calculator/mvhagcn3tt

Except now it's hand waving. I don't have any way to tie the parameters of the moving average model to the parameters for a sensible advection/diffusion model. I have no idea what the size of the moving average window ought to be. I can try to match it to fit observed max/min temp for earth, but what's to say that would apply to other planets?

More to come...

Owner

davidson16807 commented Nov 29, 2017

My first pass at a mechanistic temperature model was to calculate black body equilibrium temp for each grid cell given that cell's insolation:

https://www.desmos.com/calculator/uhq3bzb5du

Obviously, this doesn't work because the poles reach absolute zero, but I figure it could be a starting point.

Next, I took black body equilibrium and tried to simulate some advection/diffusion.

https://www.desmos.com/calculator/goiuvmxstl

This resulted in an unrealistic graph with jagged peaks and valleys. This is because I started with black body equilibrium, which features a very sharp drop off where the poles reach absolute zero. The jagged drop off propogrates from advection/diffusion. Advection and diffusion need to operate on temperature fields where advection and diffusion have already been taken into consideration, so its sort of a a chicken-egg problem.

The best solution I found was to take a moving average of black body equilibrium.

https://www.desmos.com/calculator/mvhagcn3tt

Except now it's hand waving. I don't have any way to tie the parameters of the moving average model to the parameters for a sensible advection/diffusion model. I have no idea what the size of the moving average window ought to be. I can try to match it to fit observed max/min temp for earth, but what's to say that would apply to other planets?

More to come...

@davidson16807 davidson16807 changed the title from Precip/Evap Model Research to Climate Model Research Nov 29, 2017

@davidson16807

This comment has been minimized.

Show comment
Hide comment
@davidson16807

davidson16807 Nov 29, 2017

Owner

Next I tried a simpler model. It couldn't get any simpler: two parcels of air, one at the poles, the other at the equator. It looks a lot like two black bodies in equilibrium, except now there's an additional energy flux that represents how convection moves heat away from the equator and towards poles.

text4247

I played with it for a while. I found I could predict temperature given wind velocity, but then I'm still left with finding velocity. I found I couldn't get anywhere without introducing additional constraints.

https://www.desmos.com/calculator/oxcyyjdmyw

The diagram above was so simple it reminded me of one of those diagrams for a heat engine, and it turned out the analogy is pretty common knowledge. The atmosphere is a kind of heat engine that transforms heat into mechanical work in the form of wind. You could use this fact to predict the maximum theoretical efficiency of the heat engine, and that would give you an upper bound for how strong the wind blows, but that's not quite what I want.

At some point while researching heat engines, I stumbled on this paper:

Lorenz2001.pdf

Wouldn't you know, that's the same diagram I drew!

So here's what the paper is saying: just like any other heat engine, the atmosphere produces entropy. Entropy is the ratio between output energy and temperature. In this case, "output energy" is convective energy flux, and temperature is whatever occurs in a given parcel of air. The difference in entropy between the two air parcels should be a good indicator for the entropy produced.

Now here's the kicker: for a complex system like an atmosphere, there's a lot of degrees of freedom, and whenever you have a lot of degrees of freedom, its postulated that the system tends to maximize the entropy produced.

So if you find a W that maximizes entropy production, then that's probably the W you see in real life. If you know what W is, then you also know what max/min temp is. If max/min temp is known, you can scale your temperature model to the correct values. You could either start with black body equilibrium temp and scale it to correct values (might still have jaggies), or better yet, you can parameterize a moving average diffusion/advection model given known min/max temp.

It's still a form of hand waving, but now the hand waving is peer reviewed!

MEP models are freaky accurate. The paper cited above accurately predicts temperature ranges for Titan and Mars. By comparison, the next best model we have for Titan is off by several orders of magnitude. It's also been used to predict temperatures for Earth
and an extrasolar hot jupiter. The model should be applicable for any terrestrial planet with atmosphere, but it falls apart when dealing with ordinary gas giants or any circumstance where internal heating is non-negligible. It's generally not certain what cut offs exist where the model loses its applicability. Planetary rotation should have a major effect on convective energy flux, and there could exist a cutoff concerning it, but so far nothing has been found to suggest one exists (much to the consternation of climatologists).

Why does MEP work? The best explanation found so far uses statistical mechanics. A highly unconstrained system simply has a higher probability of assuming a state where entropy production is maximized. Dewar (2005) explains why this must necessarily be the case, but I haven't read it yet.

Here's what I implemented using an MEP model:

https://www.desmos.com/calculator/n4bwzkwsx3

The vertical line indicates the convective energy flux at which MEP occurs. The intercept between it and the other lines indicate the max and min temp of earth. Curiously, it only gives sensible max/min values for earth when insolation is doubled. I'm not sure why, but I think its a problem with my implementation. It might have something to do with representing only a single hemisphere. I think this was brought up at some point during my reading. I'll have to go back and check.

Owner

davidson16807 commented Nov 29, 2017

Next I tried a simpler model. It couldn't get any simpler: two parcels of air, one at the poles, the other at the equator. It looks a lot like two black bodies in equilibrium, except now there's an additional energy flux that represents how convection moves heat away from the equator and towards poles.

text4247

I played with it for a while. I found I could predict temperature given wind velocity, but then I'm still left with finding velocity. I found I couldn't get anywhere without introducing additional constraints.

https://www.desmos.com/calculator/oxcyyjdmyw

The diagram above was so simple it reminded me of one of those diagrams for a heat engine, and it turned out the analogy is pretty common knowledge. The atmosphere is a kind of heat engine that transforms heat into mechanical work in the form of wind. You could use this fact to predict the maximum theoretical efficiency of the heat engine, and that would give you an upper bound for how strong the wind blows, but that's not quite what I want.

At some point while researching heat engines, I stumbled on this paper:

Lorenz2001.pdf

Wouldn't you know, that's the same diagram I drew!

So here's what the paper is saying: just like any other heat engine, the atmosphere produces entropy. Entropy is the ratio between output energy and temperature. In this case, "output energy" is convective energy flux, and temperature is whatever occurs in a given parcel of air. The difference in entropy between the two air parcels should be a good indicator for the entropy produced.

Now here's the kicker: for a complex system like an atmosphere, there's a lot of degrees of freedom, and whenever you have a lot of degrees of freedom, its postulated that the system tends to maximize the entropy produced.

So if you find a W that maximizes entropy production, then that's probably the W you see in real life. If you know what W is, then you also know what max/min temp is. If max/min temp is known, you can scale your temperature model to the correct values. You could either start with black body equilibrium temp and scale it to correct values (might still have jaggies), or better yet, you can parameterize a moving average diffusion/advection model given known min/max temp.

It's still a form of hand waving, but now the hand waving is peer reviewed!

MEP models are freaky accurate. The paper cited above accurately predicts temperature ranges for Titan and Mars. By comparison, the next best model we have for Titan is off by several orders of magnitude. It's also been used to predict temperatures for Earth
and an extrasolar hot jupiter. The model should be applicable for any terrestrial planet with atmosphere, but it falls apart when dealing with ordinary gas giants or any circumstance where internal heating is non-negligible. It's generally not certain what cut offs exist where the model loses its applicability. Planetary rotation should have a major effect on convective energy flux, and there could exist a cutoff concerning it, but so far nothing has been found to suggest one exists (much to the consternation of climatologists).

Why does MEP work? The best explanation found so far uses statistical mechanics. A highly unconstrained system simply has a higher probability of assuming a state where entropy production is maximized. Dewar (2005) explains why this must necessarily be the case, but I haven't read it yet.

Here's what I implemented using an MEP model:

https://www.desmos.com/calculator/n4bwzkwsx3

The vertical line indicates the convective energy flux at which MEP occurs. The intercept between it and the other lines indicate the max and min temp of earth. Curiously, it only gives sensible max/min values for earth when insolation is doubled. I'm not sure why, but I think its a problem with my implementation. It might have something to do with representing only a single hemisphere. I think this was brought up at some point during my reading. I'll have to go back and check.

@astrographer

This comment has been minimized.

Show comment
Hide comment
@astrographer

astrographer Nov 30, 2017

So, this
temp_gradient3.m.txt is a little matlab/octave script I modified from an original found here, http://www.gps.caltech.edu/classes/ese148b/files/temp_gradient.m

As you have already observed, the poles experience zero insolation at least some time during the year leading to absolute zero black body temperatures. This problem was disguised in the original script by hardwiring a 23º axial obliquity into the model and only returning the annual average temperatures. As soon as you access the daily temperatures or set the obliquity to zero, the problem becomes obvious. I've had a good deal of trouble figuring out how to do anything like a diffusion in Matlab. To start with, at least, I'd go with something relatively simple, like a convolution kernel.

For best results, I'd do some kind of diffusion both in time and in space. The spatial diffusion, of course would represent movement of warm air over the globe, whereas, the much small temporal diffusion would represent storage effects. Making it two way in time would be a bit odd…

The next step would be to increase the grid to three dimensions(latitude, longitude and day). I'm pretty sure Matlab can handle this, but again my ignorance is showing badly. Or very well…
At this point, we could get interesting land/sea variations by giving the sea a higher diffusion parameter in time(hopefully only with the previous days) to represent storage and a higher diffusion parameter in space to represent movement of heat in ocean currents(another place that would benefit from fluid mechanical models of some sort).

Your last post started out pretty understandable to me and bounded straight over my head from there. I've been pretty much incapable of figuring out temperatures from your MEP model.

Your moving average model had the best results, so I think you're onto something there. While it's not rigorous as is, I wonder if you've tried using that moving average model as an initial input to a more rigorous advection/diffusion model. You might at least avoid the jagged chicken-egg problem. While not precisely rigorous, the blended model is almost certainly closer to a valid solution than the straight up black body solution. At least near the poles.

One problem I have with sanity checking your moving average is that it displaces the minimum temperature to the right of the pole. Part of the solution, I think is to force the output of all latitudes outside of the range -90<=latitude<=90 to zero. The other, and more important part will be to make sure that the point of interest is at the center of the window. I have the feeling that at present it's simply averaging with the four cells to the right(perhaps weighted, I'm not sure). If the moving average is centered on the cell of interest(average of itself, the two cells to the right and the two cells to the left), the result should have a minimum centered on the pole, even with the repeat.

I figure if you repeated the advection/diffusion procedure several times, the jagged artifacts would work themselves out, but the procedure would be slow or unsatisfactory. I think starting with a reasonably massaged initial state would give best results. Although there is no way that storage(even in water) would damp changes over million year time steps, it would probably pay off, for our purposes, to use final temperatures from the previous time step to initialize the model. That might however fail the sanity check if there are large differences between land and sea temperatures being propagated. Perhaps re-initializing from a moving average black body temperature with each time step might be preferable.

Finally, in response to your question about precipitation and pressures…
"Yeah, how should I fix that? Right now it's pressure = -cos(5lat). Should it be 6lat as with precip?"
I think the answer would be yes. That precipitation model, such as it is, basically reflects those pressure bands. Which brings up another point. You should be able to use those pressures as a first approximation to precipitation as it varies over the year. Second approximation would be to have the total amplitude of the rainfall fall off as you approach the poles. A third approximation would be to have precipitation fall off as you move inland. Optimally, that inland falloff would be relegated to an advection/diffusion model along wind paths and weighted by wind strengths. Besides a summer/winter slider for precipitation, I'd like to see an annual average option.

One very nice thing, once we can generate winter and summer extremes for temperature and precipitation, we can create better biome maps. Still better if we can model the flow of water over the land surface. If you have even a simple model of upland moisture concentrations for erosion already, then you can probably adapt that to generating something like that. It'd be awesome to have Egyptian or Mesopotamian areas of increased fertility winding through drier areas fed by moister uplands.

astrographer commented Nov 30, 2017

So, this
temp_gradient3.m.txt is a little matlab/octave script I modified from an original found here, http://www.gps.caltech.edu/classes/ese148b/files/temp_gradient.m

As you have already observed, the poles experience zero insolation at least some time during the year leading to absolute zero black body temperatures. This problem was disguised in the original script by hardwiring a 23º axial obliquity into the model and only returning the annual average temperatures. As soon as you access the daily temperatures or set the obliquity to zero, the problem becomes obvious. I've had a good deal of trouble figuring out how to do anything like a diffusion in Matlab. To start with, at least, I'd go with something relatively simple, like a convolution kernel.

For best results, I'd do some kind of diffusion both in time and in space. The spatial diffusion, of course would represent movement of warm air over the globe, whereas, the much small temporal diffusion would represent storage effects. Making it two way in time would be a bit odd…

The next step would be to increase the grid to three dimensions(latitude, longitude and day). I'm pretty sure Matlab can handle this, but again my ignorance is showing badly. Or very well…
At this point, we could get interesting land/sea variations by giving the sea a higher diffusion parameter in time(hopefully only with the previous days) to represent storage and a higher diffusion parameter in space to represent movement of heat in ocean currents(another place that would benefit from fluid mechanical models of some sort).

Your last post started out pretty understandable to me and bounded straight over my head from there. I've been pretty much incapable of figuring out temperatures from your MEP model.

Your moving average model had the best results, so I think you're onto something there. While it's not rigorous as is, I wonder if you've tried using that moving average model as an initial input to a more rigorous advection/diffusion model. You might at least avoid the jagged chicken-egg problem. While not precisely rigorous, the blended model is almost certainly closer to a valid solution than the straight up black body solution. At least near the poles.

One problem I have with sanity checking your moving average is that it displaces the minimum temperature to the right of the pole. Part of the solution, I think is to force the output of all latitudes outside of the range -90<=latitude<=90 to zero. The other, and more important part will be to make sure that the point of interest is at the center of the window. I have the feeling that at present it's simply averaging with the four cells to the right(perhaps weighted, I'm not sure). If the moving average is centered on the cell of interest(average of itself, the two cells to the right and the two cells to the left), the result should have a minimum centered on the pole, even with the repeat.

I figure if you repeated the advection/diffusion procedure several times, the jagged artifacts would work themselves out, but the procedure would be slow or unsatisfactory. I think starting with a reasonably massaged initial state would give best results. Although there is no way that storage(even in water) would damp changes over million year time steps, it would probably pay off, for our purposes, to use final temperatures from the previous time step to initialize the model. That might however fail the sanity check if there are large differences between land and sea temperatures being propagated. Perhaps re-initializing from a moving average black body temperature with each time step might be preferable.

Finally, in response to your question about precipitation and pressures…
"Yeah, how should I fix that? Right now it's pressure = -cos(5lat). Should it be 6lat as with precip?"
I think the answer would be yes. That precipitation model, such as it is, basically reflects those pressure bands. Which brings up another point. You should be able to use those pressures as a first approximation to precipitation as it varies over the year. Second approximation would be to have the total amplitude of the rainfall fall off as you approach the poles. A third approximation would be to have precipitation fall off as you move inland. Optimally, that inland falloff would be relegated to an advection/diffusion model along wind paths and weighted by wind strengths. Besides a summer/winter slider for precipitation, I'd like to see an annual average option.

One very nice thing, once we can generate winter and summer extremes for temperature and precipitation, we can create better biome maps. Still better if we can model the flow of water over the land surface. If you have even a simple model of upland moisture concentrations for erosion already, then you can probably adapt that to generating something like that. It'd be awesome to have Egyptian or Mesopotamian areas of increased fertility winding through drier areas fed by moister uplands.

@davidson16807

This comment has been minimized.

Show comment
Hide comment
@davidson16807

davidson16807 Nov 30, 2017

Owner

One problem I have with sanity checking your moving average is that it displaces the minimum temperature to the right of the pole

Good point. I started out modeling air as a particle moving through the atmosphere before switching to modeling a stationary air parcel. The min temp displacement made a lot more sense in the original model, where the x axis was time.

I have the feeling that at present it's simply averaging with the four cells to the right(perhaps weighted, I'm not sure)

actually, that's also an artifact of the earlier model. I really didn't do a good job transitioning that!

Perhaps re-initializing from a moving average black body temperature with each time step might be preferable.

This is the approach I'd like to start with. If performance becomes a concern we can always rely on temperatures from the last time step to give a better initial guess.

I do think the proper solution going forward is to initialize with black body temp and then run some sort of space/time averaging using parameters from the MEP model. Bonus points if the space/time averaging also happens to be some sort of rigorous physics simulation.

So we want to do two types of averaging: "spatial" and "temporal". I'm pretty sure the "spatial" part can be made rigorous. I'm thinking about using the convection/diffusion equation, but we might only need the diffusion equation. From the MEP model, we can derive max temp, min temp, wind velocity, and what's known as the "meridional diffusion coefficient", D. Taking W from the diagram above, the Lorenz paper says that W=2DΔT. I think the "D" coefficient from the Lorenz paper could be the same "D" coefficient from the wiki article above. I think it's kind of like you're modeling convection as if it were conduction. Even if it isn't, we can simply run the diffusion equation a few times with some arbitrary parameters until max and min temp start to resemble the values predicted by MEP. The thought also occurred to me that the diffusion equation uses the laplacian, and the laplacian operator is sort of analogous to a "moving average" or "convolutional kernel", so that could explain why those methods worked well.

As for the temporal averaging... the only way I see it being made rigorous is to run a simulation over many time steps, but I can guarantee that's going to be too slow.

I know for certain we'll need temporal averaging, because from my own experiments I know summer time black body equilibrium temp is higher at the poles than it is at the equator. I don't think the temperatures in Antarctica ever get hotter than the temperatures at the equator, at least not on average at the same time of the year. The MEP model doesn't solve this problem, either, because it only tells what the max/min temps are, not where they're located.

You mention that temperature model worked around its problem by averaging summer and winter temps. I think something similar could be done for this model, but instead of just sampling two points across time, you could sample from many points and then have some sort of weighted average between them. Let's say you sample black body temperature T from times t=1,2,3. Your weighted average for t=3 could be (3T(3) + 2T(2) + 1T(1)) / 6, or something similar. I'd just like to see if something like this could have some more physical justification. That way, I could put more faith in the results and spend less time parameterizing. I could also potentially generalize a physically justified model to examine a planet across other cycles: day/night cycles, annual cycles, binary star system cycles, milankovitch cycles, etc.

Besides a summer/winter slider for precipitation, I'd like to see an annual average option

It's probably going to wind up being equivalent to precip during the equinoxes.

One very nice thing, once we can generate winter and summer extremes for temperature and precipitation, we can create better biome maps. Still better if we can model the flow of water over the land surface

😁

Owner

davidson16807 commented Nov 30, 2017

One problem I have with sanity checking your moving average is that it displaces the minimum temperature to the right of the pole

Good point. I started out modeling air as a particle moving through the atmosphere before switching to modeling a stationary air parcel. The min temp displacement made a lot more sense in the original model, where the x axis was time.

I have the feeling that at present it's simply averaging with the four cells to the right(perhaps weighted, I'm not sure)

actually, that's also an artifact of the earlier model. I really didn't do a good job transitioning that!

Perhaps re-initializing from a moving average black body temperature with each time step might be preferable.

This is the approach I'd like to start with. If performance becomes a concern we can always rely on temperatures from the last time step to give a better initial guess.

I do think the proper solution going forward is to initialize with black body temp and then run some sort of space/time averaging using parameters from the MEP model. Bonus points if the space/time averaging also happens to be some sort of rigorous physics simulation.

So we want to do two types of averaging: "spatial" and "temporal". I'm pretty sure the "spatial" part can be made rigorous. I'm thinking about using the convection/diffusion equation, but we might only need the diffusion equation. From the MEP model, we can derive max temp, min temp, wind velocity, and what's known as the "meridional diffusion coefficient", D. Taking W from the diagram above, the Lorenz paper says that W=2DΔT. I think the "D" coefficient from the Lorenz paper could be the same "D" coefficient from the wiki article above. I think it's kind of like you're modeling convection as if it were conduction. Even if it isn't, we can simply run the diffusion equation a few times with some arbitrary parameters until max and min temp start to resemble the values predicted by MEP. The thought also occurred to me that the diffusion equation uses the laplacian, and the laplacian operator is sort of analogous to a "moving average" or "convolutional kernel", so that could explain why those methods worked well.

As for the temporal averaging... the only way I see it being made rigorous is to run a simulation over many time steps, but I can guarantee that's going to be too slow.

I know for certain we'll need temporal averaging, because from my own experiments I know summer time black body equilibrium temp is higher at the poles than it is at the equator. I don't think the temperatures in Antarctica ever get hotter than the temperatures at the equator, at least not on average at the same time of the year. The MEP model doesn't solve this problem, either, because it only tells what the max/min temps are, not where they're located.

You mention that temperature model worked around its problem by averaging summer and winter temps. I think something similar could be done for this model, but instead of just sampling two points across time, you could sample from many points and then have some sort of weighted average between them. Let's say you sample black body temperature T from times t=1,2,3. Your weighted average for t=3 could be (3T(3) + 2T(2) + 1T(1)) / 6, or something similar. I'd just like to see if something like this could have some more physical justification. That way, I could put more faith in the results and spend less time parameterizing. I could also potentially generalize a physically justified model to examine a planet across other cycles: day/night cycles, annual cycles, binary star system cycles, milankovitch cycles, etc.

Besides a summer/winter slider for precipitation, I'd like to see an annual average option

It's probably going to wind up being equivalent to precip during the equinoxes.

One very nice thing, once we can generate winter and summer extremes for temperature and precipitation, we can create better biome maps. Still better if we can model the flow of water over the land surface

😁

@davidson16807

This comment has been minimized.

Show comment
Hide comment
@davidson16807

davidson16807 Dec 1, 2017

Owner

I read the explanation by Dewar. It feels like one of those ideas you could explain in a sentence, but every time I try it winds up being a paragraph. Here goes:

Let's say I want to guess the position of every molecule in some random parcel of air. What's more likely: the molecules are all scrunched up in one corner, or the molecules are evenly distributed? You're probably going to say they're more likely to be evenly distributed, but why? Its because there are way more states that describe molecules being equally distributed than there are states that describe molecules being scrunched in one corner. We can phrase this several different ways:

  • A scenario is more likely when it has more possible interpretations
  • A scenario is more likely when it tells us less about the true state of things.
  • A scenario is more likely when it holds less "information".
  • A scenario is more likely when it has high "entropy"

OK, now let's say I drop some creamer into coffee. What's more likely: the creamer stays scrunched up, or the molecules diffuse in the usual way? Obviously the latter, but why? It's because the latter scenario describes more possible outcomes. The statement has more possible interpretations. It contains less "information." It produces more "entropy". You can see this and start to understand how any scenario that produces more entropy is inherently going to be more likely. Once you're comfortable with this, it follows that the most likely scenario of all is the one maximizes the production of entropy.

This statement has nothing to do with physics. It rests solely on statistics, and it's been proven rigorously. If your MEP model doesn't produce the most likely outcome, it's not because the MEP principle is false. Rather, the problem lies with your implementation: what constraints do you use, how did you define entropy, etc.

If I could cut it down to a sentence, it's this: entropy is a measure for how vague a prediction is, and the most likely prediction is the vaguest one

Owner

davidson16807 commented Dec 1, 2017

I read the explanation by Dewar. It feels like one of those ideas you could explain in a sentence, but every time I try it winds up being a paragraph. Here goes:

Let's say I want to guess the position of every molecule in some random parcel of air. What's more likely: the molecules are all scrunched up in one corner, or the molecules are evenly distributed? You're probably going to say they're more likely to be evenly distributed, but why? Its because there are way more states that describe molecules being equally distributed than there are states that describe molecules being scrunched in one corner. We can phrase this several different ways:

  • A scenario is more likely when it has more possible interpretations
  • A scenario is more likely when it tells us less about the true state of things.
  • A scenario is more likely when it holds less "information".
  • A scenario is more likely when it has high "entropy"

OK, now let's say I drop some creamer into coffee. What's more likely: the creamer stays scrunched up, or the molecules diffuse in the usual way? Obviously the latter, but why? It's because the latter scenario describes more possible outcomes. The statement has more possible interpretations. It contains less "information." It produces more "entropy". You can see this and start to understand how any scenario that produces more entropy is inherently going to be more likely. Once you're comfortable with this, it follows that the most likely scenario of all is the one maximizes the production of entropy.

This statement has nothing to do with physics. It rests solely on statistics, and it's been proven rigorously. If your MEP model doesn't produce the most likely outcome, it's not because the MEP principle is false. Rather, the problem lies with your implementation: what constraints do you use, how did you define entropy, etc.

If I could cut it down to a sentence, it's this: entropy is a measure for how vague a prediction is, and the most likely prediction is the vaguest one

@astrographer

This comment has been minimized.

Show comment
Hide comment
@astrographer

astrographer Dec 1, 2017

astrographer commented Dec 1, 2017

@davidson16807

This comment has been minimized.

Show comment
Hide comment
@davidson16807

davidson16807 Apr 12, 2018

Owner

This and #15 are going to be my top priority now that #18 and #14 are done.

Next step will be to partition out functionality within the World class. There's going to be separate classes for each submodel: lithosphere, hydrosphere, atmosphere, biosphere. The World class is going to be way too convoluted if I properly implement atmospheres without any separation of concern.

Owner

davidson16807 commented Apr 12, 2018

This and #15 are going to be my top priority now that #18 and #14 are done.

Next step will be to partition out functionality within the World class. There's going to be separate classes for each submodel: lithosphere, hydrosphere, atmosphere, biosphere. The World class is going to be way too convoluted if I properly implement atmospheres without any separation of concern.

@davidson16807

This comment has been minimized.

Show comment
Hide comment
@davidson16807

davidson16807 May 2, 2018

Owner

As of commit b09df17, the model is producing temperature estimates in keeping with the Maximum Entropy Production (MEP) principle.

This is big: temperature is now being calculated in a scientifically defensible manner. It's not just some trickery like it was in the past when we rescaled mean annual insolation to earth-like temperature values.

I'm solidly pleased with the results:

  • The model is fairly accurate. Temperatures during the equinoxes are looking pretty good. More work will be needed to get solstice temperatures right, though. More on that later.

  • The model is stupid simple, and requires no other input besides max and min insolation.

  • The model is blazing fast, way better than trying to run a GCM every timestep.

  • The model is very well generalized. It should be able to approximate max and min temperatures for virtually any terrestrial planet with an atmosphere, even where the processes at play are not well understood, as with the Martian CO2 sublimation cycle discussed by Lorenz.

For a model like ours, which requires really fast temperature approximations every timestep on a world that's constantly evolving, this model is definitely the way to go.

Owner

davidson16807 commented May 2, 2018

As of commit b09df17, the model is producing temperature estimates in keeping with the Maximum Entropy Production (MEP) principle.

This is big: temperature is now being calculated in a scientifically defensible manner. It's not just some trickery like it was in the past when we rescaled mean annual insolation to earth-like temperature values.

I'm solidly pleased with the results:

  • The model is fairly accurate. Temperatures during the equinoxes are looking pretty good. More work will be needed to get solstice temperatures right, though. More on that later.

  • The model is stupid simple, and requires no other input besides max and min insolation.

  • The model is blazing fast, way better than trying to run a GCM every timestep.

  • The model is very well generalized. It should be able to approximate max and min temperatures for virtually any terrestrial planet with an atmosphere, even where the processes at play are not well understood, as with the Martian CO2 sublimation cycle discussed by Lorenz.

For a model like ours, which requires really fast temperature approximations every timestep on a world that's constantly evolving, this model is definitely the way to go.

@davidson16807

This comment has been minimized.

Show comment
Hide comment
@davidson16807

davidson16807 May 2, 2018

Owner

I've been reading more on the MEP principle and stumbled on what could be a solution to modeling seasonal temperatures. Crooks 1999 mentions that a macroscopic system doesn't just tend to maximize entropy production during a steady-state. It maximizes entropy produced across an arbitrary time interval.

Owner

davidson16807 commented May 2, 2018

I've been reading more on the MEP principle and stumbled on what could be a solution to modeling seasonal temperatures. Crooks 1999 mentions that a macroscopic system doesn't just tend to maximize entropy production during a steady-state. It maximizes entropy produced across an arbitrary time interval.

@davidson16807

This comment has been minimized.

Show comment
Hide comment
@davidson16807

davidson16807 May 4, 2018

Owner

seasonal temperature demo:

https://jsfiddle.net/16807/dcvtm08p/111/

The dots indicate the temperature at two latitudes, 90° and 45°. Temperature is reported in degrees Celcius.

The arrow indicates the heat flux due to entropy. The heat flux is reported in Watts per square meter.

Note I say the heat flux is due to entropy, not wind. Wind is just one of several processes that contribute to entropy. I would say it's generated by wind, but that would invite confusion as to what's being simulated.

All temperatures are generated strictly from physics simulation using real world parameters. The only parameters fed to the model are heat capacities, insolations, and the greenhouse gas factor.

Owner

davidson16807 commented May 4, 2018

seasonal temperature demo:

https://jsfiddle.net/16807/dcvtm08p/111/

The dots indicate the temperature at two latitudes, 90° and 45°. Temperature is reported in degrees Celcius.

The arrow indicates the heat flux due to entropy. The heat flux is reported in Watts per square meter.

Note I say the heat flux is due to entropy, not wind. Wind is just one of several processes that contribute to entropy. I would say it's generated by wind, but that would invite confusion as to what's being simulated.

All temperatures are generated strictly from physics simulation using real world parameters. The only parameters fed to the model are heat capacities, insolations, and the greenhouse gas factor.

@davidson16807

This comment has been minimized.

Show comment
Hide comment
@davidson16807

davidson16807 May 5, 2018

Owner

I think this is the correct way to model temperatures as a scalar field:

https://jsfiddle.net/16807/rk1nzk5k/7/

I'm not sure its strictly correct, but it does produce good results. Good enough to press forward.

Owner

davidson16807 commented May 5, 2018

I think this is the correct way to model temperatures as a scalar field:

https://jsfiddle.net/16807/rk1nzk5k/7/

I'm not sure its strictly correct, but it does produce good results. Good enough to press forward.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment