Twin Peaks, Mars Pathfinder Image

Site developed by Eddy Barratt, Sam Van Kooten, Than Putzig, and Matthew Perry

Funding provided by NASA MFRP Grant NNX11AM33G
PGG Grant NNX14AN09G
PGG Grant NNX17AB25G

South West Research Institute
University of Colorado, Boulder

Thermal Model Documentation


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.


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.

Users Guide

Program Flow

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:

  1. Iterate Ground Temperatures
  2. Compute IR rediative balance with the 15µm band (dust and CO2)
  3. Compute IR radiative balance outside the 15µm band (dust only)
  4. Iterate the atmospheric temperatures
  5. Test and adjust for convective instability
  6. Test and adjust for conduction
  7. 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.


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:

FieldName in config.txtName in Output FilesUnitsDefaultAllowed RangeInfo
Thermal InertiaTIItiu = J/(m2Ks0.5)256.430695 - 5000The assumed Thermal Inertia of the surface.
Albedo Of Bare GroundALGNDAgndnone0.250 - 1The albedo of bare ground.
Albedo of CO2 Surface IceALICEAicenone0.650 - 1The albedo of CO2 Surface Ice.
Effective Semi-Infinite CO2CO2ICETIceTkg/m26.00 - 1000The amount of CO2 ice above which ice properties, rather than rock, will dominate surface thermal behaviour.
Bare Ground Emissivity @ 15µmEG15GNDE15none1.00 - 1The Emissivity of Bare Ground in the 15µm band, an absorption line of CO2.
Bare Ground Emissivity NOT @ 15µmEGOGNDEoutnone1.00 - 1The Emissivity of Bare Ground outside of the 15µm band.
Surface Slope AngleSLANGSLANGDegrees (°) from horizontal0.00 - 90°The surface slope angle.
Surface Slope AzimuthSLAZISLAZIDegrees (°) East of North0.00 - 360°The surface slope azimuth.
Surface PressurePSLPmbar, hPa6.00 - 1000The surface atmospheric pressure.
Atmospheric Dust OpacityTAUDUSTTauDm2/kg0.10 - 1000The Atmospheric Dust Opacity.
Lowest (Southern Most) LatitudePHI0NAdegrees (°)-90.0-90 - 90°Lowest (Souther Most Latitude).
Latitude IncrementDELPHINAdegrees (°)5.00 - 180°Increment between latitude outputs
Latitudes Per RunJONAnone11 - 37Number of latitudes for which the output are produced. 'JO' outputs will be produced at 'DELPHI' spacing up from 'PHI0'.
Ground Time Steps Per DayNTNTnone1441 - 10000The 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 StepNFREQNFQnone31 - 99To improve performance, the atmospheric temperatures are calculated 'NFREQ' times less frequently than the ground ones.
ModeM_FLAGNAnone10 or 10 = 'Seasonal', 1 = 'Diurnal'. See below for details.
Solar LongitudeLSUBSLsdegrees (°)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.
IntervalIDAYNAdays81 - 668In 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 DayPLASTNAnoneNo0 or 1One 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 IntervalPSTEPNAnone61 - 99Number of ground steps between outputs. i.e. if NT=144, and PSTEP=6, there will be 144/6=24 outputs per day.
Print TypeP_FLAGNAnone10 or 10 = '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:

Soil Heat Capacity627.900J kg-1 K-1
IR/visible dust opacity at 15µm0.500m2/kg
IR/visible dust opacity not at 15µm0.500m2/kg
CO2 ice emissivity at 15µm0.800none
CO2 ice emissivity not at 15µm0.800none
Soil Density1.5×103kg m-3
Initial Temperature200K
Wind speed0.000m s-1
Seasonal model runs3yearsFor convergence, one year output.
Diurnal model runs10daysFor 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.


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:

0 0.0 | 0.00 0.32 215.0 7.00 23.00 6.00 144 3 0.000 0.00

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
17 8.6 | 0.00 0.32 215.0 7.00 23.00 6.00 144 3 0.00 0.00

3699. 179.95 178.52
7398. 178.10 176.94
11097. 176.45 175.53

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:

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

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

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 '-----' 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.


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:

  1. 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.
  2. 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.
  3. Temporal variations (such as seasonal, or day-to-day) in the atmospheric dust opacity are ignored.
  4. 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.
  5. Possible layering (vertical structure) within the atmosphere (e.g., dust) and within the subsurface are ignored. Obviously, this can be an important effect.
  6. Any effects of local topography, such as the small-scale spatial structure resulting from the presence of rocks, is ignored.
  7. Any local or regional topography is ignored.
  8. Any radiative transfer between differently sloped surfaces is ignored.
  9. Any variation of brightness temperature with wavelength that results from spatial structure of the surface is ignored.
So with all of the obvious oversimplifications, what is the temperature or the thermal inertia that is derived from it good for? Clearly, there is no single meaning to the question of what is the temperature of the surface. Most simply put, the temperatures (and the thermal inertia) reproduce the observed brightness temperatures at whatever wavelength, etc. They also reproduce with reasonable accuracy the subsurface temperatures, as those represent the average of the surface temperatures, convolved with the conduction to depth.

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.

model_documentation.php last modified by Than on 2016-04-04 at 13:38:15.