Risk Assessment Skill

name: risk-assessment
version: 1.0.0
category: sme
tags: [risk, probabilistic, monte-carlo, reliability, uncertainty, sensitivity-analysis, marine-safety]
created: 2026-01-06
updated: 2026-01-06
author: Claude
description: |
  Expert risk assessment and probabilistic analysis for marine and offshore
  operations. Includes Monte Carlo simulations, reliability calculations,
  uncertainty quantification, and decision-making under risk.

When to Use This Skill

Use this skill when you need to:

• Perform Monte Carlo simulations for uncertainty quantification
• Calculate system reliability and failure probabilities
• Conduct sensitivity analysis to identify critical parameters
• Create risk matrices for hazard assessment
• Perform probabilistic design and analysis
• Quantify uncertainties in marine operations
• Make decisions under uncertainty with risk metrics
• Validate designs against reliability targets

Core Knowledge Areas

1. Monte Carlo Simulation

Basic Monte Carlo framework:

import numpy as np
from scipy import stats
from dataclasses import dataclass
from typing import Callable, Dict, List, Tuple, Optional
import pandas as pd

@dataclass
class RandomVariable:
    """Statistical distribution for a random variable."""
    name: str
    distribution: str  # 'normal', 'lognormal', 'uniform', 'weibull', etc.
    parameters: dict  # Distribution-specific parameters

    def sample(self, size: int = 1) -> np.ndarray:
        """
        Generate random samples from distribution.

        Args:
            size: Number of samples

        Returns:
            Array of random samples

        Example:
            >>> rv = RandomVariable(
            ...     name='wave_height',
            ...     distribution='weibull',
            ...     parameters={'c': 2.0, 'scale': 3.5}
            ... )
            >>> samples = rv.sample(1000)
        """
        if self.distribution == 'normal':
            return np.random.normal(
                loc=self.parameters['mean'],
                scale=self.parameters['std'],
                size=size
            )
        elif self.distribution == 'lognormal':
            return np.random.lognormal(
                mean=self.parameters['mean'],
                sigma=self.parameters['std'],
                size=size
            )
        elif self.distribution == 'uniform':
            return np.random.uniform(
                low=self.parameters['min'],
                high=self.parameters['max'],
                size=size
            )
        elif self.distribution == 'weibull':
            return stats.weibull_min.rvs(
                c=self.parameters['c'],
                scale=self.parameters['scale'],
                size=size
            )
        elif self.distribution == 'exponential':
            return np.random.exponential(
                scale=self.parameters['scale'],
                size=size
            )
        else:
            raise ValueError(f"Unknown distribution: {self.distribution}")

def monte_carlo_simulation(
    model: Callable,
    random_variables: List[RandomVariable],
    n_samples: int = 10000,
    seed: Optional[int] = None
) -> Dict[str, np.ndarray]:
    """
    Perform Monte Carlo simulation.

    Args:
        model: Function that takes random variable samples and returns result
        random_variables: List of RandomVariable objects
        n_samples: Number of Monte Carlo samples
        seed: Random seed for reproducibility

    Returns:
        Dictionary with input samples and output results

    Example:
        >>> def mooring_tension_model(wave_height, current_speed):
        ...     # Simplified mooring tension model
        ...     return 1000 + 150 * wave_height + 50 * current_speed
        >>>
        >>> rv_wave = RandomVariable(
        ...     'wave_height',
        ...     'weibull',
        ...     {'c': 2.0, 'scale': 3.5}
        ... )
        >>> rv_current = RandomVariable(
        ...     'current_speed',
        ...     'normal',
        ...     {'mean': 1.0, 'std': 0.3}
        ... )
        >>>
        >>> results = monte_carlo_simulation(
        ...     model=lambda samples: mooring_tension_model(
        ...         samples['wave_height'],
        ...         samples['current_speed']
        ...     ),
        ...     random_variables=[rv_wave, rv_current],
        ...     n_samples=10000
        ... )
    """
    if seed is not None:
        np.random.seed(seed)

    # Generate samples for all random variables
    samples = {}
    for rv in random_variables:
        samples[rv.name] = rv.sample(n_samples)

    # Run model for all samples
    outputs = model(samples)

    # Combine inputs and outputs
    results = {**samples, 'output': outputs}

    return results

def calculate_statistics(
    data: np.ndarray,
    percentiles: List[float] = [5, 50, 95]
) -> dict:
    """
    Calculate statistical parameters from Monte Carlo results.

    Args:
        data: Array of simulation results
        percentiles: Percentile values to calculate

    Returns:
        Dictionary with statistics

    Example:
        >>> stats = calculate_statistics(results['output'])
        >>> print(f"Mean: {stats['mean']:.1f}")
        >>> print(f"Std: {stats['std']:.1f}")
        >>> print(f"95th percentile: {stats['p95']:.1f}")
    """
    stats_dict = {
        'mean': np.mean(data),
        'std': np.std(data),
        'min': np.min(data),
        'max': np.max(data),
        'median': np.median(data),
        'cv': np.std(data) / np.mean(data)  # Coefficient of variation
    }

    # Add percentiles
    for p in percentiles:
        stats_dict[f'p{p}'] = np.percentile(data, p)

    return stats_dict

2. Reliability Analysis

Calculate reliability and failure probability:

def calculate_reliability(
    response_data: np.ndarray,
    limit_state: float,
    mode: str = 'less_than'
) -> dict:
    """
    Calculate reliability from Monte Carlo results.

    Args:
        response_data: Array of response values
        limit_state: Limit state threshold
        mode: 'less_than' or 'greater_than'

    Returns:
        Dictionary with reliability metrics

    Example:
        >>> # Mooring tension should be < 8000 kN
        >>> reliability = calculate_reliability(
        ...     response_data=results['output'],
        ...     limit_state=8000,
        ...     mode='less_than'
        ... )
        >>> print(f"Reliability: {reliability['reliability']:.4f}")
        >>> print(f"Probability of failure: {reliability['pf']:.6f}")
    """
    n_samples = len(response_data)

    if mode == 'less_than':
        # Failure when response >= limit_state
        failures = response_data >= limit_state
    elif mode == 'greater_than':
        # Failure when response <= limit_state
        failures = response_data <= limit_state
    else:
        raise ValueError(f"Unknown mode: {mode}")

    n_failures = np.sum(failures)
    pf = n_failures / n_samples  # Probability of failure
    reliability = 1 - pf

    # Reliability index (beta)
    # β = -Φ^(-1)(Pf) where Φ is standard normal CDF
    if pf > 0 and pf < 1:
        beta = -stats.norm.ppf(pf)
    elif pf == 0:
        beta = np.inf
    else:  # pf == 1
        beta = -np.inf

    return {
        'n_samples': n_samples,
        'n_failures': n_failures,
        'pf': pf,
        'reliability': reliability,
        'beta': beta,
        'target_met': reliability >= 0.9999  # Example: 4-9's reliability
    }

def system_reliability_series(
    component_reliabilities: List[float]
) -> float:
    """
    Calculate system reliability for series system.

    Series system: All components must work for system to work.

    Args:
        component_reliabilities: List of component reliabilities

    Returns:
        System reliability

    Example:
        >>> # 3 mooring lines in series (all must hold)
        >>> R_sys = system_reliability_series([0.999, 0.9995, 0.998])
        >>> print(f"System reliability: {R_sys:.6f}")
    """
    # R_sys = Π R_i (product of all component reliabilities)
    R_sys = np.prod(component_reliabilities)
    return R_sys

def system_reliability_parallel(
    component_reliabilities: List[float]
) -> float:
    """
    Calculate system reliability for parallel system.

    Parallel system: At least one component must work for system to work.

    Args:
        component_reliabilities: List of component reliabilities

    Returns:
        System reliability

    Example:
        >>> # 3 redundant mooring lines (only 1 needs to hold)
        >>> R_sys = system_reliability_parallel([0.95, 0.95, 0.95])
        >>> print(f"System reliability: {R_sys:.6f}")
    """
    # R_sys = 1 - Π(1 - R_i) (1 minus product of all failures)
    pf_components = [1 - R for R in component_reliabilities]
    pf_sys = np.prod(pf_components)
    R_sys = 1 - pf_sys
    return R_sys

def calculate_form_reliability(
    mean: np.ndarray,
    cov_matrix: np.ndarray,
    limit_state_gradient: np.ndarray
) -> dict:
    """
    First Order Reliability Method (FORM) for reliability analysis.

    Args:
        mean: Mean values of random variables
        cov_matrix: Covariance matrix
        limit_state_gradient: Gradient of limit state function at mean point

    Returns:
        Dictionary with FORM results

    Example:
        >>> # 2 random variables: Hs, Tp
        >>> mean = np.array([5.0, 10.0])
        >>> cov = np.array([[1.0, 0.5], [0.5, 2.0]])
        >>> gradient = np.array([150.0, 50.0])  # ∂g/∂Hs, ∂g/∂Tp
        >>> form_result = calculate_form_reliability(mean, cov, gradient)
        >>> print(f"Reliability index β: {form_result['beta']:.3f}")
    """
    # Reliability index: β = μ / sqrt(∇g^T Σ ∇g)
    # where μ is mean, Σ is covariance, ∇g is gradient

    numerator = np.dot(gradient, mean)
    denominator = np.sqrt(
        np.dot(gradient, np.dot(cov_matrix, gradient))
    )

    beta = numerator / denominator

    # Probability of failure from β
    pf = stats.norm.cdf(-beta)

    return {
        'beta': beta,
        'pf': pf,
        'reliability': 1 - pf
    }

3. Sensitivity Analysis

Identify critical parameters:

def sensitivity_analysis_correlation(
    inputs: Dict[str, np.ndarray],
    output: np.ndarray
) -> pd.DataFrame:
    """
    Sensitivity analysis using correlation coefficients.

    Args:
        inputs: Dictionary of input variable samples
        output: Output variable samples

    Returns:
        DataFrame with correlation coefficients sorted by absolute value

    Example:
        >>> sensitivity = sensitivity_analysis_correlation(
        ...     inputs={
        ...         'wave_height': results['wave_height'],
        ...         'current_speed': results['current_speed'],
        ...         'wind_speed': results['wind_speed']
        ...     },
        ...     output=results['output']
        ... )
        >>> print(sensitivity)
    """
    correlations = {}

    for var_name, var_samples in inputs.items():
        # Pearson correlation coefficient
        corr = np.corrcoef(var_samples, output)[0, 1]
        correlations[var_name] = corr

    # Create DataFrame
    df = pd.DataFrame.from_dict(
        correlations,
        orient='index',
        columns=['correlation']
    )

    # Add absolute value and rank
    df['abs_correlation'] = df['correlation'].abs()
    df = df.sort_values('abs_correlation', ascending=False)
    df['rank'] = range(1, len(df) + 1)

    return df

def sensitivity_analysis_variance(
    inputs: Dict[str, np.ndarray],
    output: np.ndarray,
    n_partitions: int = 10
) -> pd.DataFrame:
    """
    Variance-based sensitivity analysis (Sobol indices approximation).

    Args:
        inputs: Dictionary of input variable samples
        output: Output variable samples
        n_partitions: Number of partitions for conditional variance

    Returns:
        DataFrame with sensitivity indices

    Example:
        >>> sensitivity = sensitivity_analysis_variance(
        ...     inputs=results,
        ...     output=results['output']
        ... )
    """
    total_variance = np.var(output)

    sensitivity_indices = {}

    for var_name, var_samples in inputs.items():
        if var_name == 'output':
            continue

        # Partition variable range
        percentiles = np.linspace(0, 100, n_partitions + 1)
        bins = np.percentile(var_samples, percentiles)

        # Calculate conditional variance
        conditional_means = []
        for i in range(n_partitions):
            mask = (var_samples >= bins[i]) & (var_samples < bins[i + 1])
            if np.sum(mask) > 0:
                conditional_means.append(np.mean(output[mask]))

        # Variance of conditional means
        variance_of_means = np.var(conditional_means)

        # First-order sensitivity index (approximation)
        S1 = variance_of_means / total_variance if total_variance > 0 else 0

        sensitivity_indices[var_name] = {
            'first_order': S1
        }

    # Create DataFrame
    df = pd.DataFrame.from_dict(sensitivity_indices, orient='index')
    df = df.sort_values('first_order', ascending=False)
    df['rank'] = range(1, len(df) + 1)

    return df

def tornado_plot_data(
    base_value: float,
    variables: Dict[str, dict],
    model: Callable
) -> pd.DataFrame:
    """
    Generate data for tornado diagram (one-at-a-time sensitivity).

    Args:
        base_value: Base case output value
        variables: Dict with variable names and {'low', 'high', 'nominal'} values
        model: Model function

    Returns:
        DataFrame with tornado plot data

    Example:
        >>> variables = {
        ...     'wave_height': {'low': 3.0, 'nominal': 5.0, 'high': 8.0},
        ...     'current_speed': {'low': 0.5, 'nominal': 1.0, 'high': 1.5}
        ... }
        >>> tornado_data = tornado_plot_data(
        ...     base_value=5000,
        ...     variables=variables,
        ...     model=mooring_model
        ... )
    """
    results = []

    for var_name, var_values in variables.items():
        # Calculate output at low and high values
        # (keeping other variables at nominal)

        inputs_low = {k: v['nominal'] for k, v in variables.items()}
        inputs_low[var_name] = var_values['low']
        output_low = model(inputs_low)

        inputs_high = {k: v['nominal'] for k, v in variables.items()}
        inputs_high[var_name] = var_values['high']
        output_high = model(inputs_high)

        # Calculate swing
        swing = abs(output_high - output_low)

        results.append({
            'variable': var_name,
            'output_low': output_low,
            'output_high': output_high,
            'swing': swing,
            'low_value': var_values['low'],
            'high_value': var_values['high']
        })

    df = pd.DataFrame(results)
    df = df.sort_values('swing', ascending=False)
    df['rank'] = range(1, len(df) + 1)

    return df

4. Risk Matrices and Hazard Assessment

from enum import Enum

class Severity(Enum):
    """Consequence severity levels."""
    NEGLIGIBLE = 1
    MINOR = 2
    MODERATE = 3
    MAJOR = 4
    CATASTROPHIC = 5

class Likelihood(Enum):
    """Event likelihood levels."""
    RARE = 1
    UNLIKELY = 2
    POSSIBLE = 3
    LIKELY = 4
    ALMOST_CERTAIN = 5

class RiskLevel(Enum):
    """Risk level categories."""
    LOW = 1
    MEDIUM = 2
    HIGH = 3
    VERY_HIGH = 4

@dataclass
class Hazard:
    """Hazard definition."""
    id: str
    description: str
    severity: Severity
    likelihood: Likelihood
    existing_controls: List[str]
    residual_risk: Optional[RiskLevel] = None

def calculate_risk_level(
    severity: Severity,
    likelihood: Likelihood
) -> RiskLevel:
    """
    Calculate risk level from severity and likelihood.

    Uses 5x5 risk matrix.

    Args:
        severity: Consequence severity
        likelihood: Event likelihood

    Returns:
        Risk level

    Example:
        >>> risk = calculate_risk_level(
        ...     Severity.MAJOR,
        ...     Likelihood.POSSIBLE
        ... )
        >>> print(risk)  # RiskLevel.HIGH
    """
    # 5x5 risk matrix
    # Rows: Likelihood (1-5)
    # Columns: Severity (1-5)
    risk_matrix = np.array([
        [1, 1, 2, 2, 3],  # Rare
        [1, 2, 2, 3, 3],  # Unlikely
        [2, 2, 3, 3, 4],  # Possible
        [2, 3, 3, 4, 4],  # Likely
        [3, 3, 4, 4, 4]   # Almost Certain
    ])

    risk_value = risk_matrix[likelihood.value - 1, severity.value - 1]

    # Map to RiskLevel
    if risk_value == 1:
        return RiskLevel.LOW
    elif risk_value == 2:
        return RiskLevel.MEDIUM
    elif risk_value == 3:
        return RiskLevel.HIGH
    else:  # risk_value == 4
        return RiskLevel.VERY_HIGH

def assess_hazards(
    hazards: List[Hazard]
) -> pd.DataFrame:
    """
    Assess list of hazards and calculate risk levels.

    Args:
        hazards: List of Hazard objects

    Returns:
        DataFrame with hazard assessment

    Example:
        >>> hazards = [
        ...     Hazard(
        ...         id='H001',
        ...         description='Mooring line failure',
        ...         severity=Severity.MAJOR,
        ...         likelihood=Likelihood.UNLIKELY,
        ...         existing_controls=['Regular inspection', 'Design factor 2.5']
        ...     ),
        ...     Hazard(
        ...         id='H002',
        ...         description='Vessel collision',
        ...         severity=Severity.CATASTROPHIC,
        ...         likelihood=Likelihood.RARE,
        ...         existing_controls=['AIS', 'Exclusion zone', 'Radar']
        ...     )
        ... ]
        >>> assessment = assess_hazards(hazards)
    """
    results = []

    for hazard in hazards:
        risk_level = calculate_risk_level(hazard.severity, hazard.likelihood)

        results.append({
            'hazard_id': hazard.id,
            'description': hazard.description,
            'severity': hazard.severity.name,
            'likelihood': hazard.likelihood.name,
            'risk_level': risk_level.name,
            'controls': '; '.join(hazard.existing_controls)
        })

    df = pd.DataFrame(results)

    # Sort by risk level (descending)
    risk_order = {
        'VERY_HIGH': 4,
        'HIGH': 3,
        'MEDIUM': 2,
        'LOW': 1
    }
    df['risk_value'] = df['risk_level'].map(risk_order)
    df = df.sort_values('risk_value', ascending=False)
    df = df.drop('risk_value', axis=1)

    return df

5. Extreme Value Analysis

from scipy.stats import genextreme

def fit_extreme_value_distribution(
    data: np.ndarray,
    method: str = 'gev'
) -> dict:
    """
    Fit extreme value distribution to data.

    Args:
        data: Array of extreme values (e.g., annual maxima)
        method: 'gev' (Generalized Extreme Value) or 'gumbel'

    Returns:
        Dictionary with fitted parameters and statistics

    Example:
        >>> # Annual maximum wave heights
        >>> annual_max_Hs = np.array([8.5, 9.2, 7.8, 10.1, 8.9, ...])
        >>> fit = fit_extreme_value_distribution(annual_max_Hs, 'gev')
        >>> print(f"Shape: {fit['shape']:.3f}")
        >>> print(f"Location: {fit['location']:.2f}")
        >>> print(f"Scale: {fit['scale']:.2f}")
    """
    if method == 'gev':
        # Fit GEV distribution
        shape, location, scale = genextreme.fit(data)

        # Calculate return values
        return_periods = [10, 50, 100, 1000, 10000]  # years
        return_values = {}

        for T in return_periods:
            # Return level for return period T
            p = 1 - 1/T
            x_T = genextreme.ppf(p, shape, loc=location, scale=scale)
            return_values[f'{T}yr'] = x_T

        return {
            'shape': shape,
            'location': location,
            'scale': scale,
            'return_values': return_values
        }

    elif method == 'gumbel':
        # Gumbel is special case of GEV with shape=0
        location, scale = stats.gumbel_r.fit(data)

        return_periods = [10, 50, 100, 1000, 10000]
        return_values = {}

        for T in return_periods:
            p = 1 - 1/T
            x_T = stats.gumbel_r.ppf(p, loc=location, scale=scale)
            return_values[f'{T}yr'] = x_T

        return {
            'shape': 0.0,  # Gumbel has shape = 0
            'location': location,
            'scale': scale,
            'return_values': return_values
        }

    else:
        raise ValueError(f"Unknown method: {method}")

Complete Examples

Example 1: Complete Mooring System Risk Assessment

import numpy as np
from pathlib import Path

def complete_mooring_risk_assessment(
    design_parameters: dict,
    environmental_parameters: dict,
    n_simulations: int = 10000
) -> dict:
    """
    Complete probabilistic risk assessment for mooring system.

    Example:
        >>> design = {
        ...     'n_lines': 8,
        ...     'line_capacity': 10000,  # kN
        ...     'safety_factor': 2.5
        ... }
        >>> environment = {
        ...     'wave_height': {'distribution': 'weibull', 'c': 2.0, 'scale': 3.5},
        ...     'wave_period': {'distribution': 'normal', 'mean': 10.0, 'std': 2.0},
        ...     'current_speed': {'distribution': 'normal', 'mean': 1.0, 'std': 0.3}
        ... }
        >>> results = complete_mooring_risk_assessment(
        ...     design,
        ...     environment,
        ...     n_simulations=10000
        ... )
    """
    print("="*70)
    print("MOORING SYSTEM RISK ASSESSMENT")
    print("="*70)

    # Step 1: Define random variables
    print("\n[Step 1/6] Defining random variables...")

    random_vars = []
    for var_name, var_params in environment_parameters.items():
        rv = RandomVariable(
            name=var_name,
            distribution=var_params['distribution'],
            parameters={k: v for k, v in var_params.items() if k != 'distribution'}
        )
        random_vars.append(rv)
        print(f"  {var_name}: {var_params['distribution']}")

    # Step 2: Define mooring tension model
    print("\n[Step 2/6] Defining mooring tension model...")

    def mooring_tension_model(samples):
        """Simplified mooring tension model."""
        Hs = samples['wave_height']
        Tp = samples['wave_period']
        V_c = samples['current_speed']

        # Simplified empirical model
        tension = (
            1000 +  # Pretension
            200 * Hs +  # Wave contribution
            50 * Hs**2 / Tp +  # Dynamic amplification
            100 * V_c  # Current drag
        )

        return tension

    # Step 3: Run Monte Carlo simulation
    print(f"\n[Step 3/6] Running Monte Carlo simulation ({n_simulations} samples)...")

    mc_results = monte_carlo_simulation(
        model=mooring_tension_model,
        random_variables=random_vars,
        n_samples=n_simulations,
        seed=42
    )

    # Calculate statistics
    stats = calculate_statistics(mc_results['output'])

    print(f"  Mean tension: {stats['mean']:.1f} kN")
    print(f"  Std deviation: {stats['std']:.1f} kN")
    print(f"  95th percentile: {stats['p95']:.1f} kN")
    print(f"  Max: {stats['max']:.1f} kN")

    # Step 4: Reliability analysis
    print("\n[Step 4/6] Performing reliability analysis...")

    # Design limit (capacity / safety factor)
    design_limit = design_parameters['line_capacity'] / design_parameters['safety_factor']

    reliability_result = calculate_reliability(
        response_data=mc_results['output'],
        limit_state=design_limit,
        mode='less_than'
    )

    print(f"  Design limit: {design_limit:.1f} kN")
    print(f"  Probability of failure: {reliability_result['pf']:.6f}")
    print(f"  Reliability: {reliability_result['reliability']:.6f}")
    print(f"  Reliability index β: {reliability_result['beta']:.3f}")

    # Step 5: Sensitivity analysis
    print("\n[Step 5/6] Performing sensitivity analysis...")

    sensitivity = sensitivity_analysis_correlation(
        inputs={k: v for k, v in mc_results.items() if k != 'output'},
        output=mc_results['output']
    )

    print(f"\nSensitivity ranking:")
    for idx, row in sensitivity.iterrows():
        print(f"  {idx}: correlation = {row['correlation']:.3f}")

    # Step 6: System reliability
    print("\n[Step 6/6] Calculating system reliability...")

    # Assume independent mooring lines
    component_reliability = reliability_result['reliability']

    # Series system (all lines must hold for 100% intact)
    R_all_intact = system_reliability_series(
        [component_reliability] * design_parameters['n_lines']
    )

    # Parallel redundancy (system fails if all lines fail)
    R_total_failure = 1 - (1 - component_reliability)**design_parameters['n_lines']

    print(f"  Single line reliability: {component_reliability:.6f}")
    print(f"  All lines intact: {R_all_intact:.6f}")
    print(f"  Avoid total failure: {R_total_failure:.10f}")

    # Summary
    summary = {
        'statistics': stats,
        'reliability': reliability_result,
        'sensitivity': sensitivity,
        'system_reliability': {
            'all_intact': R_all_intact,
            'total_failure_avoidance': R_total_failure
        }
    }

    print("\n" + "="*70)
    print("ASSESSMENT COMPLETE")
    print("="*70)

    return summary

# Run assessment
design_parameters = {
    'n_lines': 8,
    'line_capacity': 10000,
    'safety_factor': 2.5
}

environmental_parameters = {
    'wave_height': {'distribution': 'weibull', 'c': 2.0, 'scale': 3.5},
    'wave_period': {'distribution': 'normal', 'mean': 10.0, 'std': 2.0},
    'current_speed': {'distribution': 'normal', 'mean': 1.0, 'std': 0.3}
}

risk_assessment = complete_mooring_risk_assessment(
    design_parameters,
    environmental_parameters,
    n_simulations=10000
)

Well Planning Risk Context

The Actionability Gap

Standard 5×5 risk matrices capture what the risk is and how severe it is, but not who has the authority to mitigate it. In well planning, this creates a systematic bias:

• Engineers write risks they can act on (operational: mud weight, BOP procedures, casing design)
• Engineers omit risks they understand but cannot control (strategic: downhole tool contracts, rig selection, fleet capacity)
• The risk register becomes operationally complete but strategically incomplete

Risk Authority Tiers

When building risk registers for well planning, classify every risk into one of three tiers:

Structured Escalation for Strategic Risks

When a risk falls outside project authority, convert it from a vague register entry into a structured escalation:

Risk: [Formation X] requires [capability Y] that current [tool/rig/contract] does not provide.
Impact if unaddressed: [specific consequence — NPT days, sidetrack probability, well objective at risk]
Mitigation lever: [what needs to change — new tool contract, different rig class, regulatory exemption]
Decision owner: [organizational function — procurement, fleet management, regulatory affairs]
Decision timeline: [when the decision must be made relative to spud date]

This reframes “we need a better downhole tool” into an actionable request with a clear owner, deadline, and consequence.

Risk Influence Map

Categorize the project’s relationship to each risk:

• Controls: Project team can directly mitigate (operational risks)
• Influences: Project team provides input but doesn’t decide (tactical risks)
• Observes: Project team identifies and escalates but has no lever (strategic risks)

A well planning risk register should have risks in all three categories. If the register contains only operational risks, the strategic dimension is being suppressed — not because those risks don’t exist, but because the project lacks empowerment to address them.

Integration with Quantitative Assessment

The Monte Carlo and reliability tools in this skill quantify operational risks effectively. For strategic risks, the quantitative output serves a different purpose — it provides the business case for escalation:

# Example: Quantify cost impact of rig capability gap
# This output goes into the escalation template, not the project risk register
rig_gap_risk = assess_hazards([
    Hazard(
        id='STR-001',
        description='Rig hookload insufficient for 9-5/8" casing at TD',
        severity=Severity.MAJOR,          # sidetrack or redesign required
        likelihood=Likelihood.POSSIBLE,    # depends on formation pressure
        existing_controls=['Lighter mud system (reduces margin)'],
    )
])

See WRK-163 for the full well planning risk empowerment framework implementation plan.

Best Practices

1. Sample Size Selection

def determine_sample_size(
    target_pf: float,
    confidence_level: float = 0.95
) -> int:
    """
    Determine required sample size for target failure probability.

    Rule of thumb: N ≈ 10/Pf for reasonable confidence

    Args:
        target_pf: Target probability of failure
        confidence_level: Confidence level

    Returns:
        Recommended sample size

    Example:
        >>> n = determine_sample_size(target_pf=1e-4, confidence_level=0.95)
        >>> print(f"Recommended samples: {n}")
    """
    # Basic rule: need at least 10 failures
    # So N * Pf ≥ 10
    # N ≥ 10 / Pf

    n_basic = int(10 / target_pf)

    # For higher confidence, increase further
    if confidence_level >= 0.99:
        n_recommended = n_basic * 5
    elif confidence_level >= 0.95:
        n_recommended = n_basic * 3
    else:
        n_recommended = n_basic * 2

    return max(n_recommended, 1000)  # Minimum 1000 samples

2. Convergence Checking

def check_monte_carlo_convergence(
    data: np.ndarray,
    window_size: int = 1000
) -> dict:
    """
    Check if Monte Carlo simulation has converged.

    Args:
        data: Simulation output data
        window_size: Window size for moving average

    Returns:
        Dictionary with convergence metrics
    """
    n = len(data)

    # Calculate cumulative mean
    cumulative_mean = np.cumsum(data) / np.arange(1, n + 1)

    # Calculate moving coefficient of variation
    if n > window_size:
        moving_std = np.std(data[-window_size:])
        moving_mean = np.mean(data[-window_size:])
        cov = moving_std / moving_mean if moving_mean != 0 else 0
    else:
        cov = np.std(data) / np.mean(data) if np.mean(data) != 0 else 0

    # Convergence criterion: COV < 5%
    converged = cov < 0.05

    return {
        'converged': converged,
        'cov': cov,
        'final_mean': cumulative_mean[-1],
        'samples_used': n
    }

Resources

Textbooks

• Ang, A.H-S., Tang, W.H. (2007). Probability Concepts in Engineering
• Melchers, R.E., Beck, A.T. (2018). Structural Reliability Analysis and Prediction
• DNV (2021). DNVGL-RP-C205: Environmental Conditions and Environmental Loads

Standards

• DNV-RP-C205: Environmental conditions and environmental loads
• DNV-RP-C206: Fatigue methodology of offshore ships
• ISO 2394: General principles on reliability for structures
• API RP 2A: Planning, designing and constructing fixed offshore platforms

Software

• @RISK: Monte Carlo simulation add-in for Excel
• Crystal Ball: Oracle’s risk analysis software
• OpenTURNS: Open source library for uncertainty quantification
• scipy.stats: Python statistical distributions

Use this skill for: Expert probabilistic risk assessment and reliability analysis for marine and offshore systems with comprehensive uncertainty quantification.