Parametric Monte Carlo Simulation with Temporal Weather Correlation

A hybrid approach combining ensemble weather perturbation with time-varying forecast grids for maritime voyage uncertainty quantification

Abstract

Weather routing systems for commercial shipping typically produce a single deterministic optimal route based on a single weather forecast. In practice, forecast uncertainty grows with lead time and vessels encounter conditions that deviate from the forecast, making deterministic ETA and fuel estimates unreliable for operational planning. This article presents a parametric Monte Carlo method that hybridizes two approaches from the academic literature: (1) temporally correlated stochastic perturbations applied to base weather fields, and (2) ensemble-like multi-timestep weather sampling from pre-ingested forecast grids stored in a database. The method produces P10/P50/P90 confidence intervals for ETA, fuel consumption, and voyage time. The implementation runs 100 simulations in under 500 milliseconds by pre-fetching weather data from compressed PostgreSQL grids and using Cholesky-decomposed correlation matrices for efficient correlated random variate generation.

1. Introduction

Maritime weather routing aims to find the path and speed profile that minimizes fuel consumption (or time, or a weighted combination) subject to safety constraints. The problem has been studied extensively since Hanssen and James (1960) first applied dynamic programming to ocean routing. Modern systems use numerical weather prediction (NWP) forecasts — typically from GFS, ECMWF, or CMEMS — to compute wind, wave, and current conditions along candidate routes.

A fundamental limitation of deterministic routing is that NWP forecasts are uncertain. Wind speed forecast errors grow roughly linearly with lead time: GFS 10-metre wind has an RMS error of ~1.5 m/s at 24 hours, increasing to ~3.5 m/s at 120 hours. Significant wave height errors follow a similar pattern. For a 5–7 day trans-oceanic voyage, the weather encountered in the final days may differ substantially from the forecast available at departure.

Two broad approaches exist in the literature for handling forecast uncertainty in ship routing:

1. Ensemble-based methods — Use multiple NWP ensemble members (e.g., GEFS 30-member ensemble, ECMWF 51-member ensemble) as distinct weather scenarios, optimizing a route for each or computing statistics across the ensemble (Dickson et al., 2019; Mannarini et al., 2019).
2. Parametric perturbation methods — Apply stochastic perturbations to a single deterministic forecast, with the perturbation magnitude calibrated to match observed forecast error statistics (Aijjou et al., 2021; Vettor & Guedes Soares, 2016).

Each approach has trade-offs. Ensemble methods capture physically consistent weather patterns but require downloading and processing 30–50 complete forecast fields per optimization cycle, which is prohibitively expensive for real-time routing. Parametric methods are computationally inexpensive but risk generating physically inconsistent weather scenarios if perturbations are applied independently across time and space.

WindMar adopts a hybrid approach that combines elements of both:

• Time-varying base weather from multi-timestep forecast grids (ensemble-like temporal variation)
• Temporally correlated parametric perturbations (Cholesky-decomposed exponential correlation)
• Cross-parameter coupling (wind–wave physical correlation)

This produces physically plausible weather scenarios at a fraction of the computational cost of full ensemble routing, while maintaining temporal coherence that independent perturbations lack.

2. Literature Review

2.1 Ensemble Approaches

Dickson et al. (2019) provide the most comprehensive treatment of uncertainty in marine weather routing. They use the GEFS 30-member ensemble to generate 30 distinct weather scenarios, compute optimal routes for each, and analyse the distribution of fuel consumption and arrival times. Their key finding is that route choice uncertainty (which path is optimal) is often more significant than weather uncertainty along a fixed route. However, they note that ensemble-based routing is computationally expensive, requiring 30× the weather data and computation of a deterministic system.

Mannarini et al. (2019) apply ensemble methods to the VISIR routing system, using ECMWF ensemble forecasts for Mediterranean routing. They demonstrate that ensemble spread provides a useful proxy for routing confidence and that routes optimised for the ensemble mean tend to be more robust than routes optimised for the control forecast.

2.2 Parametric Perturbation Methods

Aijjou et al. (2021) propose a comprehensive approach to weather uncertainty in ship route optimization using Monte Carlo simulation with parametric perturbations. They apply multiplicative factors to wind speed and wave height drawn from calibrated distributions, and show that the resulting fuel consumption distributions closely match those obtained from historical reanalysis data. A key contribution is their use of correlated perturbations across weather parameters (e.g., wind and waves), reflecting the physical relationship between wind forcing and sea state development.

Vettor and Guedes Soares (2016) study the impact of weather uncertainty on voluntary speed reduction decisions, finding that wave height uncertainty is the dominant factor in fuel consumption variance for laden tankers, while wind uncertainty dominates for container vessels with high windage.

2.3 Temporal Correlation

A critical insight from the meteorological literature is that forecast errors are temporally correlated. If the forecast underestimates wind speed at hour 12 of a voyage, it is very likely to also underestimate wind speed at hour 15, because the same synoptic-scale weather system persists. The decorrelation time scale is on the order of 1–3 days for mid-latitude synoptic systems (Holton & Hakim, 2013).

Papers on wind power forecasting uncertainty (e.g., Pinson & Madsen, 2012) have long used temporally correlated perturbation models for wind speed, typically with exponential or Gaussian correlation kernels. The maritime routing literature has been slower to adopt this approach, with most parametric methods applying independent perturbations per leg or per time step.

2.4 Gap Addressed

The approach implemented in WindMar bridges several gaps in the existing literature:

• Temporal correlation of perturbations (following wind power literature, applied to ship routing)
• Cross-parameter coupling (wind–wave correlation, following Aijjou et al.)
• Time-varying base weather from multi-timestep forecast grids (following ensemble spirit, but from a single deterministic forecast at multiple lead times)
• Computational efficiency enabling real-time use (100 simulations < 500ms)

3. Route Discretization

The voyage route is discretized into n time slices, evenly spaced along the route duration. The number of slices adapts to the voyage length:

This produces approximately one slice per 1.2 hours, capped at 100 slices for computational efficiency. For a typical 5-day MR tanker voyage (~120 hours), this yields 100 slices. For a short coastal voyage (~24 hours), 20 slices.

Each slice i has:

• Timestamp t_i = departure + (i / (n-1)) × T
• Position (lat_i, lon_i) via linear interpolation along the cumulative distance table of the route legs

The interpolation follows rhumb-line segments between waypoints, consistent with the voyage calculator's leg geometry.

4. Weather Pre-fetch from Forecast Grids

A distinguishing feature of this implementation is the use of pre-ingested multi-timestep forecast grids stored in PostgreSQL. Rather than using a single weather snapshot for all time slices, the simulator loads wave forecast data at multiple forecast hours (0, 3, 6, …, 120h) and selects the appropriate forecast hour for each time slice.

4.1 Data Sources per Slice

4.2 Forecast Hour Selection

For each time slice with timestamp t_i, the simulator computes the offset from the forecast run time and selects the closest available forecast hour:

This means a voyage departing 6 hours after the latest forecast run will use forecast hour f006 for its first slices and f120 for slices ~114 hours later. The base weather thus naturally varies along the route, reflecting the forecast model's prediction of how conditions will evolve — a form of implicit ensemble sampling from a single deterministic forecast at different lead times.

4.3 Bilinear Interpolation

Weather values at each slice position are obtained via bilinear interpolation from the forecast grid. The GridWeatherProvider._interp() static method performs nearest-neighbour index lookup followed by bilinear weighting:

This runs in microseconds per point, enabling 100 × 100 = 10,000 weather lookups (100 simulations × ~100 slices queried by the voyage calculator) without measurable overhead.

5. Temporal Correlation Model

The central innovation of this implementation is the use of temporally correlated perturbations via Cholesky decomposition of an exponential correlation matrix. This ensures that weather perturbations at nearby time slices are strongly correlated, producing physically plausible weather evolution within each simulation.

5.1 Exponential Correlation Kernel

The correlation between perturbations at time slices i and j is modeled by an exponential kernel:

The choice of τ = 0.3 corresponds to a decorrelation time of ~30% of the voyage duration. For a 5-day (120-hour) voyage, this is ~36 hours, consistent with the characteristic lifetime of mid-latitude synoptic systems (1–3 days).

5.2 Cholesky Decomposition

The correlation matrix C is decomposed into its lower-triangular Cholesky factor:

Correlated standard normal variates are then generated as:

This is computationally efficient: the Cholesky decomposition (O(n³)) is computed once for the entire batch of simulations, and each simulation only requires a matrix-vector multiplication (O(n²)), which for n = 100 takes microseconds.

5.3 Physical Interpretation

The exponential kernel produces perturbation sequences where:

• Adjacent slices (1.2h apart) have correlation ≈ 0.96 — nearly identical perturbation
• Slices 12h apart have correlation ≈ 0.67 — similar but diverging
• Slices 36h apart have correlation ≈ 0.37 — moderately correlated
• Slices 72h apart have correlation ≈ 0.14 — nearly independent

This means that if simulation #42 experiences a wind perturbation factor of 1.4 at hour 24, it will likely experience ~1.35 at hour 27 and ~1.2 at hour 36, but might experience 0.9 at hour 96. This temporal coherence is critical for producing realistic fuel consumption distributions — a simulation where wind is alternately high and low every few hours would underestimate variance.

6. Perturbation Model

Weather parameters are perturbed using multiplicative factors drawn from log-normal distributions (for positive quantities like speed and height) and additive Gaussian offsets (for directions).

6.1 Log-Normal Multiplicative Factors

For wind speed, wave height, and current speed, the perturbation factor is:

6.2 Parameter-Specific Sigmas

Wind has the largest sigma (σ = 0.35) because GFS forecast errors grow significantly with lead time, and wind resistance is the dominant weather penalty for laden tankers (Blendermann wind resistance model). Wave height has a smaller sigma (σ = 0.20) because wave models assimilate altimetry data and exhibit lower relative errors than wind forecasts at the same lead time.

6.3 Direction Perturbations

Wind and wave directions are perturbed by additive Gaussian offsets with σ = 15°, temporally correlated using the same Cholesky mechanism:

7. Wind–Wave Coupling

In the real ocean, waves are driven by wind. An increase in wind speed leads to increased wave height after a development time that depends on fetch and duration. The perturbation model captures this physical relationship by making the wave height perturbation partially correlated with the wind speed perturbation.

The correlation coefficient ρ = 0.7 means that 49% of wave height variance is explained by wind speed variance, with the remaining 51% from independent factors (swell from remote storms, local bathymetric effects, wave model biases). This value is consistent with observational studies of the wind–sea relationship in open ocean conditions (Hasselmann et al., 1973; Young, 1999).

The coupling ensures that simulations where the vessel encounters higher-than-forecast winds also tend to experience higher-than-forecast waves, producing realistic combined resistance penalties. Without this coupling, the Monte Carlo would underestimate the variance of fuel consumption because high-wind and high-wave events would occur independently instead of together.

8. Algorithm

The complete Monte Carlo simulation proceeds as follows:

8.1 Pre-computation (once per batch)

1. Discretize route into n time slices with timestamps and positions
2. Pre-fetch base weather for all slices (DB wave grids + GFS wind + CMEMS current)
3. Compute Cholesky factor L of the n × n exponential correlation matrix

8.2 Per-simulation loop (N = 100 times)

1. Generate 4 independent standard normal vectors: ε_wind, ε_wave,indep, ε_current, ε_dir
2. Multiply each by L to obtain temporally correlated normals: z = L · ε
3. Apply wind–wave coupling: z_wave = 0.7 z_wind + 0.714 z_wave,indep
4. Convert to perturbation factors: f = exp(σ z - ½σ²)
5. Build perturbed weather provider callable:
   • For each query (lat, lon, time), find the nearest time slice (top-5 by time, then nearest spatially)
   • Return base weather at that slice, multiplied by the slice's perturbation factors
6. Run full voyage calculation with the perturbed weather provider
7. Collect total_time, total_fuel, arrival_time

8.3 Post-processing

1. Compute 10th, 50th, 90th percentiles of fuel, time, and arrival arrays
2. Return MonteCarloResult with P10/P50/P90 confidence intervals

8.4 Weather Provider Matching

The perturbed weather provider maps each voyage calculator query to the nearest pre-computed time slice using a two-stage nearest-neighbour search:

1. Find the 5 slices closest in time (using pre-computed timestamps as a numpy array)
2. Among those 5, select the one closest spatially (Euclidean distance in lat/lon space)

This hybrid time-space matching handles the case where the voyage calculator queries leg midpoints that don't exactly coincide with the pre-computed slice positions.

9. Performance

9.1 Computational Cost Breakdown

9.2 Scaling

The simulation scales linearly with the number of simulations and linearly with the number of route legs. The weather pre-fetch is amortized across all simulations and is the dominant fixed cost. Doubling the number of simulations from 100 to 200 adds ~250ms.

9.3 Comparison with Live Weather

Without the pre-ingested weather database, each simulation would need to query weather at each slice position via the unified weather provider, which triggers live CMEMS/GFS API calls. In testing, this took 78.5 seconds for 50 simulations (~1.6s per simulation, dominated by HTTP latency). The database pre-fetch architecture achieved a 200× speedup.

10. Limitations and Future Work

10.1 Current Limitations

• No spatial correlation — Perturbations are correlated in time but not in space. A proper spatial correlation model would require 2D random fields (e.g., Matérn covariance), significantly increasing computational cost.
• Isotropic perturbations — The same perturbation magnitude is applied regardless of geographic region. In practice, forecast skill varies by region (better in data-rich areas, worse in the Southern Ocean).
• Fixed sigma values — The perturbation sigmas are constant across forecast lead times. In reality, forecast error grows with lead time; an adaptive sigma schedule would improve realism.
• Single forecast scenario — The base weather comes from a single deterministic forecast (best available). Using actual NWP ensemble members would capture physically consistent large-scale pattern uncertainty that parametric perturbations cannot.
• Time-invariant ETA in mild weather — The voyage calculator only reduces speed when required power exceeds 90% MCR. In mild weather (e.g., Mediterranean summer), engine power is sufficient to maintain commanded speed, producing zero ETA variance. Only fuel consumption varies.

10.2 Future Improvements

• Lead-time-dependent sigma — Scale perturbation magnitude with forecast lead time: σ(t) = σ₀ · (1 + α · t/T), calibrated against GFS/CMEMS verification statistics.
• NWP ensemble integration — When available, use GEFS or ECMWF ensemble members as base weather scenarios instead of parametric perturbations.
• Spatial correlation via Matérn fields — Generate spatially correlated perturbation fields for routes that span multiple weather regimes.
• Adaptive grid resolution — Increase slice density in regions of high weather gradient (e.g., near frontal zones) and decrease in stable regions.

References

1. Aijjou, A.E., Bahaj, M., and Wakrime, A. (2021). “A Comprehensive Approach to Account for Weather Uncertainties in Ship Route Optimization.” Journal of Marine Science and Engineering, 9(12):1434. DOI: 10.3390/jmse9121434
2. Dickson, T., Farr, H., Sear, D., and Blake, J. (2019). “Uncertainty in marine weather routing.” arXiv preprint, arXiv:1901.03840.
3. Hamill, T.M., Whitaker, J.S., and Wei, X. (2004). “Ensemble Reforecasting: Improving Medium-Range Forecast Skill Using Retrospective Forecasts.” Monthly Weather Review, 132(6):1434–1447.
4. Hanssen, G.L. and James, R.W. (1960). “Optimum ship routing.” Journal of Navigation, 13(3):253–272.
5. Hasselmann, K. et al. (1973). “Measurements of wind-wave growth and swell decay during the Joint North Sea Wave Project (JONSWAP).” Ergänzungsheft zur Deutschen Hydrographischen Zeitschrift, Reihe A, Nr. 12.
6. Holton, J.R. and Hakim, G.J. (2013). An Introduction to Dynamic Meteorology, 5th edition. Academic Press.
7. Mannarini, G. and Carelli, L. (2019). “VISIR-1.b: Ocean surface gravity waves and currents for energy-efficient navigation.” Geoscientific Model Development, 12:3449–3480.
8. Pinson, P. and Madsen, H. (2012). “Adaptive modelling and forecasting of offshore wind power fluctuations with Markov-switching autoregressive models.” Journal of Forecasting, 31(4):281–313.
9. Vettor, R. and Guedes Soares, C. (2016). “Development of a ship weather routing system.” Ocean Engineering, 123:1–14.
10. Young, I.R. (1999). Wind Generated Ocean Waves. Elsevier Ocean Engineering Book Series, Vol. 2.

Further Reading

1. Nielsen, U.D. et al. “Wave Estimation and Forecasting Using Ships as Wave Buoys.” DTU Mechanical Engineering, Technical University of Denmark. orbit.dtu.dk
2. Heikkilä, M. et al. “Improving the Global Predictions of Shipping Emission by Modelling the Effects of Ambient Conditions.” In: Air Pollution Modeling and its Application XXIX, Springer, 2025. DOI: 10.1007/978-3-031-70424-6_34
3. ITTC (2021). “Recommended Procedures and Guidelines 7.5–02–02–04: Wave Profile Measurement and Wave Pattern Resistance Analysis.” International Towing Tank Conference. ittc.info