# Geological Simulation System

**Status**: Designed - Ready for Implementation
**Scope**: Complete planetary formation from accretion to industrial-ready geology
**Duration**: 4.65 billion years simulation in 95 cycles

## System Overview

Revolutionary geological simulation system that generates realistic planetary geology through scientifically-inspired processes. Creates diverse, coherent geology with proper resource distribution for industrial gameplay.

### Key Innovations
- **RegionalInfluence framework** adapted for geological processes
- **Tectonic regions** as simple circles with forces (not complex mesh)
- **Carbon region system** for realistic coal/oil formation
- **Dynamic sea level** affecting all geological processes
- **Optimized 24-byte tiles** supporting geological temperature ranges

## Simulation Phases

### Phase 1: Planetary Accretion (30 cycles × 100M years = 3.0 Ga)

**Initial Conditions**:
- Temperature: -100°C (space cold)
- Elevation: -30,000m (no surface)
- Empty planetary core

**Process - Meteorite Bombardment**:
```
For each wave (100M years):
    1. Generate 100-500 random meteorite impacts
    2. Each impact adds:
       - Kinetic energy → heat (up to +2000°C)
       - Mass → elevation change (+crater formation)
       - Metal composition → resource deposits
    3. Heavy metals sink to planetary core
    4. Core temperature drives volcanic eruptions
    5. Eruptions redistribute core metals to surface
    6. Gradual cooling between waves
```

**Expected Results**:
- Surface elevation: -30,000m → -15,000m
- Core composition: 80% iron, 15% nickel, 5% precious metals
- Crater-based geology with metal deposits
- Temperature stabilization around 1000-1500°C

### Phase 2: Tectonic Formation (25 cycles × 100M years = 2.5 Ga)

**Reduced Meteorite Activity**:
```cpp
// Phase 2 meteorite parameters (reduced from Phase 1)
impacts_per_wave = rng.randInt(10, 50);        // 10x reduction
meteorite_mass = rng.randFloat(10^12, 10^15);  // Smaller impacts
impact_probability = 0.3f;                     // 30% chance per cycle
```

**Tectonic Region System**:
```cpp
struct TectonicRegion {
    uint32_t region_id;
    float center_x, center_y;           // Mobile center (world coordinates)
    float radius;                       // Current radius (tiles)
    float velocity_x, velocity_y;       // Movement (tiles per cycle)
    float mass;                         // For repulsive force calculations
    RegionType type;                    // OCEANIC, CONTINENTAL, VOLCANIC

    // Evolution parameters
    float split_probability;            // Chance to divide (0.01 = 1% per cycle)
    float growth_rate;                  // Radius change per cycle (+/- tiles)
    float stability;                    // Resistance to external forces

    // Formation tracking
    int formation_cycle;                // When this region was created
    uint32_t parent_region_id;          // If split from another region
};

enum class RegionType {
    OCEANIC,        // Dense, subducts under continental
    CONTINENTAL,    // Light, forms mountains when compressed
    VOLCANIC,       // Active volcanic zone (temporary)
};
```

**Tectonic Physics Engine**:
```cpp
void updateTectonicRegions(float time_step) {
    // 1. Calculate repulsive forces between all regions
    for (auto& region1 : tectonic_regions) {
        Vector2 total_force = {0.0f, 0.0f};

        for (auto& region2 : tectonic_regions) {
            if (region1.region_id == region2.region_id) continue;

            float distance = calculateDistance(region1.center, region2.center);
            float overlap = (region1.radius + region2.radius) - distance;

            if (overlap > 0) {
                // Collision detected!
                Vector2 repulsion_direction = normalize(region1.center - region2.center);

                // Force magnitude: F = k * overlap^2 / (mass1 + mass2)
                float force_magnitude = REPULSION_CONSTANT * overlap * overlap /
                                      (region1.mass + region2.mass);

                Vector2 repulsion_force = repulsion_direction * force_magnitude;
                total_force += repulsion_force;

                // Create volcanic zone at collision point
                createVolcanicZone(region1, region2, overlap);
            }
        }

        // Apply acceleration: a = F/m
        Vector2 acceleration = total_force / region1.mass;
        region1.velocity_x += acceleration.x * time_step;
        region1.velocity_y += acceleration.y * time_step;

        // Apply velocity damping (friction)
        region1.velocity_x *= VELOCITY_DAMPING;
        region1.velocity_y *= VELOCITY_DAMPING;
    }

    // 2. Update positions and sizes
    for (auto& region : tectonic_regions) {
        // Move regions
        region.center_x += region.velocity_x * time_step;
        region.center_y += region.velocity_y * time_step;

        // Boundary wrapping (world is spherical)
        wrapCoordinates(region.center_x, region.center_y);

        // Evolve size
        region.radius += region.growth_rate * time_step;
        region.radius = std::clamp(region.radius, MIN_REGION_RADIUS, MAX_REGION_RADIUS);

        // Check for splitting
        if (region.radius > SPLIT_THRESHOLD &&
            rng.probability(region.split_probability * time_step)) {
            splitTectonicRegion(region);
        }
    }
}
```

**Volcanic Zone Creation**:
```cpp
void createVolcanicZone(TectonicRegion& r1, TectonicRegion& r2, float collision_intensity) {
    // Volcanic zone at intersection point
    Vector2 intersection = calculateIntersectionCenter(r1, r2);
    float volcanic_radius = collision_intensity * VOLCANIC_RADIUS_FACTOR;

    // Create temporary volcanic region
    TectonicRegion volcanic_zone = {
        .center_x = intersection.x,
        .center_y = intersection.y,
        .radius = volcanic_radius,
        .type = RegionType::VOLCANIC,
        .formation_cycle = current_cycle,
        .split_probability = 0.05f,  // High chance to break apart
        .growth_rate = -0.1f         // Shrinks over time
    };

    // Apply volcanic influence to affected tiles
    RegionalInfluence volcanic_influence = {
        .type = "volcanic_mountain_formation",
        .intensity = collision_intensity,
        .properties = json{
            {"elevation_bonus", collision_intensity * 500.0f},  // +500m per intensity
            {"temperature_bonus", collision_intensity * 300.0f}, // +300°C
            {"resource_deposits", json::array({"iron", "sulfur", "rare_metals"})},
            {"volcanic_features", json::array({"volcano", "hot_springs", "lava_tubes"})}
        }
    };

    applyRegionalInfluence(volcanic_influence, intersection.x, intersection.y, volcanic_radius);
    tectonic_regions.push_back(volcanic_zone);
}
```

**Region Splitting Algorithm**:
```cpp
void splitTectonicRegion(TectonicRegion& parent) {
    // Create two child regions
    float split_distance = parent.radius * 0.4f;
    Vector2 split_direction = randomUnitVector();

    TectonicRegion child1 = parent;
    TectonicRegion child2 = parent;

    // Assign new IDs
    child1.region_id = next_region_id++;
    child2.region_id = next_region_id++;
    child1.parent_region_id = parent.region_id;
    child2.parent_region_id = parent.region_id;

    // Separate positions
    child1.center_x += split_direction.x * split_distance;
    child1.center_y += split_direction.y * split_distance;
    child2.center_x -= split_direction.x * split_distance;
    child2.center_y -= split_direction.y * split_distance;

    // Reduce sizes
    child1.radius *= 0.7f;
    child2.radius *= 0.7f;

    // Opposite velocities (moving apart)
    child1.velocity_x = split_direction.x * SPLIT_VELOCITY;
    child1.velocity_y = split_direction.y * SPLIT_VELOCITY;
    child2.velocity_x = -split_direction.x * SPLIT_VELOCITY;
    child2.velocity_y = -split_direction.y * SPLIT_VELOCITY;

    // Replace parent with children
    removeRegion(parent);
    addRegion(child1);
    addRegion(child2);

    // Create rift valley between children
    createRiftValley(child1, child2);
}
```

**Expected Results**:
- Surface elevation: -15,000m → -5,000m (+10km crustal thickening)
- Formation of 15-25 stable tectonic regions
- Mountain ranges at collision zones (+1000-3000m elevation)
- Rift valleys where regions separate (-500m elevation)
- Distinct continental vs oceanic regions

### Phase 3: Hydrological Cycles (25 cycles × 20M years = 0.5 Ga)

**Dynamic Sea Level System**:
```cpp
struct OceanLevel {
    float current_sea_level;        // Global reference (meters)
    float ice_volume;               // Polar ice volume (km³)
    float thermal_expansion_factor; // Ocean thermal expansion coefficient
    float tectonic_basin_volume;    // Total oceanic basin volume

    void updateSeaLevel(float global_temperature, float volcanic_co2) {
        // 1. Ice volume changes (glaciation/melting)
        float target_ice_volume = calculateIceVolume(global_temperature);
        float ice_change = (target_ice_volume - ice_volume) * 0.1f; // 10% change per cycle
        ice_volume += ice_change;

        // Sea level change: 1 km³ ice = +2.7mm sea level
        float ice_effect = -ice_change * 0.0027f;

        // 2. Thermal expansion
        float thermal_effect = (global_temperature - 15.0f) * thermal_expansion_factor;

        // 3. Tectonic effects (basin formation/destruction)
        float tectonic_effect = (tectonic_basin_volume - baseline_basin_volume) * 0.001f;

        // 4. Update sea level
        current_sea_level += ice_effect + thermal_effect + tectonic_effect;
        current_sea_level = std::clamp(current_sea_level, -200.0f, +100.0f); // Realistic bounds
    }

    float calculateIceVolume(float global_temp) const {
        if (global_temp < -5.0f) return MAX_ICE_VOLUME;      // Ice age
        if (global_temp > 25.0f) return 0.0f;               // No ice
        return MAX_ICE_VOLUME * (25.0f - global_temp) / 30.0f; // Linear interpolation
    }
};
```

**Hydraulic Erosion System**:
```cpp
void applyHydraulicErosion(float cycle_duration_years) {
    const float HYDRAULIC_EFFICIENCY = 10.0f; // 10x more effective than tectonic

    for (int y = 0; y < map_height; ++y) {
        for (int x = 0; x < map_width; ++x) {
            WorldTile tile = world_map->getTile(x, y);
            float elevation = tile.getElevation();

            if (elevation > current_sea_level) {
                // CONTINENTAL EROSION

                // 1. Calculate water flow (slope-based)
                float average_neighbor_elevation = calculateAverageNeighborElevation(x, y);
                float slope = elevation - average_neighbor_elevation;
                float water_flow = std::max(0.0f, slope * 0.01f); // Flow intensity

                // 2. Erosion rate based on elevation and flow
                float elevation_factor = (elevation - current_sea_level) / 1000.0f; // Higher = more erosion
                float erosion_rate = water_flow * elevation_factor * HYDRAULIC_EFFICIENCY;

                // 3. Apply erosion
                float erosion_amount = erosion_rate * cycle_duration_years / 1000000.0f; // Per million years
                tile.setElevation(elevation - erosion_amount);

                // 4. Transport sediments downstream
                transportSediments(tile, x, y, erosion_amount);

            } else if (elevation > current_sea_level - 200.0f) {
                // COASTAL EROSION

                float depth_below_sea = current_sea_level - elevation;
                float coastal_erosion_rate = 50.0f; // Intense coastal erosion

                if (depth_below_sea < 50.0f) { // Shallow coastal zone
                    float erosion_amount = coastal_erosion_rate * cycle_duration_years / 1000000.0f;
                    tile.setElevation(elevation - erosion_amount);

                    // Form beaches and deltas
                    if (hasRiverInput(x, y)) {
                        formDelta(tile, x, y);
                    }
                }

            } else {
                // DEEP OCEAN - sediment deposition

                float sedimentation_rate = calculateSedimentInput(x, y);
                float deposition_amount = sedimentation_rate * cycle_duration_years / 1000000.0f;
                tile.setElevation(elevation + deposition_amount);
            }
        }
    }
}

void transportSediments(WorldTile& source_tile, int x, int y, float sediment_amount) {
    // Find steepest downhill direction
    auto [target_x, target_y] = findSteepestDescentDirection(x, y);

    if (target_x != x || target_y != y) {
        WorldTile target_tile = world_map->getTile(target_x, target_y);

        // Deposit sediments at target
        float current_elevation = target_tile.getElevation();
        target_tile.setElevation(current_elevation + sediment_amount * 0.5f); // 50% deposition

        // Create river network
        createRiverSegment(x, y, target_x, target_y);

        // Continue transport if target is above sea level
        if (target_tile.getElevation() > current_sea_level) {
            transportSediments(target_tile, target_x, target_y, sediment_amount * 0.5f);
        }
    }
}
```

**River Network Formation**:
```cpp
struct RiverNetwork {
    std::unordered_map<uint32_t, std::vector<uint32_t>> river_connections; // tile -> downstream tiles
    std::unordered_map<uint32_t, float> flow_volumes; // tile -> water volume

    void createRiverSegment(int from_x, int from_y, int to_x, int to_y) {
        uint32_t from_index = from_y * map_width + from_x;
        uint32_t to_index = to_y * map_width + to_x;

        // Add connection
        river_connections[from_index].push_back(to_index);

        // Increase flow volume
        flow_volumes[to_index] += flow_volumes[from_index] + 1.0f;

        // Mark tiles as having river features
        WorldTile from_tile = world_map->getTile(from_x, from_y);
        WorldTile to_tile = world_map->getTile(to_x, to_y);

        from_tile.setFlag(TileFlags::HAS_RIVER, true);
        to_tile.setFlag(TileFlags::HAS_RIVER, true);
    }

    void formDelta(WorldTile& tile, int x, int y) {
        // Delta formation where river meets ocean
        RegionalInfluence delta_influence = {
            .type = "river_delta",
            .intensity = flow_volumes[y * map_width + x] / 100.0f, // Intensity based on river size
            .properties = json{
                {"elevation_bonus", 10.0f}, // Slight elevation increase
                {"sediment_richness", 1.5f}, // Rich sediments
                {"features", json::array({"wetlands", "fertile_soil", "natural_harbors"})}
            }
        };

        applyRegionalInfluence(delta_influence, x, y, 5.0f); // 5-tile radius
    }
};
```

**Climate Integration**:
```cpp
void updateGlobalClimate(float cycle_duration) {
    // 1. Volcanic CO2 emissions
    float volcanic_co2 = calculateVolcanicCO2Emissions();
    atmospheric_co2 += volcanic_co2;

    // 2. Weathering CO2 absorption
    float weathering_absorption = calculateWeatheringRate() * cycle_duration;
    atmospheric_co2 -= weathering_absorption;

    // 3. Temperature from CO2 (greenhouse effect)
    global_temperature = BASE_TEMPERATURE + log(atmospheric_co2 / BASELINE_CO2) * CLIMATE_SENSITIVITY;

    // 4. Update sea level based on new temperature
    ocean_level.updateSeaLevel(global_temperature, volcanic_co2);

    // 5. Regional climate effects
    updateRegionalClimate();
}

void updateRegionalClimate() {
    for (int y = 0; y < map_height; ++y) {
        for (int x = 0; x < map_width; ++x) {
            WorldTile tile = world_map->getTile(x, y);
            float elevation = tile.getElevation();

            // Temperature varies with elevation (lapse rate: -6.5°C/km)
            float local_temperature = global_temperature - (elevation / 1000.0f) * 6.5f;

            // Distance from ocean affects temperature (continental effect)
            float ocean_distance = calculateDistanceToOcean(x, y);
            float continental_effect = ocean_distance / 1000.0f * 2.0f; // +2°C per 1000km inland
            local_temperature += continental_effect;

            tile.setTemperature(local_temperature);

            // Humidity based on proximity to water and temperature
            float humidity = calculateHumidity(ocean_distance, local_temperature, elevation);
            tile.setHumidity(humidity);
        }
    }
}
```

**Expected Results**:
- Surface elevation: -5,000m → 0m (modern sea level)
- Formation of permanent river networks draining to oceans
- Coastal features: deltas, beaches, fjords, cliffs
- Climate zones established based on elevation and ocean proximity
- Sedimentary basins filled with eroded material

### Phase 4: Carboniferous Period (35 cycles × 10M years = 0.35 Ga)

**Dynamic Forest System**:
```cpp
// Forest seeds are continuously created and destroyed by geological events
void manageDynamicForests(GMap& map, int cycle) {
    // Always seed new forests with base probability
    seedNewForests(map);

    // Existing forests attempt to expand
    expandExistingForests(map);

    // Geological events destroy forests
    destroyForestsByGeologicalEvents(map, cycle);
}

void seedNewForests(GMap& map) {
    for (each_tile) {
        float seed_probability = base_forest_seed_rate;

        // Environmental modifiers
        if (tile.getTemperature() >= 15.0f && tile.getTemperature() <= 35.0f) {
            seed_probability *= temperature_bonus;  // Optimal temperature
        }
        if (tile.getHumidity() > 0.4f) {
            seed_probability *= humidity_bonus;     // Sufficient moisture
        }
        if (tile.getElevation() < 2000.0f) {
            seed_probability *= elevation_bonus;    // Below treeline
        }

        if (random() < seed_probability) {
            createForestSeed(tile);
        }
    }
}

void expandExistingForests(GMap& map) {
    for (each_forest_tile) {
        for (each_neighbor) {
            if (isViableForForest(neighbor) && random() < forest_expansion_rate) {
                expandForestTo(neighbor);
            }
        }
    }
}

void destroyForestsByGeologicalEvents(GMap& map, int cycle) {
    // Volcanic eruptions destroy forests
    for (auto& volcanic_event : getCurrentVolcanicEvents(cycle)) {
        destroyForestsInRadius(volcanic_event.center, volcanic_event.destruction_radius);
    }

    // Major floods destroy forests
    for (each_tile_with_extreme_water_level) {
        if (tile.getWaterLevel() > flood_destruction_threshold) {
            destroyForest(tile);
        }
    }

    // Extreme erosion destroys forests
    for (each_tile) {
        if (getCurrentErosionRate(tile) > erosion_destruction_threshold) {
            destroyForest(tile);
        }
    }

    // Climate extremes destroy forests
    if (tile.getTemperature() < -10.0f || tile.getTemperature() > 50.0f) {
        destroyForest(tile);
    }
}
```

**Carbon Region System**:
```cpp
struct CarbonRegion {
    uint32_t region_id;
    float center_x, center_y;           // Position (world coordinates)
    float radius;                       // Affected area (tiles)
    float carbon_mass;                  // Total carbon content (megatons)
    int formation_cycle;                // When this region was created
    CarbonType type;                    // Current state: COAL, OIL, GAS

    // Tectonic attachment
    uint32_t attached_tectonic_region_id;    // Moves with tectonic plate
    float attachment_strength;               // 0.8-1.0 (strong attachment)

    // Conversion tracking
    float original_coal_mass;               // Initial coal amount
    float conversion_progress;              // 0.0-1.0 (coal→oil conversion)
    int cycles_underwater;                  // How long submerged
};

enum class CarbonType {
    COAL,           // Terrestrial formation
    OIL,            // Marine conversion
    NATURAL_GAS     // Late-stage oil maturation
};
```

**Forest Growth and Coal Formation**:
```cpp
void processForestGrowth(int cycle) {
    for (int y = 0; y < map_height; ++y) {
        for (int x = 0; x < map_width; ++x) {
            WorldTile tile = world_map->getTile(x, y);

            if (isForestSuitable(tile, x, y)) {
                float biomass_potential = calculateBiomassPotential(tile);
                float forest_density = std::min(1.0f, biomass_potential);

                if (forest_density > 0.3f) { // Minimum density for coal formation
                    // Create or enhance carbon region
                    CarbonRegion* existing_region = findNearestCarbonRegion(x, y, 5.0f);

                    if (existing_region) {
                        // Add to existing region
                        enhanceCarbonRegion(*existing_region, forest_density);
                    } else {
                        // Create new carbon region
                        CarbonRegion new_region = createCarbonRegion(x, y, forest_density, cycle);
                        carbon_regions.push_back(new_region);
                    }
                }
            }
        }
    }

    // Merge nearby carbon regions
    mergeCarbonRegions();
}

bool isForestSuitable(const WorldTile& tile, int x, int y) {
    float elevation = tile.getElevation();
    float temperature = tile.getTemperature();
    float humidity = tile.getHumidity();

    return (elevation > current_sea_level + 10.0f) &&     // Above sea level
           (elevation < current_sea_level + 1000.0f) &&   // Not too mountainous
           (temperature >= 10.0f && temperature <= 30.0f) && // Temperate range
           (humidity > 0.4f) &&                           // Sufficient moisture
           (!tile.getFlag(TileFlags::VOLCANIC_ACTIVE));   // Not volcanically active
}

CarbonRegion createCarbonRegion(int x, int y, float forest_density, int cycle) {
    CarbonRegion region;
    region.region_id = next_carbon_region_id++;
    region.center_x = x;
    region.center_y = y;
    region.radius = forest_density * 8.0f; // Dense forests = larger regions
    region.carbon_mass = forest_density * 100.0f; // Base carbon mass (megatons)
    region.formation_cycle = cycle;
    region.type = CarbonType::COAL;

    // Attach to nearest tectonic region
    TectonicRegion* nearest_tectonic = findNearestTectonicRegion(x, y);
    if (nearest_tectonic) {
        region.attached_tectonic_region_id = nearest_tectonic->region_id;
        region.attachment_strength = 0.9f; // Strong attachment
    }

    return region;
}
```

**Carbon Region Movement and Merging**:
```cpp
void updateCarbonRegionMovement(float time_step) {
    for (auto& carbon_region : carbon_regions) {
        TectonicRegion* tectonic = getTectonicRegion(carbon_region.attached_tectonic_region_id);

        if (tectonic) {
            // Move with tectonic plate
            carbon_region.center_x += tectonic->velocity_x * carbon_region.attachment_strength * time_step;
            carbon_region.center_y += tectonic->velocity_y * carbon_region.attachment_strength * time_step;

            // Handle tectonic collisions
            if (tectonic->isInCollision()) {
                float collision_stress = tectonic->getCollisionIntensity();

                if (collision_stress > 2.0f) {
                    // High stress can redistribute carbon deposits
                    redistributeCarbonDeposit(carbon_region, collision_stress);
                }
            }
        }
    }
}

void mergeCarbonRegions() {
    for (size_t i = 0; i < carbon_regions.size(); ++i) {
        for (size_t j = i + 1; j < carbon_regions.size(); ++j) {
            float distance = calculateDistance(carbon_regions[i].center, carbon_regions[j].center);
            float merge_threshold = (carbon_regions[i].radius + carbon_regions[j].radius) * 0.8f;

            if (distance < merge_threshold) {
                // Merge regions
                CarbonRegion merged;
                merged.region_id = next_carbon_region_id++;
                merged.center_x = (carbon_regions[i].center_x + carbon_regions[j].center_x) / 2.0f;
                merged.center_y = (carbon_regions[i].center_y + carbon_regions[j].center_y) / 2.0f;
                merged.radius = sqrt(carbon_regions[i].radius * carbon_regions[i].radius +
                                   carbon_regions[j].radius * carbon_regions[j].radius);
                merged.carbon_mass = carbon_regions[i].carbon_mass + carbon_regions[j].carbon_mass;
                merged.formation_cycle = std::min(carbon_regions[i].formation_cycle, carbon_regions[j].formation_cycle);
                merged.type = CarbonType::COAL; // Combined regions start as coal

                // Remove original regions and add merged
                carbon_regions.erase(carbon_regions.begin() + j);
                carbon_regions.erase(carbon_regions.begin() + i);
                carbon_regions.push_back(merged);

                break; // Start over after modification
            }
        }
    }
}
```

**Coal to Oil Conversion**:
```cpp
void processCoalToOilConversion(int current_cycle) {
    const float CONVERSION_RATE_PER_CYCLE = 0.05f; // 5% per cycle (10M years)
    const float NATURAL_GAS_RATIO = 0.15f;         // 15% of oil becomes gas

    for (auto& coal_region : carbon_regions) {
        if (coal_region.type == CarbonType::COAL && isRegionUnderwater(coal_region)) {
            coal_region.cycles_underwater++;

            // Only convert if underwater for multiple cycles (pressure + time)
            if (coal_region.cycles_underwater >= 3) {
                float conversion_amount = coal_region.carbon_mass * CONVERSION_RATE_PER_CYCLE;

                if (conversion_amount > 1.0f) { // Minimum threshold
                    // Create oil region at same location
                    CarbonRegion oil_region = coal_region;
                    oil_region.region_id = next_carbon_region_id++;
                    oil_region.type = CarbonType::OIL;
                    oil_region.carbon_mass = conversion_amount;
                    oil_region.radius *= 0.8f; // Oil regions are more concentrated
                    oil_region.formation_cycle = current_cycle;

                    // Create natural gas as byproduct
                    if (conversion_amount > 10.0f) {
                        CarbonRegion gas_region = oil_region;
                        gas_region.region_id = next_carbon_region_id++;
                        gas_region.type = CarbonType::NATURAL_GAS;
                        gas_region.carbon_mass = conversion_amount * NATURAL_GAS_RATIO;
                        gas_region.radius *= 1.2f; // Gas spreads more widely

                        carbon_regions.push_back(gas_region);
                    }

                    // Reduce original coal region
                    coal_region.carbon_mass -= conversion_amount;
                    coal_region.conversion_progress += CONVERSION_RATE_PER_CYCLE;

                    // Add oil region
                    carbon_regions.push_back(oil_region);

                    // Remove depleted coal regions
                    if (coal_region.carbon_mass < 1.0f) {
                        // Mark for removal
                        coal_region.carbon_mass = 0.0f;
                    }
                }
            }
        }
    }

    // Remove depleted regions
    carbon_regions.erase(
        std::remove_if(carbon_regions.begin(), carbon_regions.end(),
                      [](const CarbonRegion& region) { return region.carbon_mass < 1.0f; }),
        carbon_regions.end()
    );
}

bool isRegionUnderwater(const CarbonRegion& region) {
    // Check if majority of region is below sea level
    int underwater_tiles = 0;
    int total_tiles = 0;

    int radius_int = static_cast<int>(region.radius);
    for (int dy = -radius_int; dy <= radius_int; ++dy) {
        for (int dx = -radius_int; dx <= radius_int; ++dx) {
            float distance = sqrt(dx*dx + dy*dy);
            if (distance <= region.radius) {
                int x = static_cast<int>(region.center_x) + dx;
                int y = static_cast<int>(region.center_y) + dy;

                if (world_map->isValidCoordinate(x, y)) {
                    WorldTile tile = world_map->getTile(x, y);
                    total_tiles++;

                    if (tile.getElevation() < current_sea_level - 50.0f) { // 50m underwater minimum
                        underwater_tiles++;
                    }
                }
            }
        }
    }

    return (static_cast<float>(underwater_tiles) / total_tiles) > 0.6f; // 60% underwater
}
```

**Resource Application to Tiles**:
```cpp
void applyCarbonRegionsToTiles() {
    for (const auto& region : carbon_regions) {
        int radius_int = static_cast<int>(region.radius);

        for (int dy = -radius_int; dy <= radius_int; ++dy) {
            for (int dx = -radius_int; dx <= radius_int; ++dx) {
                float distance = sqrt(dx*dx + dy*dy);
                if (distance <= region.radius) {
                    int x = static_cast<int>(region.center_x) + dx;
                    int y = static_cast<int>(region.center_y) + dy;

                    if (world_map->isValidCoordinate(x, y)) {
                        WorldTile tile = world_map->getTile(x, y);
                        float influence = 1.0f - (distance / region.radius); // Gradient from center

                        // Apply resource based on carbon type
                        switch (region.type) {
                            case CarbonType::COAL:
                                addResourceToTile(tile, "coal", region.carbon_mass * influence);
                                break;
                            case CarbonType::OIL:
                                addResourceToTile(tile, "petroleum", region.carbon_mass * influence);
                                break;
                            case CarbonType::NATURAL_GAS:
                                addResourceToTile(tile, "natural_gas", region.carbon_mass * influence);
                                break;
                        }
                    }
                }
            }
        }
    }
}
```

**Expected Results**:
- **Dynamic forest evolution**: Forests continuously seed, expand, and get destroyed by geological events
- **Multiple forest generations**: Layers of carbon deposits from different geological periods
- **Realistic forest patterns**: Forests avoid active volcanic zones, extreme climates, and flood-prone areas
- **Coal deposit diversity**: Various ages and qualities of coal from different forest cycles
- **Oil formation**: Submerged forest areas naturally convert to petroleum over time
- **Geological storytelling**: Each coal seam represents a specific period of forest growth and burial

### Phase 5: Pre-Faunal Stabilization (15 cycles × 10M years = 0.15 Ga)

**Maturation-Only Phase - No New Formation**:
```cpp
void processStabilizationCycle(int cycle) {
    // 1. HYDROCARBON MATURATION (no new formation)
    continueHydrocarbonMaturation();

    // 2. PETROLEUM MIGRATION TO TRAPS
    migratePetroleumToGeologicalTraps();

    // 3. FINAL EROSION PHASE
    applyFinalErosion();

    // 4. SOIL LAYER DEVELOPMENT
    developSoilLayers();

    // 5. CLIMATE STABILIZATION
    stabilizeClimate();

    // 6. GEOLOGICAL STRUCTURE FINALIZATION
    finalizeGeologicalStructures();
}
```

**Continued Hydrocarbon Maturation**:
```cpp
void continueHydrocarbonMaturation() {
    const float REDUCED_CONVERSION_RATE = 0.03f; // Slower than Carboniferous (3% per cycle)

    for (auto& coal_region : carbon_regions) {
        if (coal_region.type == CarbonType::COAL && isRegionUnderwater(coal_region)) {
            // Continue coal→oil conversion at reduced rate
            float conversion_amount = coal_region.carbon_mass * REDUCED_CONVERSION_RATE;

            if (conversion_amount > 0.5f) { // Lower threshold for final conversion
                processCoalToOilConversion(coal_region, conversion_amount);
            }
        }
    }

    // Oil→Gas maturation for old oil deposits
    for (auto& oil_region : carbon_regions) {
        if (oil_region.type == CarbonType::OIL) {
            int oil_age = current_cycle - oil_region.formation_cycle;

            if (oil_age > 10 && oil_region.carbon_mass > 5.0f) { // Mature oil deposits
                float gas_conversion = oil_region.carbon_mass * 0.02f; // 2% oil→gas

                createNaturalGasFromOil(oil_region, gas_conversion);
            }
        }
    }
}
```

**Petroleum Migration to Geological Traps**:
```cpp
struct GeologicalTrap {
    float center_x, center_y;
    float capacity;                    // Maximum hydrocarbon storage
    float current_fill;               // Current hydrocarbon content
    TrapType type;                    // ANTICLINE, FAULT, SALT_DOME, STRATIGRAPHIC
    float formation_cycle;            // When trap was formed
};

enum class TrapType {
    ANTICLINE,        // Upward fold in rock layers
    FAULT_TRAP,       // Hydrocarbon blocked by fault
    SALT_DOME,        // Hydrocarbon trapped around salt intrusion
    STRATIGRAPHIC     // Trapped between different rock types
};

void migratePetroleumToGeologicalTraps() {
    // 1. Identify geological traps from tectonic history
    std::vector<GeologicalTrap> traps = identifyGeologicalTraps();

    // 2. Migrate oil and gas to nearest suitable traps
    for (auto& hydrocarbon_region : carbon_regions) {
        if (hydrocarbon_region.type == CarbonType::OIL || hydrocarbon_region.type == CarbonType::NATURAL_GAS) {
            GeologicalTrap* nearest_trap = findNearestTrap(hydrocarbon_region, traps);

            if (nearest_trap && nearest_trap->current_fill < nearest_trap->capacity) {
                float migration_amount = calculateMigrationAmount(hydrocarbon_region, *nearest_trap);

                // Move hydrocarbons to trap
                nearest_trap->current_fill += migration_amount;
                hydrocarbon_region.carbon_mass -= migration_amount;

                // Concentrate hydrocarbons in trap location
                concentrateHydrocarbonsAtTrap(*nearest_trap, hydrocarbon_region.type, migration_amount);
            }
        }
    }
}

std::vector<GeologicalTrap> identifyGeologicalTraps() {
    std::vector<GeologicalTrap> traps;

    // Find anticlines (upward folds) from tectonic compression
    for (const auto& tectonic_region : tectonic_regions) {
        if (tectonic_region.type == RegionType::CONTINENTAL) {
            // Look for elevation highs within the region (potential anticlines)
            auto high_points = findElevationHighsInRegion(tectonic_region);

            for (const auto& point : high_points) {
                GeologicalTrap anticline = {
                    .center_x = point.x,
                    .center_y = point.y,
                    .capacity = calculateAnticlineCapacity(point),
                    .current_fill = 0.0f,
                    .type = TrapType::ANTICLINE,
                    .formation_cycle = tectonic_region.formation_cycle
                };
                traps.push_back(anticline);
            }
        }
    }

    // Find fault traps from tectonic collisions
    for (const auto& collision_zone : historical_tectonic_collisions) {
        GeologicalTrap fault_trap = {
            .center_x = collision_zone.center_x,
            .center_y = collision_zone.center_y,
            .capacity = collision_zone.intensity * 50.0f, // Capacity based on collision intensity
            .current_fill = 0.0f,
            .type = TrapType::FAULT_TRAP,
            .formation_cycle = collision_zone.formation_cycle
        };
        traps.push_back(fault_trap);
    }

    return traps;
}
```

**Final Erosion and Surface Formation**:
```cpp
void applyFinalErosion() {
    const float FINAL_EROSION_RATE = 5.0f; // Reduced from active hydrological phase

    for (int y = 0; y < map_height; ++y) {
        for (int x = 0; x < map_width; ++x) {
            WorldTile tile = world_map->getTile(x, y);
            float elevation = tile.getElevation();

            if (elevation > current_sea_level) {
                // Expose coal seams through erosion
                exposeCoalSeams(tile, x, y);

                // Carve final river valleys
                carveRiverValleys(tile, x, y);

                // Form final coastal features
                if (elevation < current_sea_level + 100.0f) { // Coastal zone
                    formCoastalFeatures(tile, x, y);
                }
            }
        }
    }
}

void exposeCoalSeams(WorldTile& tile, int x, int y) {
    // Check if there are coal deposits below current elevation
    CarbonRegion* coal_region = findCoalRegionAt(x, y);

    if (coal_region && coal_region->type == CarbonType::COAL) {
        float erosion_depth = calculateErosionDepth(tile);

        if (erosion_depth > 50.0f) { // Significant erosion
            // Expose coal seam at surface
            tile.setFlag(TileFlags::SURFACE_COAL, true);

            // Add surface coal resource
            addResourceToTile(tile, "surface_coal", coal_region->carbon_mass * 0.1f);
        }
    }
}
```

**Soil Development System**:
```cpp
void developSoilLayers() {
    for (int y = 0; y < map_height; ++y) {
        for (int x = 0; x < map_width; ++x) {
            WorldTile tile = world_map->getTile(x, y);

            if (tile.getElevation() > current_sea_level + 5.0f) { // Above sea level
                SoilType soil = determineSoilType(tile, x, y);
                float soil_depth = calculateSoilDepth(tile, x, y);

                // Apply soil influence
                RegionalInfluence soil_influence = {
                    .type = "soil_development",
                    .intensity = soil_depth / 10.0f, // Depth in meters / 10
                    .properties = json{
                        {"soil_type", getSoilTypeName(soil)},
                        {"fertility", calculateSoilFertility(soil, tile)},
                        {"drainage", calculateSoilDrainage(soil, tile)},
                        {"ph_level", calculateSoilPH(soil, tile)}
                    }
                };

                applyRegionalInfluence(soil_influence, x, y, 1.0f);
            }
        }
    }
}

enum class SoilType {
    CLAY,           // Heavy, nutrient-rich, poor drainage
    SAND,           // Light, good drainage, low nutrients
    LOAM,           // Balanced, ideal for agriculture
    PEAT,           // Organic-rich, acidic, wetland areas
    ROCKY,          // Thin soil, mountainous areas
    ALLUVIAL        // River deposits, very fertile
};

SoilType determineSoilType(const WorldTile& tile, int x, int y) {
    float elevation = tile.getElevation();
    float temperature = tile.getTemperature();
    float humidity = tile.getHumidity();

    // River deltas and floodplains
    if (tile.getFlag(TileFlags::HAS_RIVER) && elevation < current_sea_level + 50.0f) {
        return SoilType::ALLUVIAL;
    }

    // Wetland areas
    if (humidity > 0.8f && elevation < current_sea_level + 20.0f) {
        return SoilType::PEAT;
    }

    // Mountainous areas
    if (elevation > current_sea_level + 1000.0f) {
        return SoilType::ROCKY;
    }

    // Climate-based soil formation
    if (temperature > 20.0f && humidity < 0.3f) {
        return SoilType::SAND; // Arid regions
    } else if (temperature < 10.0f && humidity > 0.6f) {
        return SoilType::CLAY; // Cold, wet regions
    } else {
        return SoilType::LOAM; // Temperate regions
    }
}
```

**Final Climate Stabilization**:
```cpp
void stabilizeClimate() {
    // Reduce climate variability and volcanic activity
    volcanic_activity_level *= 0.8f; // 20% reduction per cycle
    climate_stability += 0.1f;       // Increase stability

    // Establish final climate zones
    for (int y = 0; y < map_height; ++y) {
        for (int x = 0; x < map_width; ++x) {
            WorldTile tile = world_map->getTile(x, y);

            ClimateZone zone = determineClimateZone(tile, x, y);
            applyClimateZone(tile, zone);
        }
    }
}

enum class ClimateZone {
    ARCTIC,         // < -10°C, low precipitation
    SUBARCTIC,      // -10 to 0°C, moderate precipitation
    TEMPERATE,      // 0 to 20°C, variable precipitation
    SUBTROPICAL,    // 20 to 30°C, high precipitation
    TROPICAL,       // > 30°C, very high precipitation
    ARID,           // Any temperature, very low precipitation
    MEDITERRANEAN   // Warm, dry summers, mild wet winters
};
```

**Final Geological State**:
```cpp
struct FinalGeologicalState {
    // TERRAIN
    ✅ stable_continents;              // Continental masses established
    ✅ ocean_basins;                   // Deep ocean basins
    ✅ mountain_ranges;                // Various ages, realistic erosion
    ✅ river_networks;                 // Mature drainage systems
    ✅ coastal_features;               // Beaches, cliffs, deltas, fjords

    // RESOURCES
    ✅ coal_deposits;                  // Continental basins, exposed seams
    ✅ oil_fields;                     // Sedimentary basins, geological traps
    ✅ natural_gas;                    // Associated with oil, separate fields
    ✅ metal_deposits;                 // From meteorite impacts and volcanism

    // CLIMATE
    ✅ climate_zones;                  // Stable temperature/precipitation patterns
    ✅ soil_types;                     // Mature soil development
    ✅ seasonal_patterns;              // Established weather cycles

    // READY FOR INDUSTRIAL GAMEPLAY
    ✅ resource_accessibility;         // Surface coal, shallow oil, metal ores
    ✅ transportation_routes;          // Rivers, coastal access, terrain variety
    ✅ strategic_locations;            // Resource clusters, defensive positions
    ✅ environmental_challenges;       // Climate zones, terrain obstacles
};
```

**Expected Results**:
- **Complete geological maturity**: All major processes stabilized
- **Industrial-ready resources**: Coal seams exposed, oil in accessible traps
- **Realistic geography**: Mountain ranges, river valleys, coastal plains
- **Climate diversity**: Multiple biomes and environmental conditions
- **Strategic complexity**: Resource distribution creates interesting gameplay choices
- **Historical coherence**: Geological features tell the story of planetary formation

## Phase 6: Climate Simulation and Biome Generation

### Overview
After geological stabilization, run advanced climate simulation using mobile WindRegions and Inter-Tropical Convergence Zones (ITCZ) to establish realistic temperature, humidity, and wind patterns. Creates emergent weather patterns through wind region interactions rather than predefined climate zones.

**Key Innovation**: Revolutionary climate simulation using **mobile wind regions** that spawn, evolve, and interact to create emergent weather patterns. Solves the "Sahara vs Congo" problem through **Inter-Tropical Convergence Zones (ITCZ)** and **planetary rotation bands** without complex 3D atmospheric physics.

### Climate Simulation Architecture

**Phase 6.1: Landmass Analysis and ITCZ Generation**

Uses existing TectonicRegions from Phase 2:
```cpp
// Analyze continental masses from existing tectonic data
std::vector<Continent> continents = groupTectonicRegions();

// Generate water masses by inverse analysis
std::vector<WaterMass> oceans = detectOceanBasins(continents);

// Place ITCZ zones on qualifying equatorial landmasses
for (auto& continent : continents) {
    if (continent.latitude >= 0.45 && continent.latitude <= 0.55 &&
        continent.area > MIN_LANDMASS_SIZE) {
        createITCZ(continent.center, sqrt(continent.area));
    }
}
```

**ITCZ Requirements:**
- **Latitude Band**: 45-55% of map height (equatorial)
- **Minimum Landmass**: 10,000+ tiles
- **Ocean Proximity**: Within 800km for moisture source
- **Continental Heating**: Large thermal mass for convection

**Phase 6.2: WindRegion Spawning System**

Mobile WindRegions spawn from ocean masses:
```json
"wind_spawn_system": {
  "spawn_locations": "ocean_masses_only",
  "spawn_frequency": "water_mass_size / 1000",
  "initial_strength": "water_temperature * evaporation_factor",
  "spawn_distribution": "random_within_water_region_biased_toward_center"
}
```

**WindRegion Mobile Entities:**
```json
"wind_region": {
  "position": [x, y],
  "wind_strength": 1.0,        // Base intensity (decays over time)
  "wetness": 0.0,              // Moisture content (gained over ocean)
  "velocity": [vx, vy],        // Movement vector
  "wind_tokens": 100,          // Distributed to tiles
  "decay_per_move": 0.02,      // -2% strength per movement
  "decay_per_tile": 0.01       // -1% strength per tile crossed
}
```

**Phase 6.3: Movement and Planetary Circulation**

Planetary rotation bands based on real atmospheric data:
```json
"planetary_circulation": {
  "polar_jet_north": {
    "latitude_range": [0.10, 0.25],
    "direction": "west_to_east",
    "strength": 1.8
  },
  "equatorial_trades": {
    "latitude_range": [0.40, 0.60],
    "direction": "east_to_west",
    "strength": 1.5
  },
  "polar_jet_south": {
    "latitude_range": [0.75, 0.95],
    "direction": "west_to_east",
    "strength": 1.8
  }
}
```

**Movement Calculation:**
```cpp
Vector2 movement =
  planetary_rotation_band * 0.6 +     // 60% planetary circulation
  equator_to_pole_bias * 0.2 +        // 20% thermal circulation
  terrain_deflection * 0.1 +          // 10% topographic influence
  random_variation * 0.1;             // 10% chaos
```

**Phase 6.4: WindRegion Evolution and Interactions**

Dynamic evolution with ITCZ gravitational effects:
```cpp
void updateWindRegion(WindRegion& region) {
    // Gain moisture over water
    if (currentTile.isOcean()) {
        region.wetness += OCEAN_MOISTURE_GAIN;
    }

    // Repulsion from other regions = acceleration
    // NOTE: Not physical repulsion - proxy for spatial competition and turbulence
    // Prevents region stacking while creating realistic dispersion patterns
    // CRITIQUE POINT: May cause "force field" effect around ITCZ zones where regions
    // oscillate/scatter instead of converging due to attraction vs repulsion conflict.
    // Alternative approaches: density-based drift, no interaction, or collision division.
    // TODO: Implement as configurable algorithm options for empirical testing.
    float repulsion = calculateRepulsionForce(region, nearbyRegions);
    region.wind_strength += repulsion * ACCELERATION_FACTOR;

    // Movement decay
    region.wind_strength *= (1.0 - DECAY_PER_MOVE);

    // Die when too weak
    if (region.wind_strength < MINIMUM_THRESHOLD) {
        destroyRegion(region);
    }
}
```

**ITCZ Gravitational Effects:**
```cpp
void applyITCZGravity(WindRegion& region) {
    for (auto& itcz : active_itcz_zones) {
        float distance = calculateDistance(region.position, itcz.center);

        if (distance < itcz.gravitational_range) {
            // Attraction force (inverse square law)
            // NOTE: "Gravitational" metaphor for influence strength, not literal physics
            // Like saying someone has "gravitas" - clear semantic meaning for developers
            float attraction = itcz.mass / (distance * distance);
            Vector2 pull_direction = normalize(itcz.center - region.position);

            // Apply attraction
            region.velocity += pull_direction * attraction;

            // Amplification effect as region approaches
            float proximity = (itcz.range - distance) / itcz.range;
            float amplification = 1.0 + (itcz.max_amplification * proximity);

            region.wind_strength *= amplification;
            region.wetness *= amplification;
        }
    }
}
```

**Phase 6.5: Token Distribution and Climate Zone Formation**

Climate zone classification using token accumulation:
```cpp
void distributeTokens(WindRegion& region) {
    // Basic climate tokens for all regions
    int wind_tokens = static_cast<int>(region.wind_strength * 10);
    int rain_tokens = static_cast<int>(region.wetness * 10);

    WorldTile& tile = world_map.getTile(region.position);
    tile.addTokens("wind", wind_tokens);
    tile.addTokens("rain", rain_tokens);

    // Special climate zone tokens for extreme weather
    if (isHighWindZone(region)) {
        tile.addTokens("highWind", 1);  // Hostile to forests
    }
    if (isFloodZone(region)) {
        tile.addTokens("flood", 1);     // Forces wetlands/marshes
    }
    if (isHurricaneZone(region)) {
        tile.addTokens("hurricane", 1); // Specialized hurricane biome
    }
}
```

### Climate Zone Effects on Biome Generation

**Token-Based Biome Classification:**
```cpp
BiomeType classifyBiome(const WorldTile& tile) {
    int total_rain = tile.getAccumulatedTokens("rain");
    int total_wind = tile.getAccumulatedTokens("wind");
    int highWind_tokens = tile.getAccumulatedTokens("highWind");
    int flood_tokens = tile.getAccumulatedTokens("flood");
    int hurricane_tokens = tile.getAccumulatedTokens("hurricane");

    // Special climate zones override normal biome classification
    if (hurricane_tokens > 0) {
        return BiomeType::HURRICANE_ZONE;     // Specialized storm-resistant vegetation
    }
    if (flood_tokens > FLOOD_THRESHOLD) {
        return BiomeType::WETLANDS;           // Forced marshes/swamps
    }
    if (highWind_tokens > STORM_THRESHOLD) {
        // High wind prevents forest growth
        if (total_rain > 300) {
            return BiomeType::STORM_PRAIRIE;  // Grasslands that can handle wind
        } else {
            return BiomeType::BADLANDS;       // Sparse, wind-resistant vegetation
        }
    }

    // Normal biome classification using basic rain/wind tokens
    if (total_rain > 500) {
        return BiomeType::TROPICAL_RAINFOREST;
    } else if (total_rain < 50) {
        return BiomeType::HOT_DESERT;
    }
    // ... additional normal biome logic
}
```

**Climate Zone Characteristics:**
- **Hurricane Zones** → Storm-resistant palms, specialized coastal vegetation
- **Flood Zones** → Wetlands, marshes, swamp vegetation mandatory
- **High Wind Zones** → No forests allowed, prairie/badlands only
- **Normal Zones** → Standard biome classification by rain/temperature

### Geographic Climate Patterns

**Realistic Climate Formation Examples:**

**Congo Basin (Rainforest):**
```
1. Large African landmass → Strong ITCZ at equator
2. Atlantic wind regions spawn → Move east via trade winds
3. ITCZ aspiration → Convergence at Congo → Amplification ×3
4. Super-humid storms → Massive rain token distribution
5. Result: Dense rainforest biome
```

**Sahara Desert:**
```
1. Sahara latitude (25-35°N) → Outside ITCZ band
2. No convergence zone → Wind regions pass through
3. Continental distance → Low initial moisture
4. Subtropical high pressure → Air descends (simulated via movement patterns)
5. Result: Minimal rain tokens → Desert biome
```

### Configuration Integration

**Climate Configuration (Hot-Reloadable):**
```json
{
  "climate_simulation": {
    "wind_spawn_system": {
      "base_spawn_rate": 0.1,
      "ocean_size_factor": 0.001,
      "max_concurrent_regions": 200
    },
    "planetary_circulation": {
      "trade_winds_strength": 1.5,
      "jet_stream_strength": 1.8,
      "calm_zone_chaos": 0.3
    },
    "itcz_system": {
      "latitude_band": [0.45, 0.55],
      "min_landmass_size": 10000,
      "max_ocean_distance": 800,
      "amplification_max": 3.0
    },
    "storm_thresholds": {
      "high_wind_min": 2.0,
      "flood_wetness_min": 1.5,
      "hurricane_wind_min": 2.5,
      "hurricane_rain_min": 2.0
    }
  }
}
```

**Note**: All parameters are hot-reloadable via the modular configuration system. Magic numbers are intentionally externalizable for real-time tuning during development - adjust values, save config, see immediate results without recompilation.

### Performance Characteristics

**Computational Complexity:**
- **Wind Regions**: O(n) for n active regions (~50-200 simultaneously)
- **ITCZ Calculations**: O(m) for m convergence zones (~5-15 globally)
- **Token Distribution**: O(tiles_visited) per region movement
- **Total per cycle**: O(n × average_movement_distance)

**Memory Usage:**
- **WindRegion**: 32 bytes per region
- **ITCZ Zone**: 24 bytes per zone
- **Token accumulation**: Uses existing tile data structure
- **Estimated total**: <5MB for global weather simulation

**Generation Time:**
- **Landmass analysis**: 1-2 seconds (one-time setup)
- **Per simulation cycle**: 10-50ms for 100-200 wind regions
- **Full climate stabilization**: 100-500 cycles → 10-30 seconds total

## Phase 7: Budget Assignment and Natural Features

### Random Budget Assignment (Normal Distribution)

After climate and hydrology stabilization, assign budget scores to each tile using a bell curve distribution:

```cpp
void assignBudgetScores(GMap& map) {
    std::random_device rd;
    std::mt19937 gen(rd());
    std::normal_distribution<float> budget_dist(0.0f, 3.0f); // Mean=0, σ=3

    for (int y = 0; y < map.getHeight(); y++) {
        for (int x = 0; x < map.getWidth(); x++) {
            float budget_value = budget_dist(gen);
            int8_t budget = static_cast<int8_t>(std::clamp(budget_value, -10.0f, 10.0f));

            WorldTile tile = map.getTile(x, y);
            tile.setTargetBudgetScore(budget);
        }
    }
}
```

**Budget Distribution:**
- **68%** of tiles: budget score -3 to +3
- **95%** of tiles: budget score -6 to +6
- **Rare extremes**: -10/-9 and +9/+10 scores for unique locations

### Natural Features Placement (Uniform Random)

Place natural geological features randomly across the map with uniform distribution:

```cpp
void placeNaturalFeatures(GMap& map, FeatureManager& feature_manager) {
    std::random_device rd;
    std::mt19937 gen(rd());
    std::uniform_real_distribution<float> prob_dist(0.0f, 1.0f);

    const float FEATURE_CHANCE = 0.05f; // 5% of tiles get features

    for (int y = 0; y < map.getHeight(); y++) {
        for (int x = 0; x < map.getWidth(); x++) {
            if (prob_dist(gen) < FEATURE_CHANCE) {
                // Natural features from gameData configuration
                std::vector<std::string> available_features = {
                    "cave", "hot_spring", "canyon", "plateau",
                    "marsh", "oasis", "geyser", "cliff", "gorge",
                    "natural_bridge", "sinkhole", "spring"
                };

                std::uniform_int_distribution<int> feature_dist(0, available_features.size() - 1);
                std::string feature = available_features[feature_dist(gen)];

                feature_manager.placeFeature(feature, x, y);
            }
        }
    }
}
```

## Phase 8: Resource Region Conversion

### Regional Influence to Resource Features

Convert abstract regional influences from geological simulation into concrete resource features on the map:

```cpp
void convertRegionsToFeatures(GMap& map, FeatureManager& feature_manager) {
    std::random_device rd;
    std::mt19937 gen(rd());

    for (auto& region : map.getRegionalInfluences()) {
        // Semi-random feature count per region (1-8 features)
        std::uniform_int_distribution<int> feature_count_dist(1, 8);
        int feature_count = feature_count_dist(gen);

        for (int i = 0; i < feature_count; i++) {
            // Random position within region radius
            auto [x, y] = getRandomPositionInRegion(region);

            // Region influence strength determines feature quality
            float distance = getDistanceFromCenter(region, x, y);
            float influence_strength = region.getInfluenceAt(distance);

            // Place appropriate features based on region type
            placeRegionalFeature(region, x, y, influence_strength, feature_manager);
        }
    }
}

void placeRegionalFeature(const RegionalInfluence& region, int x, int y,
                         float strength, FeatureManager& feature_manager) {

    if (region.getType() == "tectonic_plate") {
        // Tectonic regions: metal deposits
        if (strength > 0.8f) {
            feature_manager.placeFeature("rich_iron_ore", x, y);
        } else if (strength > 0.6f) {
            feature_manager.placeFeature("copper_deposit", x, y);
        } else if (strength > 0.4f) {
            feature_manager.placeFeature("tin_deposit", x, y);
        } else {
            feature_manager.placeFeature("stone_quarry", x, y);
        }
    }

    else if (region.getType() == "carbon_region") {
        // Carbon regions: coal and oil deposits
        if (strength > 0.7f) {
            feature_manager.placeFeature("rich_coal_seam", x, y);
        } else if (strength > 0.5f) {
            feature_manager.placeFeature("oil_well", x, y);
        } else if (strength > 0.3f) {
            feature_manager.placeFeature("coal_outcrop", x, y);
        } else {
            feature_manager.placeFeature("peat_bog", x, y);
        }
    }

    else if (region.getType() == "volcanic_zone") {
        // Volcanic regions: geothermal and rare minerals
        if (strength > 0.8f) {
            feature_manager.placeFeature("geothermal_vent", x, y);
        } else if (strength > 0.6f) {
            feature_manager.placeFeature("sulfur_deposit", x, y);
        } else if (strength > 0.4f) {
            feature_manager.placeFeature("obsidian_field", x, y);
        } else {
            feature_manager.placeFeature("hot_spring", x, y);
        }
    }

    else if (region.getType() == "recent_meteorite_crater") {
        // Recent meteorite impacts: rare metals and exotic materials
        // NOTE: Different from Phase 1 planetary accretion meteorites
        if (strength > 0.9f) {
            feature_manager.placeFeature("platinum_deposit", x, y);
        } else if (strength > 0.7f) {
            feature_manager.placeFeature("rare_earth_deposit", x, y);
        } else if (strength > 0.5f) {
            feature_manager.placeFeature("gold_vein", x, y);
        } else {
            feature_manager.placeFeature("impact_glass", x, y);
        }
    }
}

std::pair<int, int> getRandomPositionInRegion(const RegionalInfluence& region) {
    std::random_device rd;
    std::mt19937 gen(rd());

    // Random angle and distance within region radius
    std::uniform_real_distribution<float> angle_dist(0.0f, 2.0f * M_PI);
    std::uniform_real_distribution<float> radius_dist(0.0f, region.getRadius());

    float angle = angle_dist(gen);
    float distance = radius_dist(gen);

    int x = static_cast<int>(region.getCenterX() + distance * cos(angle));
    int y = static_cast<int>(region.getCenterY() + distance * sin(angle));

    return {x, y};
}
```

### Resource Feature Characteristics

**Influence Strength Distribution:**
- **0.8-1.0**: Premium resources (rich deposits, rare materials)
- **0.6-0.8**: High-quality resources (standard industrial deposits)
- **0.4-0.6**: Medium-quality resources (adequate for basic industry)
- **0.2-0.4**: Low-quality resources (marginal extraction)
- **0.0-0.2**: Minimal resources (trace amounts only)

**Regional Resource Mapping:**
- **Tectonic Regions**: Iron, copper, tin, stone → Industrial base materials
- **Carbon Regions**: Coal, oil, peat → Energy resources
- **Volcanic Zones**: Geothermal, sulfur, obsidian → Specialized materials
- **Recent Meteorite Craters**: Platinum, rare earths, gold → Advanced technology materials

**Feature Distribution:**
- **1-8 features per region** (semi-random count)
- **Position**: Random within region radius
- **Quality**: Determined by distance from region center (closer = better)
- **Type**: Fixed by region type, quality by influence strength

## Planetary Mass Conservation System

### Core-Based Mass Conservation

To maintain realistic mass conservation throughout geological simulation, the system uses a planetary core that absorbs eroded material and releases it through volcanic activity:

```cpp
struct PlanetaryCore {
    float core_mass;                    // Accumulated eroded material
    float surface_mass;                 // Current surface terrain mass
    float max_core_capacity;            // Core saturation threshold
    float volcanic_overflow_rate;       // Rate of volcanic material expulsion (0.1-0.3)
    float total_planetary_mass;         // Constant after Phase 1 (meteorite bombardment)

    // Derived values
    float core_pressure_ratio;          // core_mass / max_core_capacity
    int pending_volcanic_events;        // Queued volcanic eruptions
};

void updateMassConservation(PlanetaryCore& core) {
    // Validation: Total mass conservation
    assert(core.core_mass + core.surface_mass == core.total_planetary_mass);

    // Core overflow triggers volcanic activity
    if (core.core_mass > core.max_core_capacity) {
        float overflow = core.core_mass - core.max_core_capacity;
        float volcanic_expulsion = overflow * core.volcanic_overflow_rate;

        // Transfer mass from core to pending volcanic events
        core.core_mass -= volcanic_expulsion;

        // Queue volcanic events proportional to overflow
        int new_volcanic_events = static_cast<int>(volcanic_expulsion / volcanic_event_mass_threshold);
        core.pending_volcanic_events += new_volcanic_events;

        // Store remaining material for future volcanism
        core.pending_volcanic_material += volcanic_expulsion;
    }
}
```

### Erosion to Core Transfer

All erosion processes transfer material to the planetary core rather than redistributing on surface:

```cpp
void erodeToCore(WorldTile& tile, float erosion_amount, PlanetaryCore& core) {
    // Remove material from surface
    float current_elevation = tile.getElevation();
    tile.setElevation(current_elevation - erosion_amount);

    // Transfer to planetary core
    core.core_mass += erosion_amount;
    core.surface_mass -= erosion_amount;

    // Track erosion for geological history
    tile.addFlag(TileFlags::ERODED_THIS_CYCLE);
}

// Apply to all erosion processes
void riverErosion(WorldTile& tile, float water_flow, PlanetaryCore& core) {
    if (water_flow > erosion_threshold) {
        float erosion_amount = water_flow * erosion_rate_factor;
        erodeToCore(tile, erosion_amount, core);
    }
}

void glacialErosion(WorldTile& tile, float ice_thickness, PlanetaryCore& core) {
    float erosion_amount = ice_thickness * glacial_erosion_factor;
    erodeToCore(tile, erosion_amount, core);
}
```

### Volcanic Overflow System

Core overflow creates realistic volcanic activity that returns material to surface:

```cpp
void processVolcanicOverflow(GMap& map, PlanetaryCore& core) {
    while (core.pending_volcanic_events > 0) {
        // Select volcanic location based on geological factors
        auto [x, y] = selectVolcanicLocation(map);

        // Calculate eruption magnitude
        float eruption_material = core.pending_volcanic_material / core.pending_volcanic_events;

        // Create volcanic eruption
        createVolcanicEruption(map, x, y, eruption_material);

        // Update core state
        core.surface_mass += eruption_material;
        core.pending_volcanic_material -= eruption_material;
        core.pending_volcanic_events--;
    }
}

std::pair<int, int> selectVolcanicLocation(const GMap& map) {
    // Prefer locations with:
    // - Active tectonic boundaries
    // - Existing volcanic history
    // - High core pressure influence

    std::vector<std::pair<int, int>> candidate_locations;

    for (auto& tectonic_region : map.getTectonicRegions()) {
        if (tectonic_region.getActivity() > volcanic_activity_threshold) {
            // Add boundary tiles as candidates
            auto boundary_tiles = tectonic_region.getBoundaryTiles();
            candidate_locations.insert(candidate_locations.end(),
                                     boundary_tiles.begin(), boundary_tiles.end());
        }
    }

    // Weight by distance from core pressure points
    return selectWeightedRandom(candidate_locations);
}

void createVolcanicEruption(GMap& map, int center_x, int center_y, float material_amount) {
    // Deposit material in volcanic pattern
    int eruption_radius = static_cast<int>(sqrt(material_amount / volcanic_density_factor));

    for (int dy = -eruption_radius; dy <= eruption_radius; dy++) {
        for (int dx = -eruption_radius; dx <= eruption_radius; dx++) {
            float distance = sqrt(dx*dx + dy*dy);
            if (distance <= eruption_radius) {
                int x = center_x + dx;
                int y = center_y + dy;

                if (map.isValidCoordinate(x, y)) {
                    // Volcanic deposition decreases with distance
                    float deposition_factor = 1.0f - (distance / eruption_radius);
                    float local_deposition = material_amount * deposition_factor / (M_PI * eruption_radius * eruption_radius);

                    WorldTile tile = map.getTile(x, y);
                    tile.setElevation(tile.getElevation() + local_deposition);

                    // Mark as volcanic terrain
                    tile.setFlag(TileFlags::VOLCANIC_DEPOSIT, true);
                }
            }
        }
    }
}
```

### Mass Conservation Benefits

**Perfect Conservation:**
- Total planetary mass remains constant after Phase 1
- All eroded material is accounted for in core
- Volcanic activity returns material to surface
- No material created or destroyed

**Realistic Geological Activity:**
- Core pressure drives continuing volcanism
- Volcanic activity decreases as core pressure reduces
- Natural equilibrium between erosion and volcanic deposition
- Geological activity persists throughout simulation

**Simplified Implementation:**
- No complex sediment transport calculations
- No surface redistribution algorithms
- Single mass conservation equation
- Volcanic activity emerges naturally from core overflow

**Gameplay Benefits:**
- Ongoing geological activity creates dynamic world
- Volcanic regions provide unique resource opportunities
- Erosion/volcanism balance creates varied topography
- Long-term geological processes affect industrial planning

## Technical Architecture

### WorldTileData Structure (32 bytes)
```cpp
struct WorldTileData {
    // Terrain (11 bytes) - Always accessed, geological simulation ready
    uint16_t terrain_type_id;       // 65k terrain types
    uint16_t biome_type_id;         // Biome classification
    uint16_t elevation;             // -11km to +32km range
    int16_t temperature;            // -3276°C to +3276°C (0.1°C precision)
    uint8_t humidity;               // 0-100% (0.4% precision)
    uint8_t wind_data;              // 4 bits direction + 4 bits intensity
    uint8_t water_level;            // Accumulated water for river formation

    // Metadata (21 bytes) - Generation and gameplay
    int8_t target_budget_score;     // -10 to +10
    uint32_t regional_influence_id; // → Regional influence data
    uint8_t influence_strength;     // 0-255
    uint32_t tile_flags;            // State flags
    uint32_t feature_set_id;        // → Feature collection
    uint8_t padding2[7];            // Future expansion space
};
```

### RegionalInfluence System
- **TectonicRegions**: Mobile circular regions with physics
- **CarbonRegions**: Carbon deposit tracking with tectonic attachment
- **VolcanicZones**: Created dynamically at tectonic collisions
- **SeaLevelInfluence**: Global parameter affecting all processes

### FeatureManager Integration
- Simple helper interface for placing geological features
- Handles feature set creation and tile updates automatically
- Used by generation algorithms, doesn't do generation itself

## Performance Characteristics

### Memory Usage
- **1M tiles**: 32MB core tile data (increased from 24MB for hydrology)
- **Feature sets**: Sparse, shared between similar tiles
- **Geological regions**: ~10-50 regions vs millions of tiles
- **Climate wind regions**: ~50-200 mobile regions during simulation
- **Total estimated**: <125MB for complete planetary geology + advanced climate simulation

### Computational Complexity
- **Tectonic simulation**: O(n²) for n regions (~25 regions = 625 operations)
- **Meteorite impacts**: O(k) for k impacts per wave
- **Sea level effects**: O(tiles) single pass
- **Carbon processes**: O(regions) sparse operations
- **Climate simulation**: O(n × movement_distance) for n wind regions (~50-200 regions)

### Generation Time
- **Geological phases**: 10-20 seconds for complete 4.65 billion year simulation
- **Climate simulation**: 10-30 seconds for 100-500 climate cycles
- **Total estimated**: 20-50 seconds for complete world generation
- **Progressive**: Can display intermediate results during generation
- **Deterministic**: Same seed produces identical geology and climate

## Implementation Priority

### Phase 1: Core Architecture (1-2 weeks)
1. RegionalInfluence and RegionalInfluenceManager classes
2. TectonicRegion basic structure and physics
3. Simple meteorite impact system
4. Basic tile elevation/temperature updates

### Phase 2: Full Geology (2-3 weeks)
1. Complete tectonic interaction system
2. Dynamic sea level integration
3. Carbon region formation and conversion
4. Volcanic zone creation

### Phase 3: Polish (1-2 weeks)
1. Parameter tuning for realistic results
2. Performance optimization
3. Generation progress reporting
4. Validation and debugging tools

## Scientific Accuracy vs Gameplay

### Scientifically Inspired Elements
- ✅ Late Heavy Bombardment period
- ✅ Planetary differentiation (metals sink to core)
- ✅ Tectonic processes creating mountains/valleys
- ✅ Carboniferous period forest→coal formation
- ✅ Marine conditions for oil formation
- ✅ Sea level variations affecting geology
- ✅ ITCZ formation from continental heating
- ✅ Planetary circulation bands (trade winds, jet streams)
- ✅ Realistic climate differentiation (Congo vs Sahara)

### Gameplay Simplifications
- ⚠️ Tectonic regions as circles (not realistic plate shapes)
- ⚠️ Fixed map size instead of planetary expansion
- ⚠️ Simplified core-mantle dynamics
- ⚠️ Compressed timescales for practical generation
- ⚠️ 2D wind simulation instead of 3D atmospheric layers
- ⚠️ Discrete token system instead of continuous climate fields
- ⚠️ Mobile wind regions as simplified weather systems

### Result
**Plausible geology and climate** that creates **interesting gameplay** with **emergent weather patterns** without requiring PhD in atmospheric science to understand or debug.

## Integration Points

### WorldGenerationModule
- Orchestrates all geological phases
- Manages simulation state and progression
- Provides query interface for generated geology

### FeatureManager
- Places geological features based on simulation results
- Handles feature set optimization and tile updates
- Simple interface for placement algorithms

### RegionalInfluenceManager
- Core system managing all regional effects
- Handles tectonic regions, carbon regions, volcanic zones
- Provides influence application to tiles

### Future Extensions
- **Dynamic climate**: Real-time weather systems during gameplay
- **Mineral resource modeling**: Detailed ore deposit formation
- **Advanced erosion**: River network evolution during gameplay
- **Geological time compression**: Speed up/slow down specific phases
- **Seasonal climate**: Monthly/yearly climate variations
- **Climate change**: Long-term climate evolution from industrial activity

---

**Status**: System designed and ready for implementation. All major components specified with clear interfaces and realistic performance targets. Scientific accuracy balanced with implementation complexity for practical game development.