No-Atmosphere Climate Model

xholonWorkbook.xml

<?xml version="1.0" encoding="UTF-8"?>
<!--Xholon Workbook http://www.primordion.com/Xholon/wb/ (C) Ken Webb Thu May 03 2012 09:56:03 GMT-0400 (EDT)-->
<XholonWorkbook>

<Notes><![CDATA[
Xholon
------
Title: No-Atmosphere Climate Model
Description: implementation of Robinson, Walter A. (2001) Modeling Dynamic Climate Systems. section 2.1 "No-Atmosphere Climate Model"
Url: http://www.primordion.com/Xholon/wb/xholon
InternalName: 
YoutubeId: 
Keywords: 

My Notes
--------

How to run the climate model implemented on this web page ::
  Click the Run button above.
  Click the Pause button in the overlay to start (unpause) the app.
  Watch the Celsius and Kelvin temperature values, and the 4 flux values, as they are written to the textarea.
  Click Pause again, when the numbers stop changing.
  The final temperature value should be -18.2830417193793 °Celsius.

This page is about the simplest climate model in the Robinson book. Solar energy flows to the surface of the Earth where it is absorbed, gradually changing the temperature to a stable value. There is no atmosphere. The surface consists of a volume of water.

He implemented the model using STELLA, while I am implementing it as a self-contained Xholon Workbook. All of the app-specific code is contained within editors on this page.

I have previously written complete Java-based implementations of this model, using Xholon jvmEdition. That source code is available at SourceForge. There are two slightly different versions, one of which uses the JScience library ::
  http://xholon.cvs.sourceforge.net/viewvc/xholon/xholon/src/org/primordion/user/app/climatechange/mdcs/m2_1mp/
  http://xholon.cvs.sourceforge.net/viewvc/xholon/xholon/src/org/primordion/user/app/climatechange/mdcs/m2_1jscience/
The m2_1mp version can be run from the following page. Click on the launch button under Climate Change models ::
  http://www.primordion.com/Xholon/jnlp/
The Java-based implementations are complete in the sense that they reproduce the results, and the temperature-vs-time graph, found in Robinson's book. This Xholon Workbook version of the model, is based on the earlier m2_1jscience version.

Final stable values in the original STELLA implementation, and in this Xholon Workbook implementation, and in the Java-based implementation, are ::
  Celsius: -18.2830417193793  Kelvin: 254.86695828062068
  SunSpc_sw: 1367 W/m^2
  SpcSrf_sw: 341.75 W/m^2
  SrfSpc_sw: 102.52499999999999 W/m^2
  SrfSpc_lw: 239.22500000000105 W/m^2
Units of the radiative fluxes are Watts per square meter (W/m^2).
The name of each flux includes the sender, the receiver, and the type of radiation. For example, "SunSpc_sw" is the amount of shortwave (sw) radiation that the Sun radiates into Space (the amount of sunlight in Watts per square meter that reaches the distance of the Earth), while "SrfSpc_lw" is the amount of longwave heat (lw) radiation that the Surface of the Earth radiates back into Space.

Note that when this equilibrium is reached, when none of the values change anymore, the incoming and outgoing radiation exactly balance each other ::
  SpcSrf_sw = SrfSpc_sw + SrfSpc_lw = 102.525 + 239.225 = 341.75 W/m^2

Values after 10 years (10 time steps) are ::
  Year: 9
  Celsius: -18.232368712936363  Kelvin: 254.9176312870636
  SunSpc_sw: 1367
  SpcSrf_sw: 341.75
  SrfSpc_sw: 102.52499999999999
  SrfSpc_lw: 239.56431138367597

STELLA, the software used in the Robinson book, is a modeling and simulation tool based on system dynamics. STELLA uses ordinary differential equations (ODE). ::
  http://www.iseesystems.com/softwares/Education/StellaSoftware.aspx
  http://en.wikipedia.org/wiki/System_dynamics

Instead of using ordinary differential equations (ODE), Xholon Workbooks use the simplest form of numerical integration (Euler method). Each time step is executed multiple times, as specified by a JavaScript TimeStepMultiplier parameter. This workbook has a TimeStepMultiplier of 10, which means that each time step is broken into 10 smaller steps. The workbook also uses a JavaScript dt variable, meaning "delta time". This is used as part of the calculation of energy each time step. Scroll down to the act() function in the Surfacebehavior editor. ::
  dt = 1 / TimeStepMultiplier
According to Calculus, the smaller the value of dt, the more accurate the result should be. If TimeStepMultiplier = 1, then dt = 1/1 = 1. If TimeStepMultiplier = 10, then dt = 1/10 = 0.1 . A dt of 0.1 should produce a more accurate value than a dt of 1 .
See ::
  http://en.wikipedia.org/wiki/Euler_method

How to change the value of TimeStepMultiplier and dt ::
  Scroll down to the TheSystembehavior editor.
  Change the value of TimeStepMultiplier from 10 to 1, and rerun the app to see what happens to the temperature and other values each time step.
  The end result is the same, but the initial changes are more abrupt.

How to find a good value for TimeStepMultiplier and dt ::
  Start with a relatively large initial value of TimeStepMultiplier, anything up to 128.
  Try smaller values until you find a value that does not effect the accuracy by too much.
You don't need to know anything about Calculus to do this. Just use trial and error.


Description of what's happening in this model
---------------------------------------------
Every second the Sun releases shortwave (sw) energy in all directions into Space.
A tiny amount of this energy (about 1367.0 Watts per square meter) reaches the circular area in Space where the Earth is located.
Because half of the Surface is always in nighttime, and because the higher latitudes receive less sunlight, the incoming solar energy is divided by 4 to derive a value (about 341.75 Watts per square meter) averaged over 24 hours and all latitudes.
In this model there is no Atmosphere, so all of the solar energy reaches the Surface.
The Surface consists of Water to a depth of 50 meters.
The Surface is given an arbitrary initial temperature of 0°C, and a corresponding amount of energy.
The Surface reflects 30% of the solar energy (albedo = 0.3), and absorbs the other 70%.
The Surface also emits some of its energy as heat (longwave (lw) energy).
Initially the Surface cools relatively rapidly, but stabilizes at around -18°C once the incoming shortwave solar energy balances the outgoing longwave energy.
Surface energy is measured as temperature, and the changing temperature is plotted on a graph.

Basically, the model calculates the actual surface temperature, given an initial guess of 0°C. The graph shows the changing temperature as it moves from the initial guess to the actual value.
]]></Notes>

<script implName="lang:python:inline:"><![CDATA[

]]></script>

<script implName="lang:javascript:inline:"><![CDATA[

]]></script>

<_-.XholonClass>

  <!-- types of domain objects -->
  <TheSystem/>
  
  <SolarSystem/>
  <Sun/>
  <Planet>
    <Earth/>
  </Planet>
  <Space/>
  
  <Surface/> <!-- "the planet is represented by an average square meter of its surface"; ocean -->
  <Water/> <!-- 70% of the Earth's surface is water; for the model we assume 100% of the surface is water -->
  <MixedLayer/> <!-- the top layer of the ocean, where significant mixing of heat/etc occurs -->
  
  <!-- fluxes -->
  <Flux/>
  <Fluxes/>
  
  <!-- constants -->
  <Constant>
    <StefanBoltzmannConstant/>
    <SecondsPerTimeStep/>
  </Constant>
  <Constants/>
  
  <!-- quantities -->
  <Depth superClass="Length"/>
  <Energy/>
  <Albedo superClass="Dimensionless"/>
  <Temperature> <!-- K -->
    <Celsius/> <!-- degrees C -->
  </Temperature>
  <Density superClass="VolumetricDensity"/>
  <SpecificHeat/>
  <HeatCapacity/>
  
</_-.XholonClass>

<xholonClassDetails>
  <Sun xhType="XhtypePureActiveObject">
    <port name="space" connector="#xpointer(ancestor::SolarSystem/Space)"/>
  </Sun>
  <Space xhType="XhtypePureActiveObject">
    <port name="planet" connector="#xpointer(ancestor::SolarSystem/Earth/Surface)"/>
    <port name="solarConstant" connector="#xpointer(ancestor::TheSystem/Fluxes/Flux[1])"/>
    <port name="longwaveRadiation" connector="#xpointer(ancestor::TheSystem/Fluxes/Flux[4])"/>
    <port name="reflectedShortwaveRadiation" connector="#xpointer(ancestor::TheSystem/Fluxes/Flux[3])"/>
  </Space>
  <Surface xhType="XhtypePureActiveObject">
    <port name="space" connector="#xpointer(ancestor::SolarSystem/Space)"/>
    <port name="energy" connector="#xpointer(Energy)"/>
    <port name="albedo" connector="#xpointer(Albedo)"/>
    <port name="temperature" connector="#xpointer(Temperature)"/>
    <port name="water" connector="#xpointer(Water)"/>
    <port name="stefanBoltzmannConstant" connector="#xpointer(ancestor::TheSystem/Constants/StefanBoltzmannConstant)"/>
    <port name="secondsPerTimeStep" connector="#xpointer(ancestor::TheSystem/Constants/SecondsPerTimeStep)"/>
    <port name="solarIn" connector="#xpointer(ancestor::TheSystem/Fluxes/Flux[2])"/>
  </Surface>
  <Water xhType="XhtypePureActiveObject">
    <port name="depth" connector="#xpointer(MixedLayer/Depth)"/>
    <port name="heatCapacity" connector="#xpointer(HeatCapacity)"/>
    <port name="density" connector="#xpointer(Density)"/>
    <port name="specificHeat" connector="#xpointer(SpecificHeat)"/>
  </Water>
  <Celsius xhType="XhtypeBehFgsxxx">
    <port name="kelvin" connector="#xpointer(..)"/>
  </Celsius>
</xholonClassDetails>

<TheSystem>
  
  <Constants>
    <StefanBoltzmannConstant>5.6696E-8</StefanBoltzmannConstant>
    <SecondsPerTimeStep>31536000.0</SecondsPerTimeStep> <!-- SecondsPerYear 365*24*60*60 -->
    <!--<SecondsPerTimeStep units="s">86400.0</SecondsPerTimeStep>--> <!-- SecondsPerDay -->
    <!--<SecondsPerTimeStep units="s">2592000.0</SecondsPerTimeStep>--> <!-- SecondsPerMonth -->
  </Constants>
  
  <SolarSystem>
    <Sun/>
    <Space/>
    <Earth>
      <Surface>
        <Energy>0.0 J</Energy>
        <Albedo>0.3</Albedo> <!-- Earth -->
        <!--<Albedo>0.136</Albedo>--> <!-- Moon -->
        <Temperature>273.15 K
          <Celsius>0.0 ℃</Celsius>
        </Temperature> <!-- initial temperature -->
        <Water>
          <Density>1000.0 kg/m^3</Density>
          <SpecificHeat>4218.0 J/kg/K</SpecificHeat>
          <MixedLayer>
            <Depth>50.0 m</Depth>
          </MixedLayer>
          <HeatCapacity>0.0 J/m^3/K</HeatCapacity>
        </Water>
      </Surface>
    </Earth>
  </SolarSystem>
  
  <Fluxes> <!-- with optional initial values -->
    <Flux roleName="SunSpc_sw">0.0 W/m^2</Flux> <!-- solar constant -->
    <Flux roleName="SpcSrf_sw">341.75 W/m^2</Flux> <!-- solar -->
    <Flux roleName="SrfSpc_sw">0.0 W/m^2</Flux> <!-- reflected -->
    <Flux roleName="SrfSpc_lw">0.0 W/m^2</Flux> <!-- infrared (IR) -->
  </Fluxes>
  
</TheSystem>

<Blockbehavior implName="lang:python:inline:"><![CDATA[

]]></Blockbehavior>

<Blockbehavior implName="lang:javascript:inline:"><![CDATA[

]]></Blockbehavior>

<TheSystembehavior implName="lang:webEditionjs:inline:"><![CDATA[
  function postConfigure() {
    $("textarea#btstrp_console").css("font-size", "12px");
    $("textarea#btstrp_console").attr("rows", "5");
    $("textarea#btstrp_console").attr("cols", "50");
    this.application('setConfigParam', 'TimeStepMultiplier', 10); // 10
    dt = 1 / this.application('getConfigParam', 'TimeStepMultiplier');
  }
  
  function preAct() {
    print("\n\nYear: " + this.application('getTimeStep'));
  }
]]></TheSystembehavior>

<Sunbehavior implName="lang:webEditionjs:inline:"><![CDATA[
  function postConfigure() {
    this.parent().attr("solarConstant", 1367.0);
  }
  
  function preAct() {
    this.bindPorts(this.parent());
    this.space.sendMessage("SIG_SW", this.parent().attr("solarConstant"), this.parent());
  }
]]></Sunbehavior>

<Spacebehavior implName="lang:webEditionjs:inline:"><![CDATA[
  function processReceivedMessage(msg) {
    this.bindPorts(this.parent());
    switch (msg.signal) {
    case "SIG_SW":
      this.solarConstant.attr("val", msg.data);
      
      // this is how divisor is calculated
      // Average the solar constant over one day, and over all latitudes
      //var radiusOfUnitCircle = 1.0; 
      //var areaOfCircle = Math.PI * Math.pow(radiusOfUnitCircle, 2); // cross-section of Earth as it appears to sunlight, is a circle
      //var areaOfSphere = 4.0 * Math.PI * Math.pow(radiusOfUnitCircle, 2); // Earth is a sphere
      //var divisor = areaOfSphere / areaOfCircle; // result: 4.0
      
      var divisor = 4.0;
      this.planet.sendMessage("SIG_SW", msg.data / divisor, this);
      break;
    case "SIG_LW":
      this.longwaveRadiation.attr("val", msg.data);
      break;
    case "SIG_SW_REFLECTED":
      this.reflectedShortwaveRadiation.attr("val", msg.data);
      break;
    default: break;
    }
  }
]]></Spacebehavior>

<Waterbehavior implName="lang:webEditionjs:inline:"><![CDATA[
  function postConfigure() {
    this.bindPorts(this.parent());
    this.heatCapacity.attr("val", this.density.attr("val") * this.depth.attr("val") * this.specificHeat.attr("val"));
  }
]]></Waterbehavior>

<Surfacebehavior implName="lang:webEditionjs:inline:"><![CDATA[
  function postConfigure() {
    var $this = this;
    this.bindPorts(this.parent());
    this.energy.attr("val", this.temperature.attr("val") * this.water.children(".HeatCapacity").attr("val"));
    this.parent().attr("solar", this.solarIn.attr("val") * (1.0 - this.albedo.attr("val")) * this.secondsPerTimeStep.attr("val"));
    this.parent().attr("infrared", this.stefanBoltzmannConstant.attr("val") * Math.pow(this.temperature.attr("val"), 4) * this.secondsPerTimeStep.attr("val"));
  }
  
  function act() {
    this.bindPorts(this.parent());
    this.energy.attr("val", this.energy.attr("val") - 0 + (this.parent().attr("solar") - this.parent().attr("infrared")) * dt);
    this.temperature.attr("val", this.energy.attr("val") / this.water.children(".HeatCapacity").attr("val"));
    this.parent().attr("solar", this.solarIn.attr("val") * (1.0 - this.albedo.attr("val")) * this.secondsPerTimeStep.attr("val"));
    this.parent().attr("infrared", this.stefanBoltzmannConstant.attr("val") * Math.pow(this.temperature.attr("val"), 4) * this.secondsPerTimeStep.attr("val"));
  }
  
  function postAct() {
    this.bindPorts(this.parent());
    this.space.sendMessage("SIG_SW_REFLECTED", this.solarIn.attr("val") * this.albedo.attr("val"), this);
    this.space.sendMessage("SIG_LW", this.parent().attr("infrared") / this.secondsPerTimeStep.attr("val"), this);
  }
  
  function processReceivedMessage(msg) {
    this.bindPorts(this.parent());
    switch (msg.signal) {
    case "SIG_SW":
      this.solarIn.attr("val", msg.data);
      break;
    default: break;
    }
  }
]]></Surfacebehavior>

<Fluxbehavior implName="lang:webEditionjs:inline:"><![CDATA[
  function postAct() {
    print("\n" + this.parent().attr("rolename") + ": " + this.parent().attr("val") + " " + this.parent().attr("unit"));
  }
]]></Fluxbehavior>

<Celsiusbehavior implName="lang:webEditionjs:inline:"><![CDATA[
  function postConfigure() {
    this.bindPorts(this.parent());
    this.parent().attr("val", this.kelvin.attr("val") - 273.15);
  }
  
  function act() {
    // draw a graph
    var oldCelsius = this.parent().attr("val");
    var newCelsius = this.parent().parent().attr("val") - 273.15;
    this.parent().attr("val", newCelsius);
    var dCelsius = oldCelsius - newCelsius;
    var path = this.parent().data("path");
    if (!path) {
      var svg = $($('div#mySVGDiv > object')[0].contentDocument.getElementsByTagNameNS(svgns, 'svg')[0]);
      path = svg.children("g").children("g#Graph").children("path");
      this.parent().data("path", path);
      this.parent().data("d", path.attr("d"));
    }
    var d = this.parent().data("d") + "l4 " + (dCelsius * 10);
    this.parent().data("d", d);
    path.attr("d", d);
  }
  
  function postAct() {
    this.bindPorts(this.parent());
    var kelvin = this.kelvin.attr("val");
    var celsius = kelvin - 273.15;
    print("\nCelsius: " + celsius + "  Kelvin: " + kelvin);
  }
]]></Celsiusbehavior>

<Blockbehavior implName="lang:bsh:inline:"><![CDATA[

]]></Blockbehavior>

<Blockbehavior implName="lang:jruby:inline:"><![CDATA[

]]></Blockbehavior>

<Blockbehavior implName="lang:groovy:inline:"><![CDATA[

]]></Blockbehavior>

<SvgClient><Attribute_String roleName="svgUri"><![CDATA[data:image/svg+xml,
<svg width="400" height="250" xmlns="http://www.w3.org/2000/svg">
 <g>
  <g id="Space">
    <title>No-Atmosphere Climate Model - Space</title>
    <rect id="TheSystem/SolarSystem/Space" fill="black" height="50" width="400"/>
    <g>
      <title>Sun</title>
      <circle id="TheSystem/SolarSystem/Sun" fill="yellow" stroke="black" cx="15" cy="10" r="5"/>
    </g>
    <g>
      <title>Surface with Water</title>
      <circle id="TheSystem/SolarSystem/Earth/Surface/Water" fill="#3CB371" cx="370" cy="25" r="22"/>
    </g>
    <g>
      <title>Earth</title>
      <circle id="TheSystem/SolarSystem/Earth" fill="#CD853F" cx="370" cy="25" r="20"/>
    </g>
  </g>
  <g id="Graph">
    <title>Graph of Earth Surface Temperature vs Time</title>
    <rect fill="black" y="50" height="200" width="400"/>
    <path d="M-4 50" stroke="white" stroke-width="1"/>
  </g>
 </g>
</svg>
]]></Attribute_String><Attribute_String roleName="setup">${MODELNAME_DEFAULT},${SVGURI_DEFAULT}</Attribute_String></SvgClient>

</XholonWorkbook>