Thermal Inertia for Mars Analogs (TIMA) Model

TIMA is a finite element numerical model of soil heat transfer that that considers time-variability in thermal conductivity as a function of temperature, moisture, and composition. The detailed parametrization of thermal conductivity is meant to enable direct translation between soil thermophysical properties on Earth and Mars. This repository is maintained by Ari Koeppel. Please cite Koeppel et al., (2024) when using the model.

Note: the 2D implementation of the model using raster data is requiring some reconstruction and will be available here soon!

References:

Koeppel, A.H., Edwards, C.S., Edgar, L.A., Nowicki, S., Bennett, K.A., Gullikson, A., Piqueux, S., Eifert, H., Chapline, D. and Rogers, A.D., 2024. A novel surface energy balance method for thermal inertia studies of terrestrial analogs. Earth and Space Science, 11(9), p.e2023EA003259.

Overview:

The original implementation of the model was designed to derive a characteristic dry soil thermal conductivity at 300K. The procedure achieves this by inputing micrometeorological data into a surface energy balance forward model and adjusting thermal conductivity (along with 5 other modifying parameters) to fit surface skin temperature data (typically obtained through radiometer observation). The preferred fitting protocol involves Markov Chain Monte Carlo simulations to identify the most probable set of constants and associated uncertainties. Model updates have included multilayer parametrizations, including consideration of subsurface ice. Thus, the model can be used to derive soil physical properties, layering and the depth to subsurface material transitions based entirely on surface temperature and weather observations.

The required fitting parameters are: 1. kdrystdupper: [W/mK] bulk dry thermal conductivty of upper layer at Tstd (scalar) 2. m: [unitless] pore-network-connectivity-parameter, typically 0-1.3 (scalar) 3. CH: [Unitless] resistance to sensible heat flux coefficient, similar to the aerodynamic scaling factor rhoair*Cpair/(log(z1/z0)^2/Kv^2) 1/CH should be 0.0028-0.0075 (or CH~100-400) for smooth to roughly open soils on Davenport Scale, CH is larger (more resistence) with rougher topography or larger vegetation (scalar) 4. CH: [Unitless] resistance to sensible heat flux coefficient, similar to the aerodynamic scaling factor rhoair*Cpair/(log(z1/z0)^2/Kv^2) 1/CH should be 0.0028-0.0075 (or CH~100-400) for smooth to roughly open soils on Davenport Scale, CH is larger (more resistence) with rougher topography or larger vegetation (scalar) 5. thetak: [0-1, fraction by volume] conductivity soil moisture inflection point - in theory this should be similar to saturation value or porosity (scalar) 6. thetaE: [0-1, fraction by volume] latent heat soil moisture inflection point (scalar)

and optional parameters:

1. 'depth_transition': [m] Depth of major transition between upper material and lower material. (scalar)
2. 'kdrystdlower': [W/mK] bulk dry thermal conductivty of lower layer at Tstd (scalar)

Forward Model:

[T_surf_C,T_sub_C,q_latent,k_eff_dt,q_conv,q_rad,q_G] = tima_heat_transfer(k_dry_std_upper,m,CH,CE,theta_k,theta_E,rho_dry_upper,dt,T_std,air_temp_C,r_short_upper,r_short_lower,r_long_upper,  windspeed_horiz,T_deep,initial_temps,layer_size,VWC_column,evap_depth,RH,emissivity,pressure_air_pa,varargin)

Input Parameters ('single quotes' means optional varargin)

### Timeseries data airTempC: [C] near surface air temperature, typically at 3 m AGL. (1D vector) 'albedo': [0-1] time variant albedo (e.g., due to wetting) for fitting surface. (1D vector) evapdepth: [m] Depth of the evaporation front (e.g., as function of VWC). (1D vector) 'fdiff': [0-1] Fraction of Global Horizontal Irradiance (GHI) or rshortupper that is diffuse. (1D vector) RH: [0-1] array of near surface relative humidity values, typically at 3m AGL. (1D vector) pressureairpa: [Pa] station pressure, typically at 3m AGL. (1D vector) rlongupper: [W/m^2] Integrated longwave radiation (4.5 to 42 μm) incident on flat surface.(1D vector) rshortlower: [W/m^2] Integrated upwelling shortwave radiation (305 to 2800 nm) from flat surface. (1D vector) rshortupper: [W/m^2] Integrated shortwave radiation (305 to 2800 nm) incident on flat surface. (1D vector) 'shadowtimeind': Index of shadowdata corresponding to timeseries. (1D vector) 'solarazimuthcwfromS': [degrees] Solar azimuth in degrees clockwise from South, typically -180:180 or 0:360 (1D vector) 'solarzenithapparent': [degrees] Solar zenith in degrees, corrected for atmospheric refraction. (1D vector) VWCcolumn: [0-1, decimal fraction by volume] array of volumetric water content for each model layer with each time step, typically interpolated. (2D vector) windspeed_horiz: [m/s] Near surface horizontal wind speed, typically at 3m AGL. (1D vector)

### Constants: 'aspectcwfromS': [degrees] slope aspect clockwise from South. (scalar) CE: [Unitless] resistance to latent heat flux coefficient, similar to the aerodynamic scaling factor rhoairCpair/(log(z1/z0)^2/Kv^2) (scalar) CH: [Unitless] resistance to sensible heat flux coefficient, similar to the aerodynamic scaling factor rhoairCpair/(log(z1/z0)^2/Kv^2) 1/CH should be 0.0028-0.0075 (or CH~100-400) for smooth to roughly open soils on Davenport Scale, CH is larger (more resistence) with rougher topography or larger vegetation (scalar) 'depthtransition': [m] Depth of major transition between upper material and lower material. (scalar) dt: [s] Time step (scalar) emissivity: [0-1] Weighted thermal emissivity over wavelength range of sensor. (scalar) 'efxn': spectral emissivity function that determines emissivity as function of wavelength in um. (function) initialtemps: [K] Initial condition, list of center temperatures for each layer at start of simulation (1D vector) 'kdrystdlower': [W/mK] bulk dry thermal conductivty of lower layer at Tstd (scalar) 'kdrystdmantle': [W/mK] bulk dry thermal conductivty of topmost mantling layer at Tstd (scalar) kdrystdupper: [W/mK] bulk dry thermal conductivty of upper layer at Tstd (scalar) layersize: [m] List of vertical thickness for each layer from top to bottom. (1D vector) m: [unitless] pore-network-connectivity-parameter, typically 0-1.3 (scalar) 'mantlethickness': [m] thickness of topmost mantling layer. (scalar) 'material': ['basalt' 'amorphous' 'granite' 'clay' 'salt' 'ice'] primary mineralogy at the surface (char) 'materiallower': ['basalt' 'amorphous' 'granite' 'clay' 'salt' 'ice'] primary mineralogy at depth (char) NDAYS: [unitless] # of days to run equilib model for: More = closer to equilib, fewer = faster (scalar) 'rhodrylower': [kg/m^3] Value for density of lower layer soil. (scalar) rhodryupper: [kg/m^3] Value for density of soil. (scalar) 'shadowdata': [0-1] Fraction of ROI/pixel in shadow corresponding to time. (1D Vector) 'slopeangle': [degrees] angle of slope. (scalar) thetak: [0-1, fraction by volume] conductivity soil moisture inflection point - in theory this should be similar to saturation value or porosity (scalar) thetaE: [0-1, fraction by volume] latent heat soil moisture inflection point (scalar) Tdeep: [K] Lower boundary condition, fixed temperature (scalar) Tstd: [K] Standard temperature; typically 300 (scalar) 'Tadj1': [index, temperature K] pair used to force column temperature change at a given time point due to wetting (1D vector) 'T_adj2': [index, temperature K] pair used to force a second column temperature change at a given time point due to wetting (1D vector)

Output Parameters:

T_Surf_C = [C] Surface temperature time series (1D vector)
T_sub_C = [C] Suburface temperature time series for all layers (2D vector)
q_latent = [W/m^2] Latent heat flux time series for all layers (2D vector)
k_eff_dt = [W/mK] Effective thermal conductivity time series for all layers (2D vector)
q_conv = [W/m^2] Surface sensible heat flux time series (1D vector)
q_rad = [W/m^2] Surface radiative heat flux time series (1D vector)
q_G = [W/m^2] Heat flux between top and second subsurface layer time series (1D vector)

Fitting Process Input/Output Setup:

### Time data - struct of timeseries data variables TData.airTempC: [C] near surface air temperature, typically at 3 m AGL (1D vector) TData.DF: [decimal fraction] Fraction of Global Horizontal Irradiance (GHI) or rshortupper that is diffuse (1D vector) TData.evapdepth: [m] Depth of the evaporation front. (1D vector) TData.tempstofitinterp: [C] Surface temperature values to be used for fitting with any data gaps filled in. (1D vector) TData.timedalbedo: [decimal fraction] time variant albedo (e.g., due to wetting) for fitting surface. (1D vector) TData.TIMESTAMP: datetime array associated with each table row spaced by Mdata.dt s TData.humidity: [decimal fraction] array of near surface relative humidity values, typically at 3m AGL TData.rlongupper: [W/m^2] Integrated longwave radiation (4.5 to 42 μm) incident on flat surface (1D vector) TData.pressureairPa: [Pa] station pressure, typically at 3m AGL (1D vector) TData.rlongupper: [W/m^2] Integrated longwave radiation (4500 to 42000 nm) incident on flat surface (1D vector) TData.rshortlower: [W/m^2] Integrated upwelling shortwave radiation (305 to 2800 nm) from flat surface (1D vector) TData.rshortupper: [W/m^2] Integrated shortwave radiation (305 to 2800 nm) incident on flat surface (1D vector) TData.solarazimuthcwfromS: [degrees] Solar azimuth in degrees clockwise from South, typically -180:180 or 0:360 (1D vector) TData.solarzenithapparent: [degrees] Solar zenith in degrees, corrected for atmospheric refraction. (1D vector) TData.VWCcolumn: [decimal fraction by volume] array of volumetric water content for each model layer with each time step, typically interpolated. (2D vector) TData.windspeedhorizms: [m/s] Near surface horizontal wind speed, typically at 3m AGL. (1D vector)

### Model Data - Struct of static and model format variables MData.burninfit: initial time length (s) to ignore in fitting (vector, default=0) MData.burninmcmc: fraction of the chain that should be removed. (scalar, default=0) MData.density: [kg/m^3] Value for density of soil beneath tower. (scalar) MData.dt: [s] Time step (scalar) MData.emissivity: [0-1] Weighted thermal emissivity over wavelength range of sensor. (scalar) MData.erf: Uncertainty as function of observed temperature (functionhandle) MData.fitind: Indecies of tempstofitinterp in which to apply fitting to (scalar) MData.kdrystdmantle: [W/mK] bulk dry thermal conductivty of topmost mantling layer at Tstd (scalar) MData.layersize: [m] List of vertical thickness for each layer from top to bottom. (1D vector) MData.mantlethickness: [m] thickness of topmost mantling layer. (scalar) MData.material: ['basalt' 'amorphous' 'granite' 'clay' 'salt' 'ice'] primary mineralogy at the surface (char) MData.materiallower: ['basalt' 'amorphous' 'granite' 'clay' 'salt' 'ice'] primary mineralogy at depth (char) MData.minit: Vector of initialized test variables with as many randomized samples as desired for fitting, 50 is good such that vector is nvarx50 (2D vector) MData.notes: Details to record in data structure (string) MData.nstep: Number of iterations for curve fitting, 250 is good (scalar) MData.nvars: Number of variables being fit for (scalar) MData.nwalkers: Number of walkers in MCMC ensemble (scalar) MData.ThinChain: MCMC data reduction by thinning all output chains by only storing every N'th step (scalar, default=10) MData.Tadj1: [index, temperature K] pair used to force column temperature change at a given time point due to wetting (1D vector, optional) MData.Tadj2: [index, temperature K] pair used to force a second column temperature change at a given time point due to wetting (1D vector, optional) MData.Tdeep: [K] Lower boundary condition, fixed temperature (scalar) MData.Tstart: [K] Initial condition, list of center temperatures for each layer at start of simulation (scalar) MData.Tstd: [K] Standard temperature; typically 300 (scalar) MData.varsinit: [k-upper [W/mK], Pore network con. par. (mk) [unitless],... Surf. ex. coef. (CH) [unitless], Surf. ex. coef. (CE) [unitless], Soil Moist. Infl. (thetak) [% by volume],... Soil Moist. Infl. (thetaE) [% by volume], (Transition Depth [m]), (k-lower [W/mK])] List of 6-8 inputs for variables to serve as either initial or fixed values. (1D vector)

Output directory:

  out_DIR: Full path to directory for outputs (string)

Input data:

Example input files are located in the Example Data folder.

Run Fitting Process

[models,names] = tima_TI_Earth_Tower(TData,MData,outDIR,'Mode','1layer','Parallel',true,'SaveModels',true);

Visualize Results

tima_plot_results(TData,MData,models,names,varargin)