Thermal Model Documentation
Introduction
MARSTHERM is a Mars surface thermal model developed and modified for use in the determination of the thermal inertia of the Martian surface from MGS TES temperature data. For a given set of thermophysical parameters (thermal inertia, albedo, dust opacity, etc.) the model calculates the surface temperature of a location on Mars, as a function of the time. The calculated surface temperatures are then compared to observed temperatures to find the thermal inertia. This model was used as the basis for calculating Mars thermal inertia as presented by Mellon et al. (2000), Putzig et al. (2005), and Putzig and Mellon (2007).
The model employs a standard Mars subsurface thermal model, basic visible and thermal infrared radiative transfer through atmospheric CO2 and dust, and a coupling of these two model components as part of the surface boundary condition.
The code can be run in one of two modes - seasonal or diurnal. Both modes incorporate the same physics and use the same input files. The seasonal model computes the surface temperatures for a full Martian year (after a period of a few years for convergence) and the diurnal model computes temperatures for a single day (after a period of several days for convergence).
Details of the theory incorporated and a user's guide are given below.
Theory
The theory incorporated in the MARSTHERM model is primarily described in Haberle and Jakosky [1991], where the precurser model was used to evaluate the effects of atmospheric thermal radiation on the surface energy budget. Components of the atmospheric model originated from the Ames Mars GCM and are described in further detail by Pollack et al [1990] and references therein dating as far back as 1979. The surface temperature is modeled by a standard thermal model with modifications to explicitly calculate the contribution of atmospheric thermal radiation. While the model is similar to that of Haberle and Jakosky, it contains numerous improvements to enhance speed and to correct physics and numerical logic - these corrections and changes include, but are not restricted to: proper inclusion of subsurface conduction in the surface energy budget, iteration of subsurface temperatures when CO2 frost is present (with modifications for transitional behavior up to a user-specifiable threshold of frost on the ground), implementation of multiple latitude looping, and extension to a seasonal model. Also, a local surface slope angle and azimuth can be included which corrects for the sky angle and slant path through the atmosphere, but does not account for any radiative tranfer between differently sloped surfaces. While these modifications represent a substantial improvement, the methodology described in Haberle and Jakosky [1991] is the same. In addition, it is believed that these changes and corrections do not impact the validity of the results of Haberle and Jakosky [1991], but that are critical for proper TES thermal inertia determination.
Standard Thermal Model
Surface and subsurface temperatures are calculated by a standard thermal model, where the thermal diffusion equation is solved for subsurface temperatures with the appropriate boundary conditions. The lower boundary is below the depth of penetration of the seasonal thermal wave, and is assumed non-conducting, and the upper boundary is the surface. The surface temperature is calculated by a balance of absorbed solar radiation, subsurface conduction, seasonal CO2 condensation, thermal emission, and absorption of thermal radiation from the atmosphere (described below).
Incorporating these components, the surface and subsurface temperatures are iterated throughout each day and throughout the Martian year. Integral to solving the surface temperature is the atmospheric model described below.
Atmospheric Model
The atmospheric model incorporates solar and thermal infrared radiative transfer of airborne dust and CO2. These effects include attenuation of solar radiation incident on the surface, emission of thermal radiation (downward onto the surface and upward out the top of the atmosphere) and attenuation of emitted radiation from the surface to space as it passes through the atmosphere to space. The model also includes the potential effects of the condensation of CO2 within the atmosphere, convective instability, and sensible heat exchange with the surface.From this model, atmospheric temperatures are interated throughout each day and throughout the Martian year. These temperatures are used to calculate the thermal energy budget for the next time step as well as the thermal radiation incident on the surface. In addition to the absorption of surface emission by atmospheric components, emission by these same atmospheric components, which adds to the total thermal flux exiting the top ot the atmosphere, is calculated.
Both seasonal and diurnal models of MARSTHERM operate by a fully explicit forward-time centered-space finite-difference solution to the thermal diffusion equation for the surface and subsurface temperatures. A similar forward-time solution to the atmospheric temperatures is employed using numerous lookup tables to speed the radiative calculations.
General program flow proceeds as follows:
- Iterate Ground Temperatures
- Compute IR rediative balance with the 15µm band (dust and CO2)
- Compute IR radiative balance outside the 15µm band (dust only)
- Iterate the atmospheric temperatures
- Test and adjust for convective instability
- Test and adjust for conduction
- Compute sensible heat transfer
The sequence is repeated in a time-marching manner. To minimize CPU requirements while providing the necessary numerical accuracy, the surface and atmospheric temperatures are iterated at two different time stepping rates (By default the surface and subsurface temperatures are iterated 3 times more often than the atmosphere - though this is adjustable with the field 'Ground Time Steps Per Atmos Step', care should be exercised to maintain numerical accuracy and stability).
MARSTHERM can calculate thermal behavior for up to 37 latitudes simultaneously (set using the field 'Latitudes Per Run'. In doing so it iterates step 1 above through all the latitudes, then separately steps 2 through 7 for each latitude. At the end of each time step the code decides if output is desired and writes it to the appropriate output file accordingly.
Output is written only during the final year of the seasonal model or the final day of the diurnal model after an appropriate period of convergence which has been predetermined from numerical tests and tuning.
Inputs
The MARSTHERM executable code requires a text formatted file containing model parameters. This file is created from the values entered into the thermal-model interface page and once a run is complete it can be downloaded as the 'config.txt' file along with the other model output.
A description of the required elements is as follows:
| Field | Name in config.txt | Name in Output Files | Units | Default | Allowed Range | Info |
|---|---|---|---|---|---|---|
| Thermal Inertia | TI | I | tiu = J/(m2Ks0.5) | 256.43069 | 5 - 5000 | The assumed Thermal Inertia of the surface. |
| Albedo Of Bare Ground | ALGND | Agnd | none | 0.25 | 0 - 1 | The albedo of bare ground. |
| Albedo of CO2 Surface Ice | ALICE | Aice | none | 0.65 | 0 - 1 | The albedo of CO2 Surface Ice. |
| Effective Semi-Infinite CO2 | CO2ICET | IceT | kg/m2 | 6.0 | 0 - 1000 | The amount of CO2 ice above which ice properties, rather than rock, will dominate surface thermal behaviour. |
| Bare Ground Emissivity @ 15µm | EG15GND | E15 | none | 1.0 | 0 - 1 | The Emissivity of Bare Ground in the 15µm band, an absorption line of CO2. |
| Bare Ground Emissivity NOT @ 15µm | EGOGND | Eout | none | 1.0 | 0 - 1 | The Emissivity of Bare Ground outside of the 15µm band. |
| Surface Slope Angle | SLANG | SLANG | Degrees (°) from horizontal | 0.0 | 0 - 90° | The surface slope angle. |
| Surface Slope Azimuth | SLAZI | SLAZI | Degrees (°) East of North | 0.0 | 0 - 360° | The surface slope azimuth. |
| Surface Pressure | PSL | P | mbar, hPa | 6.0 | 0 - 1000 | The surface atmospheric pressure. |
| Atmospheric Dust Opacity | TAUDUST | TauD | m2/kg | 0.1 | 0 - 1000 | The Atmospheric Dust Opacity. |
| Lowest (Southern Most) Latitude | PHI0 | NA | degrees (°) | -90.0 | -90 - 90° | Lowest (Souther Most Latitude). |
| Latitude Increment | DELPHI | NA | degrees (°) | 5.0 | 0 - 180° | Increment between latitude outputs |
| Latitudes Per Run | JO | NA | none | 1 | 1 - 37 | Number of latitudes for which the output are produced. 'JO' outputs will be produced at 'DELPHI' spacing up from 'PHI0'. |
| Ground Time Steps Per Day | NT | NT | none | 144 | 1 - 10000 | The number of timesteps for which the model is calculated for ground temperatures during each day. 144 is roughly once every 10 minutes. |
| Ground Time Steps Per Atmos Step | NFREQ | NFQ | none | 3 | 1 - 99 | To improve performance, the atmospheric temperatures are calculated 'NFREQ' times less frequently than the ground ones. |
| Mode | M_FLAG | NA | none | 1 | 0 or 1 | 0 = 'Seasonal', 1 = 'Diurnal'. See below for details. |
| Solar Longitude | LSUBS | Ls | degrees (°) | 180.0° | 0 - 360° | Position in orbit of Mars, signifying season or date. Degrees (forward) from Northern Spring Equinox. i.e. 90° signifies the summer solstace. Diurnal mode only, ignored otherwise. |
| Interval | IDAY | NA | days | 8 | 1 - 668 | In seasonal mode the output will be generated as a table of diurnal temperatures produced for every IDAY-th day of the year. Seasonal mode only, ignored otherwise. |
| Print Last Day | PLAST | NA | none | No | 0 or 1 | One Mars year is around 669 Mars Days. If there is not an integer number of 'IDAY' dividers in one Mars Year, then setting 'PLAST' to yes (1=yes, 0=no) will force the output of the final day's output. Seasonal mode only, ignored otherwise. |
| Output Interval | PSTEP | NA | none | 6 | 1 - 99 | Number of ground steps between outputs. i.e. if NT=144, and PSTEP=6, there will be 144/6=24 outputs per day. |
| Print Type | P_FLAG | NA | none | 1 | 0 or 1 | 0 = 'Research', 1 = 'Production'. See below for details. |
Hard Wired Parameters
As well as those parameters that can be changed by users, the following parameters used in the generation of the output files are hard wired into the code:
| Field | Value | Units | Info |
|---|---|---|---|
| Soil Heat Capacity | 627.900 | J kg-1 K-1 | |
| IR/visible dust opacity at 15µm | 0.500 | m2/kg | |
| IR/visible dust opacity not at 15µm | 0.500 | m2/kg | |
| CO2 ice emissivity at 15µm | 0.800 | none | |
| CO2 ice emissivity not at 15µm | 0.800 | none | |
| Soil Density | 1.5×103 | kg m-3 | |
| Initial Temperature | 200 | K | |
| Wind speed | 0.000 | m s-1 | |
| Seasonal model runs | 3 | years | For convergence, one year output. |
| Diurnal model runs | 10 | days | For convergence, one day output. |
Many radiative properties of the atmosphere cannot be varied with the existing code as much of these calculations were performed off line and have been incorporated simply as lookup tables.
Output
Model output is written in two text formats. One format is termed 'production' output and contains output required for generating a lookup table for 'on the fly' thermal inertia determination from TES data. The second format is termed 'research' and contains a zeroth order set of output that may provide more information under some circumstances, including details about atmospheric temperatures.
One or the other output type is written depending on the value of the field 'Print Type'. Production output is written as a text formatted file. Because MARSTHERM can run up to 37 latitudes simultaneously, one file is written for each latitude in a single run with the lowest unit number corresponding to the southern most latitude. For example, if the code is to run 3 latitudes (5, 30, and 55 degrees), output will be written to files output.production.01, output.production.02 and output.production.03 corresponding to each latitude respectively. Research output works identically except that 'production' is replaced with 'research' in the file name.
Production Output
Production output format contains a limited set of results needed for thermal inertia lookup tables. In diurnal mode the output is written as a diurnal cycle at uniform increments throughout the day. In seasonal mode the output consists of a series of diurnal outputs at uniform increments through the Martian year.
Each day's output is preceded by a set of basic parameters. An example of this format is:
-----
DAY Ls | TauD Agnd I P Lat IceT NT NFQ SLANG SLAZI
0 0.0 | 0.00 0.32 215.0 7.00 23.00 6.00 144 3 0.000 0.00
TIME T_S(K) T_B(K)
3699. 179.39 177.97
7398. 177.56 176.42
11097. 175.93 175.02
14796. 174.45 173.75
18495. 173.11 172.59
22194. 171.87 171.53
25893. 186.22 183.39
29592. 207.68 201.84
33291. 227.64 219.65
36990. 243.70 234.65
40689. 255.07 246.15
44388. 262.11 251.70
48086. 264.17 253.68
51785. 261.89 251.58
55484. 255.40 245.47
59183. 244.70 235.52
62882. 229.70 221.78
66581. 210.12 204.29
70280. 199.85 195.34
73979. 194.17 190.47
77678. 189.98 186.93
81377. 186.65 184.11
85076. 183.88 181.77
88775. 181.50 179.76
-----
DAY Ls | TauD Agnd I P Lat IceT NT NFQ SLANG SLAZI
17 8.6 | 0.00 0.32 215.0 7.00 23.00 6.00 144 3 0.00 0.00
TIME T_S(K) T_B(K)
3699. 179.95 178.52
7398. 178.10 176.94
11097. 176.45 175.53
etc
etc
For Diurnal runs only one day of data is produced. For Seasonal runs output for different days is seperated by a line that reads '-----'.
Note that time is listed in seconds since midnight and the Mars day is 88775 seconds in length. The surface temperature T_S is the surface kinetic temperature as computed by the model. The brightness temperature T_B is the temperature that corresponds to the exiting infrared flux from the top of the the Martian atmosphere (as would be measured from orbit). The brightness temperature accounts for absorption of the thermal energy emitted from the surface, as well as the addition of atmospheric thermal emission. In calculating the brightness temperature from the emitted flux, an emissivity of 1 was assumed. All temperatures are in Kelvin.
Day is the number of the day since the Mars spring equinox; hard wired to 10 for diurnal runs, variable in seasonal runs. Other values are as defined above.
Research Output
Research output is similar but a bit more verbose, in that it includes additional basic parameters and additional diurnal results pertaining to the atmospheric temperature structure. An example of this format is:
-----
In this case more input parameter information is given, including the subsolar declination (SDEC) and the orbital distance in au (RAU). Pressures are given in mb, IRTOP is the total exiting infrared flux from the top of the atmosphere in W/m^2, IRBOT is the total downwelling flux at the bottom of the atmosphere in W/m^2, IR15NL is the 15 um band downwelling flux at the bottom of the atmosphere in W/m^2, and time is now the time of day in hours of a 24 hour Martian day. Everything else is as for the production output, incluiding the '
DAY Ls | SDEC RAU | TauD Agnd I P Lat E15 Eout Aice IceT NT NFQ SLANG SLAZI
0 0.0 | 0.00 1.5580 | 0.10 0.25 256.4 5.0000 30.00 1.00 1.00 0.65 6.00 144 3 0.00 0.00
ATM LAYER # 16 15 14 13 12 11 10 9 8 7 6 5 4 3 2 1
PRESSURE 4.8709 4.5163 3.9571 3.2416 2.4975 1.8316 1.2939 0.8889 0.5983 0.3968 0.2605 0.1699 0.1106 0.0721 0.0474 0.0316
TIME T_SURF IRTOP IRBOT IR15NL T_ATM ...
1.00 183.34 63.36 17.23 16.05 204.0 216.5 213.9 206.7 199.8 193.0 187.2 180.2 175.2 168.2 164.6 157.7 156.1 148.6 149.6 151.7
2.00 181.57 61.27 16.58 15.41 201.9 215.0 213.1 206.3 199.5 192.7 187.0 179.9 175.0 167.7 164.4 157.1 155.8 148.1 149.4 151.6
3.00 179.98 59.44 16.00 14.84 200.1 213.6 212.2 205.8 199.1 192.4 186.7 179.5 174.7 167.3 164.1 156.6 155.5 147.7 149.2 151.2
4.00 178.54 57.80 15.47 14.33 198.3 212.3 211.3 205.4 198.8 192.0 186.4 179.1 174.5 166.9 163.9 156.2 155.1 147.3 148.9 150.7
etc
etc
-----' line to seperate different day's outputs.
Plotted Output
To aid fast interpretation of output, .png plots of surface and brightness tempertatures (for production output) or surface and atmospheric temperatures plus IR flux (for research output) are also produced for download, or viewing in a web browser.Uncertainties
An important but difficult factor to assess in such a model is the uncertainty in the output. Two sources of uncertainty exist: 1) uncertainty in the type of physics included in the model and how the physics is included, and 2) the inherent inaccuracy in the numerical techniques employed.
To assess the uncertainty in the physical model, consider that the model assumes a semi-infinite half-space with uniform properties, and a 1-D approach to understanding the physics. Anything that causes the surface to depart from these assumptions introduces uncertainties or errors into the meaning of the model temperatures. As these errors are difficult to quantify, and are obviously the subject of ongoing analysis, we only list some of the issues here to make the novice aware of them:
- The atmospheric temperatures are derived assuming radiative exchange of energy. Any effects of atmospheric dynamics (beyond vertical convection) on the temperatures are ignored. These can be substantial, especially higher up in the atmosphere.
- Seasonal variations in atmospheric pressure resulting from either condensation/sublimation in the polar regions or atmospheric dynamics are ignored. These are not very significant effects.
- Temporal variations (such as seasonal, or day-to-day) in the atmospheric dust opacity are ignored.
- Emissivity and albedo are cosine-tapered from ground to ice conditions according to the user-specified threshold quantity of CO2 ice (in kg/m^2). This represents the minimum amount of CO2 ice which is effectively semi-infinite with respect to albedo and emissivity. Other radiative effects of atmospheric condensates are ignored.
- Possible layering (vertical structure) within the atmosphere (e.g., dust) and within the subsurface are ignored. Obviously, this can be an important effect.
- Any effects of local topography, such as the small-scale spatial structure resulting from the presence of rocks, is ignored.
- Any local or regional topography is ignored.
- Any radiative transfer between differently sloped surfaces is ignored.
- Any variation of brightness temperature with wavelength that results from spatial structure of the surface is ignored.
Inferring the temperature of a given component of the surface, or using the wavelength dependence of observed brightness temperature to infer surface rock abundance, can be done. However, it requires an understanding of the physics and energy transfer of the surface. Also, it is clear from the Viking results that the derived thermal inertia has a strong connection to the geology of the surface and the properties of the near-surface layer, although those connections are not well understood or unique.
The second of the two types of uncertainty is easier to evaluate based on sensitivity to numerical parameters and comparison with simplified analytical solutions. The surface and subsurface temperature calculations are estimated to give a maximum numerical error in the surface temperature of about 0.25 K or less. Numerical error induced by the atmospheric component of the model is larger, at about 1K or less (being largest during the mid-day). This error can be reduced by decreasing the length of the atmospheric time steps, but at significant cost in run time.