loading...
running...
How to use
Pick a preset, then change anything. The coloured lines in the graphs update live; pause with Space and reset with R. Try: (1) Plateau winter: no PV for weeks, watch the diesel trace and the fuel counter; (2) raise insulation from 8 to 30 cm and read the fourth graph: the marginal saving per centimetre falls fast; (3) set Snowfall to 4 m/yr, clear auto-jack, and let a year pass at 10 days per second until the drifts bury the floor; then press Jack up; (4) lower the design snow load to 1 kPa or choose IPE 100 joists and the roof check turns red, then press Clear roof; (5) in Deep-core campaign watch the average rate fall with depth because the hoist takes longer; (6) switch on the thermal camera in the 3D view to see which surfaces leak the most (windows, then the unsnowed roof). Checks at the bottom compare the code against independent finite-difference and closed-form results.
The 3D page puts the same model on a Blender ice sheet: the sun, sky and aurora follow the latitude and date, snow piles up as a deformable drift around the camp and as a layer on the roofs, the platform rises on its screw jacks, the turbine, anemometer, winch and PV tilt move with the model, and the thermal camera colours each surface by its computed outside temperature. Drag to orbit, wheel to zoom, use the camera buttons or the guided tour. The station, power plant, drill and vehicles are Blender 5.2 models with baked PBR textures.
The science
What this lab models. A research station on an ice sheet has four coupled problems: the snow that falls and drifts on it, the heat it loses, the electricity it can make, and the ice core it exists to drill. All four are simulated together on a clock; the weather is a deterministic synthetic series (seasonal and diurnal cycles plus slow synoptic waves and storms), so a run is reproducible.
Sun and solar power. The solar elevation follows from latitude φ, declination δ and hour angle H:
sin h = sinφ sinδ + cosφ cosδ cos H, δ = 23.44° sin(2π(d - 81)/365)
Direct beam uses the Meinel clear-sky formula at 2.5 km altitude, DNI = 1353 [(1-0.14a) 0.7^(AM^0.678) + 0.14a] with air mass AM = 1/sin h, plus isotropic sky diffuse and a snow-reflected component (albedo 0.8). Clouds scale the result by 1 - 0.75 c^3.4. PV power is η A G [1 - 0.004 (T_cell - 25)], so cold panels are slightly better, but there is no sun at all in polar night.
Wind. P = ½ ρ A C_p v^3 between cut-in (3 m/s) and the rated speed (11 m/s), flat above it, zero beyond the 25 m/s cut-out. Cold air is denser (ρ = p/RT with the pressure of a 2.5 km plateau), the rotor area is sized so rated power is reached at 11 m/s.
Heat loss. Each element has a series-resistance U-value, with the outside film depending on wind:
U = (1 + b) / (R_si + t/k + R_snow + 1/(5.7 + 3.8 v))
Snow on the roof adds R_snow = d/k_snow, with the Sturm conductivity k = 0.138 - 1.01e-3 ρ + 3.233e-6 ρ^2: a snow blanket both loads the roof and insulates it. Ventilation loss is ρ c_p Q (1 - η_HRV) (T_in - T_out) with the larger of infiltration (ACH50/20, raised by wind) and the fresh-air need (8 L/s per person). The building is a single thermal mass (6 MJ/K) with a thermostat; internal gains are the crew (100 W each) plus every electric load that ends up as heat. Heat pump COP is a 42% fraction of the Carnot value between inside and outside, limited to 1 to 4.
Electric bus. Critical loads (base + heating) are served first from PV, wind, generator and battery; the drill gets what is left. The generator starts below 25% and stops above 80% charge and burns fuel on the Willans line F = 0.246 P_out + 0.08145 P_rated (L/h, kW). The battery is 95% efficient each way, its usable capacity falls in the cold, and it never discharges below empty: unmet demand first throttles the heater, then the base loads.
Snow and the roof. Snow accumulates as water equivalent (1 mm = 1 kg/m² = 9.81 Pa); the roof load is s = μ C_e C_t s_g (ISO 4355 form, C_t = 1 here). Steep roofs shed snow (μ = 0.8 up to 30°, to 0 at 60°). Drifts around a body are deeper than the open snow depth, and a raised platform lets wind scour underneath: the model uses d_drift = d (1 + 1.8 (1 - e^(-v/6)) (1 - 0.8 min(1, h/4))). When the drift reaches the floor the station must be jacked up (screw jacks, as on real stations). Roof joists are checked as simply supported beams: M = wL^2/8, σ = M/W, δ = 5wL^4/(384 E I), limit L/250.
Ice core. With uniform strain the Nye model gives the age at depth z in ice of thickness H and accumulation a (m of ice per year):
age(z) = (H/a) ln(H / (H - z))
Annual layers thin as a (1 - z/H). The temperature profile is the Robin solution of k T'' = -ρ c (a/H) z' T' with geothermal flux G at the bed:
T(z) = T_s + (G/k) (√π/2) l [erf(H/l) - erf((H - z)/l)], l = √(2kH/(ρ c a))
capped by the pressure-melting line (-0.67 K per km). Drilling is a cycle: cut one 3 m run, then hoist and lower the sonde, so the average rate is 3 / (t_cut + 2z/v + t_handling) and slows with depth; the winch needs (m g v)/η with a buoyancy-reduced cable weight.
Honest simplifications. Weather is synthetic, not a climate record. The building is one node; wind direction is fixed; snow density does not evolve; Robin assumes steady-state, vertically uniform strain and no sliding; no firn layer, no borehole closure, no real drilling fluid physics, and PV snow cover is ignored. The numbers are plausible engineering estimates for teaching, not a design basis. Not a model of any real station.