diff --git a/CLAUDE.md b/CLAUDE.md
index 3d91d0e..084680e 100644
--- a/CLAUDE.md
+++ b/CLAUDE.md
@@ -4,7 +4,7 @@ This file provides guidance to Claude Code (claude.ai/code) when working with co
 
 ## Project Overview
 
-LiDAR archaeological processing pipeline that generates 18 terrain visualizations from LAZ/LAS point clouds. Runs in Docker with optional NVIDIA GPU acceleration (CuPy). Designed for French LiDAR HD data in Lambert 93 (EPSG:2154).
+LiDAR archaeological processing pipeline that generates 17 terrain visualizations from LAZ/LAS point clouds. Runs in Docker with optional NVIDIA GPU acceleration (CuPy). Designed for French LiDAR HD data in Lambert 93 (EPSG:2154).
 
 ## Commands
 
@@ -17,6 +17,7 @@ All commands run inside Docker. Use `./run.sh` as the primary interface.
 ./run.sh --test                                 # Run unit tests
 ./run.sh -g --file LHD_FXX_1000_6882_PTS_LAMB93_IGN69.copc  # Single file
 ./run.sh --ground-classification pmf            # Force PMF ground classification
+./run.sh -g --keep-tif                          # Keep intermediate TIFF files
 ./run.sh                                        # Print help (no args)
 ```
 
@@ -33,7 +34,7 @@ docker run --rm --gpus all -v $(pwd)/input:/data/input:ro -v $(pwd)/output:/data
 - **`cli.py`** — argparse + logging setup. Entry point via `python -m lidar_pipeline`.
 - **`pipeline.py`** — `LidarArchaeoPipeline` orchestrator. `VIZ_STEPS` registry maps names to generate functions. `FilePrefixFilter` for parallel logging.
 - **`dtm.py`** — PDAL ground classification (SMRF/PMF/CSF + auto-detection) and DTM generation via scipy `binned_statistic_2d`.
-- **`visualizations.py`** — 16 `generate_*` functions + 2 IGN overlay lambdas. All take `(dem_file, basename, vis_dir, resolution)` and return a TIF path or None.
+- **`visualizations.py`** — 15 `generate_*` functions + 2 IGN overlay lambdas. All take `(dem_file, basename, vis_dir, resolution)` and return a TIF path or None.
 - **`gpu.py`** — CuPy/numpy abstraction: `HAS_GPU`, `to_gpu()`, `to_cpu()`, `xp_gaussian_filter()`, `xp_uniform_filter()`, `xp_minimum_filter()`, `gpu_cleanup()`. Falls back to CPU gracefully.
 - **`ign.py`** — IGN WMTS tile download + overlay generation for orthophoto and topographic maps.
 - **`rendering.py`** — `COLORMAPS` dict maps filename keywords to (cmap, title, legend, description). `tif_to_png()` converts TIF→WebP with legend/scale/north arrow. `generate_pdf_report()` creates A3 PDF.
@@ -56,7 +57,7 @@ Override with `--ground-classification {auto,smrf,pmf,csf}`.
 
 ### NaN handling
 
-DTM zones without LiDAR data are kept as NaN (no interpolation). Visualization functions use `_fill_nans()` and `_filter_nanaware()` to avoid NaN propagation through filters.
+DTM small gaps (< 1m from existing data) are filled using `rasterio.fill.fillnodata`. Large gaps remain as NaN. Visualization functions use `_fill_nans()` and `_filter_nanaware()` to avoid NaN propagation through filters.
 
 ### Parallel processing
 
@@ -67,5 +68,6 @@ Uses `ProcessPoolExecutor` with `'spawn'` start method (required for CUDA). Each
 - **Language**: UI messages and comments in French. Code identifiers in English.
 - **Logging**: Use `logger = logging.getLogger("lidar")`. Prefix per-file logs via `_file_filter.basename`.
 - **GPU pattern**: `arr_gpu = to_gpu(arr)` → compute → `result = to_cpu(arr_gpu)` → `gpu_cleanup()` between visualizations.
-- **Output format**: Visualizations saved as WebP (not PNG). PDF reports use `PILImage.open().convert('RGB')`.
+- **Output format**: Visualizations saved as WebP (not PNG). TIFF intermediates deleted unless `--keep-tif`. PDF reports use `PILImage.open().convert('RGB')`.
+- **Flow accumulation**: Uses numba JIT for D8 accumulation loop. Falls back to pure Python if numba unavailable.
 - **Tests**: Run only inside Docker via `./run.sh --test`. Synthetic DEM fixture in `tests/conftest.py`.
\ No newline at end of file
diff --git a/Dockerfile b/Dockerfile
index 7270bbf..91af8f8 100644
--- a/Dockerfile
+++ b/Dockerfile
@@ -31,7 +31,8 @@ RUN pip3 install --no-cache-dir \
     scipy \
     tqdm \
     Pillow \
-    pytest
+    pytest \
+    numba
 
 # Install CuPy for GPU acceleration (optional - will fallback to numpy if not available)
 RUN pip3 install --no-cache-dir cupy-cuda12x || echo "CuPy not available - GPU acceleration disabled"
diff --git a/lidar_pipeline/cli.py b/lidar_pipeline/cli.py
index 7042ece..92b5aaa 100644
--- a/lidar_pipeline/cli.py
+++ b/lidar_pipeline/cli.py
@@ -6,6 +6,7 @@ Handles argument parsing, logging configuration, and entry point.
 import argparse
 import logging
 import os
+import shutil
 import signal
 import sys
 
@@ -110,6 +111,11 @@ Exemples:
         action="store_true",
         help="Reclassifier le sol même si le fichier .las existe déjà"
     )
+    parser.add_argument(
+        "--keep-tif",
+        action="store_true",
+        help="Conserver les fichiers TIFF intermédiaires (sinon supprimés après conversion WebP)"
+    )
     parser.add_argument(
         "--ground-classification",
         choices=["auto", "smrf", "pmf", "csf"],
@@ -163,7 +169,8 @@ Exemples:
             workers=args.workers,
             force=args.force,
             ground_method=args.ground_classification,
-            force_classify=args.force_classification
+            force_classify=args.force_classification,
+            keep_tif=args.keep_tif
         )
 
         # If --file is specified, process only matching files
@@ -199,6 +206,18 @@ Exemples:
                 logger.info(f"  → {laz_file.name}")
             for laz_file in unique_files:
                 pipeline.process_file(laz_file)
+
+            # Clean up temporary files
+            logger.info("Nettoyage des fichiers temporaires...")
+            try:
+                if pipeline.temp_dir.exists():
+                    shutil.rmtree(pipeline.temp_dir)
+                temp_base = pipeline.output_dir / "temp"
+                if temp_base.exists():
+                    shutil.rmtree(temp_base)
+                logger.info("  ✓ Fichiers temporaires supprimés")
+            except Exception as e:
+                logger.warning(f"  Note: Impossible de supprimer les fichiers temporaires: {e}")
         else:
             pipeline.process_all()
     except Exception as e:
diff --git a/lidar_pipeline/dtm.py b/lidar_pipeline/dtm.py
index 433fc98..d76bf27 100644
--- a/lidar_pipeline/dtm.py
+++ b/lidar_pipeline/dtm.py
@@ -289,13 +289,25 @@ def create_dtm_fast(las_file, basename, dtm_dir, resolution):
         dtm = stat.statistic.T
         dtm = dtm[::-1, :]  # Flip Y so north is at top
 
-        # Keep NaN for zones without LiDAR data — no interpolation
+        # Fill small gaps (< 1m from existing data) while keeping large gaps as NaN
         nan_count = np.count_nonzero(np.isnan(dtm))
         if nan_count > 0:
             total = dtm.size
             nan_pct = 100.0 * nan_count / total
             logger.info(f"    {nan_count:,} pixels sans données ({nan_pct:.1f}%)")
 
+            max_gap_pixels = max(1, int(1.0 / resolution))
+            from rasterio.fill import fillnodata
+            valid_mask = ~np.isnan(dtm)
+            dtm_filled = fillnodata(dtm, mask=valid_mask, max_search_distance=max_gap_pixels)
+            small_gap_mask = np.isnan(dtm) & ~np.isnan(dtm_filled)
+            filled_count = np.count_nonzero(small_gap_mask)
+            if filled_count > 0:
+                dtm = np.where(small_gap_mask, dtm_filled, dtm)
+                logger.info(f"    {filled_count:,} petits trous comblés (< {max_gap_pixels}px)")
+            remaining = np.count_nonzero(np.isnan(dtm))
+            logger.info(f"    {remaining:,} pixels restent sans données (grands écarts)")
+
         # Save as GeoTIFF
         output_tif = dtm_dir / f"{basename}_dtm.tif"
         transform = from_bounds(min_x, min_y, max_x, max_y, width, height)
diff --git a/lidar_pipeline/pipeline.py b/lidar_pipeline/pipeline.py
index bbd498f..6a35387 100644
--- a/lidar_pipeline/pipeline.py
+++ b/lidar_pipeline/pipeline.py
@@ -3,7 +3,7 @@
 LidarArchaeoPipeline coordinates the full processing chain:
 1. Ground classification (PDAL/SMRF)
 2. DTM generation
-3. Visualization generation (18 products)
+3. Visualization generation (17 products)
 4. Rendering (WebP + PDF report)
 """
 
@@ -57,7 +57,7 @@ from .dtm import classify_ground, create_dtm_fast
 from .visualizations import (
     generate_hillshade, generate_slope, generate_aspect, generate_curvature,
     generate_lrm, generate_svf, generate_openness,
-    generate_mslrm, generate_tpi, generate_depressions, generate_sailore,
+    generate_mslrm, generate_tpi, generate_sailore,
     generate_roughness, generate_anomalies, generate_wavelet,
     generate_flow,
 )
@@ -80,7 +80,6 @@ VIZ_STEPS = [
     ('neg_open',   lambda d, b, v, r: generate_openness(d, b, v, r, positive=False)),
     ('mslrm',      generate_mslrm),
     ('tpi',        generate_tpi),
-    ('depressions', generate_depressions),
     ('sailore',    generate_sailore),
     ('roughness',  generate_roughness),
     ('anomalies',  generate_anomalies),
@@ -106,7 +105,7 @@ VIZ_STEPS = [
 class LidarArchaeoPipeline:
     """Orchestrates the LiDAR archaeological analysis pipeline."""
 
-    def __init__(self, input_dir, output_dir, resolution=0.5, workers=1, force=False, ground_method='auto', force_classify=False):
+    def __init__(self, input_dir, output_dir, resolution=0.5, workers=1, force=False, ground_method='auto', force_classify=False, keep_tif=False):
         self.input_dir = Path(input_dir)
         self.output_dir = Path(output_dir)
         self.resolution = resolution
@@ -114,6 +113,7 @@ class LidarArchaeoPipeline:
         self.force = force
         self.ground_method = ground_method
         self.force_classify = force_classify
+        self.keep_tif = keep_tif
         self.temp_dir = self.output_dir / "temp"
 
         if not self.input_dir.exists():
@@ -137,6 +137,7 @@ class LidarArchaeoPipeline:
         logger.info(f"  Force       : {'OUI' if self.force else 'non (skip existing)'}")
         logger.info(f"  Classification sol : {self.ground_method}")
         logger.info(f"  Force classif.: {'OUI' if self.force_classify else 'non'}")
+        logger.info(f"  Keep TIFF    : {'OUI' if self.keep_tif else 'non'}")
 
     def find_laz_files(self):
         """Find all LAZ/LAS files in input directory."""
@@ -171,8 +172,24 @@ class LidarArchaeoPipeline:
 
         vis_results = {}
         total = len(VIZ_STEPS)
+        elapsed_times = []
 
         for idx, (name, func) in enumerate(VIZ_STEPS, 1):
+            # When --force, delete existing TIF to ensure clean regeneration
+            if self.force:
+                for tif in file_vis_dir.glob(f"{basename}_{name}.tif"):
+                    tif.unlink(missing_ok=True)
+                # Special cases for differently-named TIFs
+                if name == 'pos_open':
+                    for tif in file_vis_dir.glob(f"{basename}_positive_openness.tif"):
+                        tif.unlink(missing_ok=True)
+                elif name == 'neg_open':
+                    for tif in file_vis_dir.glob(f"{basename}_negative_openness.tif"):
+                        tif.unlink(missing_ok=True)
+                elif name == 'hillshade':
+                    for tif in file_vis_dir.glob(f"{basename}_hillshade_multi.tif"):
+                        tif.unlink(missing_ok=True)
+
             # Check if output WebP already exists (skip unless --force)
             if not self.force:
                 # Determine expected WebP filename from the viz name
@@ -199,8 +216,14 @@ class LidarArchaeoPipeline:
                 result = func(dtm_file, basename, file_vis_dir, self.resolution)
                 vis_results[name] = result
                 elapsed = time.time() - t0
+                elapsed_times.append(elapsed)
                 if result:
-                    logger.info(f"    [{idx}/{total}] ✓ {name} ({elapsed:.1f}s)")
+                    eta = ""
+                    if len(elapsed_times) > 1:
+                        avg_time = sum(elapsed_times) / len(elapsed_times)
+                        remaining = (total - idx) * avg_time
+                        eta = f" — ETA: {remaining:.0f}s"
+                    logger.info(f"    [{idx}/{total}] ✓ {name} ({elapsed:.1f}s){eta}")
                 else:
                     logger.warning(f"    [{idx}/{total}] ✗ {name} — no output ({elapsed:.1f}s)")
             except Exception as e:
@@ -214,7 +237,7 @@ class LidarArchaeoPipeline:
         logger.info("  Conversion images WebP:")
         for name, tif_file in vis_results.items():
             if tif_file and isinstance(tif_file, Path) and tif_file.suffix == '.tif' and tif_file.exists():
-                webp_file = tif_to_png(tif_file, file_vis_dir, self.resolution)
+                webp_file = tif_to_png(tif_file, file_vis_dir, self.resolution, keep_tif=self.keep_tif)
                 if webp_file:
                     logger.info(f"    ✓ {webp_file.name}")
 
@@ -293,7 +316,7 @@ class LidarArchaeoPipeline:
             logger.info(f"Fichiers: {len(files)}")
             with ProcessPoolExecutor(max_workers=self.workers) as executor:
                 future_to_file = {
-                    executor.submit(_process_file_standalone, str(laz_file), str(self.input_dir), str(self.output_dir), self.resolution, self.force, self.ground_method, self.force_classify): laz_file
+                    executor.submit(_process_file_standalone, str(laz_file), str(self.input_dir), str(self.output_dir), self.resolution, self.force, self.ground_method, self.force_classify, self.keep_tif): laz_file
                     for laz_file in files
                 }
                 done = 0
@@ -346,16 +369,16 @@ class LidarArchaeoPipeline:
         try:
             if self.temp_dir.exists():
                 shutil.rmtree(self.temp_dir)
-            # Also clean up per-file temp directories from parallel workers
-            for d in self.output_dir.glob("temp_*"):
-                if d.is_dir():
-                    shutil.rmtree(d, ignore_errors=True)
+            # Also clean up any subdirectories inside temp/
+            temp_base = self.output_dir / "temp"
+            if temp_base.exists():
+                shutil.rmtree(temp_base)
             logger.info("  ✓ Fichiers temporaires supprimés")
         except Exception as e:
             logger.warning(f"  Note: Impossible de supprimer les fichiers temporaires: {e}")
 
 
-def _process_file_standalone(laz_file_str, input_dir, output_dir, resolution, force=False, ground_method='auto', force_classify=False):
+def _process_file_standalone(laz_file_str, input_dir, output_dir, resolution, force=False, ground_method='auto', force_classify=False, keep_tif=False):
     """Standalone function for multiprocessing — creates its own pipeline instance.
 
     Each worker gets its own temp directory to avoid file conflicts.
@@ -371,9 +394,9 @@ def _process_file_standalone(laz_file_str, input_dir, output_dir, resolution, fo
         worker_logger.addHandler(handler)
     worker_logger.addFilter(_file_filter)
 
-    pipeline = LidarArchaeoPipeline(input_dir, output_dir, resolution=resolution, workers=1, force=force, ground_method=ground_method, force_classify=force_classify)
+    pipeline = LidarArchaeoPipeline(input_dir, output_dir, resolution=resolution, workers=1, force=force, ground_method=ground_method, force_classify=force_classify, keep_tif=keep_tif)
     basename = _file_basename(laz_file_str)
-    pipeline.temp_dir = pipeline.output_dir / f"temp_{basename}"
+    pipeline.temp_dir = pipeline.output_dir / "temp" / basename
     pipeline.temp_dir.mkdir(exist_ok=True)
     laz_file = Path(laz_file_str)
     result = pipeline.process_file(laz_file)
diff --git a/lidar_pipeline/rendering.py b/lidar_pipeline/rendering.py
index 32a36ef..41aae74 100644
--- a/lidar_pipeline/rendering.py
+++ b/lidar_pipeline/rendering.py
@@ -24,7 +24,6 @@ import matplotlib
 matplotlib.use('Agg')
 import matplotlib.pyplot as plt
 from matplotlib import rcParams
-from matplotlib.gridspec import GridSpec
 from matplotlib.patches import Polygon as MplPolygon, Rectangle as RectPatch
 from mpl_toolkits.axes_grid1.inset_locator import inset_axes
 
@@ -119,14 +118,6 @@ COLORMAPS = {
         'description': 'Identifie la position topographique — utile pour repérer crêtes vs vallées à grande échelle',
         'vmin_mode': 'symmetric', 'sym_pct': (2, 98),
     },
-    'depressions': {
-        'cmap': 'YlOrRd',
-        'title': 'Dépressions (Remplissage hydrologique)',
-        'legend': 'Profondeur des cuvettes détectées (m)\nTransparent = Pas de dépression\nJaune = Dépression légère | Rouge = Dépression profonde\nMax: {vmax:.2f}m',
-        'description': 'Simule le remplissage d\'eau — détecte dolines, sinkholes, cuvettes et zones inondables',
-        'vmin_mode': 'fixed', 'vmin_val': 0,
-        'vmax_mode': 'percentile', 'vmax_pct': 99,
-    },
     'sailore': {
         'cmap': 'seismic',
         'title': 'SAILORE - LRM Auto-Adaptatif',
@@ -194,8 +185,9 @@ def _apply_colormap(data, tif_file):
             info = RGB_LEGENDS[key]
             return data, None, info['title'], info['legend'], info['description'], True
 
-    # Find matching colormap
-    for key, info in COLORMAPS.items():
+    # Find matching colormap — sort by key length descending so 'mslrm' matches before 'lrm'
+    for key in sorted(COLORMAPS.keys(), key=len, reverse=True):
+        info = COLORMAPS[key]
         if key in name:
             valid_data = np.asarray(data.compressed() if hasattr(data, 'compressed') else data.flatten())
             valid_data = valid_data[~np.isnan(valid_data)]
@@ -246,13 +238,14 @@ def _apply_colormap(data, tif_file):
     return data, 'terrain', title, 'Altitude normalisée', '', False
 
 
-def tif_to_png(tif_file, vis_dir, resolution):
+def tif_to_png(tif_file, vis_dir, resolution, keep_tif=False):
     """Convert GeoTIFF to visualization WebP with GPS coordinates, legend, and scale bar.
 
     Args:
         tif_file: Path to input GeoTIFF.
         vis_dir: Output directory for the WebP file.
         resolution: Grid resolution in m/px.
+        keep_tif: If True, keep the source TIFF after conversion.
 
     Returns:
         Path to output WebP file, or None on failure.
@@ -359,18 +352,21 @@ def tif_to_png(tif_file, vis_dir, resolution):
                 saved_vmin = None
                 saved_vmax = None
 
-            # Create figure — adapt width to data resolution for sharp rendering
-            # At high res (5000+px wide), we need a larger figure to avoid downsampling artifacts
+            # Create figure with FIXED layout for consistent data area position
+            # All visualizations use the same axes positions so they can be overlaid
             fig_width = max(20, width / 150)
-            fig_width = min(fig_width, 40)  # cap at 40 inches
-            map_aspect = height / width
-            fig = plt.figure(figsize=(fig_width, fig_width * map_aspect * 0.7 + 2.5),
-                            facecolor='white')
-            gs = GridSpec(2, 1, height_ratios=[1.0, 0.06],
-                         hspace=0.04, figure=fig,
-                         left=0.06, right=0.88, top=0.93, bottom=0.08)
+            fig_width = min(fig_width, 40)
+            fig_height = fig_width * 0.7 + 2.0  # Fixed header + footer space
+            fig = plt.figure(figsize=(fig_width, fig_height), facecolor='white')
 
-            ax = fig.add_subplot(gs[0])
+            # Fixed data area position — identical for ALL visualization types
+            # This ensures overlay/superposition works across all WebP images
+            data_left = 0.08
+            data_bottom = 0.12
+            data_width_frac = 0.74
+            data_height_frac = 0.78
+
+            ax = fig.add_axes([data_left, data_bottom, data_width_frac, data_height_frac])
             if is_rgba or is_rgb:
                 im = ax.imshow(data, aspect='equal', origin='upper',
                                interpolation='bilinear')
@@ -380,22 +376,34 @@ def tif_to_png(tif_file, vis_dir, resolution):
 
             ax.set_title(f"{title}\n{description}", fontsize=15, fontweight='bold', pad=10)
 
-            if not is_rgb:
-                if is_rgba and saved_cmap is not None:
-                    # Create a ScalarMappable for the colorbar from the saved colormap
-                    sm = plt.cm.ScalarMappable(cmap=saved_cmap,
-                                               norm=plt.Normalize(vmin=saved_vmin, vmax=saved_vmax))
-                    sm.set_array([])
-                    cbar = plt.colorbar(sm, ax=ax, pad=0.02, shrink=0.85, aspect=30)
-                else:
-                    cbar = plt.colorbar(im, ax=ax, pad=0.02, shrink=0.85, aspect=30)
+            # Colorbar/legend area — always at the same position for consistent layout
+            cbar_left = data_left + data_width_frac + 0.02
+            cbar_width = 0.04
+            if is_rgb:
+                # RGB: descriptive text label instead of gradient colorbar
+                cbar_ax = fig.add_axes([cbar_left, data_bottom, cbar_width, data_height_frac])
+                cbar_ax.set_xticks([])
+                cbar_ax.set_yticks([])
+                cbar_ax.text(0.5, 0.5, legend_label, transform=cbar_ax.transAxes,
+                             fontsize=9, fontweight='bold', rotation=90,
+                             verticalalignment='center', horizontalalignment='center',
+                             wrap=True)
+                cbar_ax.set_frame_on(False)
+            elif is_rgba and saved_cmap is not None:
+                cbar_ax = fig.add_axes([cbar_left, data_bottom, cbar_width, data_height_frac])
+                sm = plt.cm.ScalarMappable(cmap=saved_cmap,
+                                           norm=plt.Normalize(vmin=saved_vmin, vmax=saved_vmax))
+                sm.set_array([])
+                cbar = plt.colorbar(sm, cax=cbar_ax)
                 cbar.ax.tick_params(labelsize=9, width=1.5)
                 cbar.outline.set_linewidth(1.5)
                 cbar.set_label(legend_label, fontsize=10, fontweight='bold')
             else:
-                ax.text(1.02, 0.5, legend_label, transform=ax.transAxes,
-                        fontsize=10, fontweight='bold', rotation=90,
-                        verticalalignment='center', horizontalalignment='left')
+                cbar_ax = fig.add_axes([cbar_left, data_bottom, cbar_width, data_height_frac])
+                cbar = plt.colorbar(im, cax=cbar_ax)
+                cbar.ax.tick_params(labelsize=9, width=1.5)
+                cbar.outline.set_linewidth(1.5)
+                cbar.set_label(legend_label, fontsize=10, fontweight='bold')
 
             # GPS coordinate ticks
             if gps_coords and 'x_ticks' in gps_coords:
@@ -444,8 +452,8 @@ def tif_to_png(tif_file, vis_dir, resolution):
             north_ax.text(0.5, 0.95, 'N', ha='center', va='top',
                          fontsize=13, fontweight='bold', color='black', zorder=11)
 
-            # Bottom info bar
-            info_ax = fig.add_subplot(gs[1])
+            # Bottom info bar — fixed position
+            info_ax = fig.add_axes([data_left, 0.02, data_width_frac + cbar_width + 0.02, 0.07])
             info_ax.axis('off')
 
             extent_km_x = (max_x - min_x) / 1000
@@ -496,12 +504,11 @@ def tif_to_png(tif_file, vis_dir, resolution):
 
             fig.patch.set_facecolor('white')
 
-            # Save as PNG then convert to WebP — use higher DPI for large data
+            # Save as PNG then convert to WebP — fixed layout, no bbox_inches='tight'
             save_dpi = 200 if width > 3000 else 150
             png_temp = vis_dir / f"{tif_file.stem}_temp.png"
             try:
-                plt.savefig(png_temp, dpi=save_dpi, bbox_inches='tight', pad_inches=0.15,
-                           facecolor='white', format='png')
+                plt.savefig(png_temp, dpi=save_dpi, facecolor='white', format='png')
             finally:
                 plt.close()
 
@@ -509,8 +516,9 @@ def tif_to_png(tif_file, vis_dir, resolution):
             img.save(str(webp_file), format='WEBP', lossless=True)
             png_temp.unlink(missing_ok=True)
 
-            # Delete source TIFF
-            tif_file.unlink(missing_ok=True)
+            # Delete source TIFF (unless --keep-tif)
+            if not keep_tif:
+                tif_file.unlink(missing_ok=True)
 
             return webp_file
 
@@ -565,7 +573,7 @@ def generate_pdf_report(basename, vis_dir, pdf_dir, resolution):
 
     # Sort analysis files by archaeological priority
     order = ['mslrm', 'svf', 'negative_openness',
-              'positive_openness', 'sailore', 'depressions', 'hillshade_multi',
+              'positive_openness', 'sailore', 'hillshade_multi',
               'lrm', 'tpi', 'slope', 'curvature', 'aspect',
               'roughness', 'anomalies', 'wavelet', 'flow']
 
diff --git a/lidar_pipeline/tests/test_cli.py b/lidar_pipeline/tests/test_cli.py
index 7fdf8fb..80faea3 100644
--- a/lidar_pipeline/tests/test_cli.py
+++ b/lidar_pipeline/tests/test_cli.py
@@ -17,6 +17,7 @@ class TestCLIParsing:
         parser.add_argument("-w", "--workers", type=int, default=1)
         parser.add_argument("-f", "--force", action="store_true")
         parser.add_argument("--file", nargs="+", type=str, default=None)
+        parser.add_argument("--keep-tif", action="store_true")
 
         args = parser.parse_args(["./input"])
         assert args.input == "./input"
@@ -25,6 +26,7 @@ class TestCLIParsing:
         assert args.workers == 1
         assert args.force is False
         assert args.file is None
+        assert args.keep_tif is False
 
     def test_file_flag_single(self):
         import argparse
@@ -65,7 +67,7 @@ class TestSetupLogging:
         assert len(logger.handlers) == 1
         # Format should not include timestamps
         fmt = logger.handlers[0].formatter._fmt
-        assert "%(asctime)" not in fmt
+        assert "%(asctime)s" not in fmt
 
     def test_verbose_logging(self):
         """Verbose logging includes timestamps."""
@@ -73,7 +75,7 @@ class TestSetupLogging:
         from lidar_pipeline.cli import setup_logging
         logger = setup_logging(verbose=True, debug=False)
         fmt = logger.handlers[0].formatter._fmt
-        assert "%(asctime)" in fmt
+        assert "%(asctime)s" in fmt
 
     def test_debug_logging(self):
         """Debug logging includes file:line info."""
@@ -82,5 +84,5 @@ class TestSetupLogging:
         logger = setup_logging(verbose=False, debug=True)
         assert logger.level == logging.DEBUG
         fmt = logger.handlers[0].formatter._fmt
-        assert "%(filename)" in fmt
-        assert "%(lineno)" in fmt
\ No newline at end of file
+        assert "%(filename)s" in fmt
+        assert "%(lineno)d" in fmt
\ No newline at end of file
diff --git a/lidar_pipeline/tests/test_gpu.py b/lidar_pipeline/tests/test_gpu.py
index a333e08..b4f22bf 100644
--- a/lidar_pipeline/tests/test_gpu.py
+++ b/lidar_pipeline/tests/test_gpu.py
@@ -11,12 +11,12 @@ def test_has_gpu_attribute():
 
 
 def test_to_gpu_returns_array():
-    """to_gpu returns a float64 array with correct values."""
+    """to_gpu returns a float32 array with correct values."""
     from lidar_pipeline.gpu import to_gpu, to_cpu, HAS_GPU
     arr = np.array([1.0, 2.0, 3.0])
     result = to_gpu(arr)
-    # On GPU: cupy.ndarray, on CPU: numpy.ndarray
-    assert result.dtype == np.float64
+    # to_gpu converts to float32 to reduce GPU memory usage
+    assert result.dtype == np.float32
     # Always bring back to CPU for comparison
     np.testing.assert_array_equal(to_cpu(result), [1.0, 2.0, 3.0])
 
diff --git a/lidar_pipeline/tests/test_pipeline.py b/lidar_pipeline/tests/test_pipeline.py
index af01343..adcdd59 100644
--- a/lidar_pipeline/tests/test_pipeline.py
+++ b/lidar_pipeline/tests/test_pipeline.py
@@ -26,9 +26,9 @@ class TestVizSteps:
         assert "solar" not in names
 
     def test_expected_visualization_count(self):
-        """Should have 19 visualizations (18 terrain + ortho + topo - solar)."""
+        """Should have 17 visualizations (15 terrain + ortho + topo)."""
         from lidar_pipeline.pipeline import VIZ_STEPS
-        assert len(VIZ_STEPS) == 19
+        assert len(VIZ_STEPS) == 17
 
     def test_ortho_and_topo_present(self):
         from lidar_pipeline.pipeline import VIZ_STEPS
diff --git a/lidar_pipeline/tests/test_visualizations.py b/lidar_pipeline/tests/test_visualizations.py
index ffe1647..177515d 100644
--- a/lidar_pipeline/tests/test_visualizations.py
+++ b/lidar_pipeline/tests/test_visualizations.py
@@ -141,13 +141,6 @@ class TestTPI:
         assert result.exists()
 
 
-class TestDepressions:
-    def test_generates_tif(self, synthetic_dem, tmp_output_dir):
-        from lidar_pipeline.visualizations import generate_depressions
-        result = generate_depressions(synthetic_dem, "test", tmp_output_dir, 5.0)
-        assert result is not None
-        assert result.exists()
-
 
 class TestSAILORE:
     def test_generates_tif(self, synthetic_dem, tmp_output_dir):
diff --git a/lidar_pipeline/visualizations.py b/lidar_pipeline/visualizations.py
index a22cc93..e6cbc25 100644
--- a/lidar_pipeline/visualizations.py
+++ b/lidar_pipeline/visualizations.py
@@ -117,7 +117,12 @@ def _filter_nanaware(arr, filter_func, *args, use_gpu=True, **kwargs):
 # ============================================================
 
 def generate_hillshade(dem_file, basename, vis_dir, resolution):
-    """Generate multi-directional hillshade (NW, NE, SW, SE) — GPU if available."""
+    """Generate multi-directional hillshade with slope shading — GPU if available.
+
+    Combines 4-direction hillshade (NW, NE, SW, SE) with slope shading
+    for improved micro-relief visibility on flat terrain.
+    Result = 0.7 * hillshade + 0.3 * cos(slope).
+    """
     gpu_tag = " [GPU]" if HAS_GPU else ""
     logger.info(f"  → Hillshade multidirectionnel{gpu_tag}...")
     t0 = time.time()
@@ -146,7 +151,10 @@ def generate_hillshade(dem_file, basename, vis_dir, resolution):
             hs = sin_alt * sin_slope + cos_alt * cos_slope * xp.cos(az_rad - aspect)
             hillshades.append(xp.clip(hs, 0, 1))
 
-        combined = xp.mean(xp.array(hillshades), axis=0)
+        combined_hillshade = xp.mean(xp.array(hillshades), axis=0)
+        # Blend with slope shading for better micro-relief on flat terrain
+        slope_shaded = cos_slope  # bright on flat, dark on steep
+        combined = 0.7 * combined_hillshade + 0.3 * slope_shaded
         _save_tif(output, to_cpu(combined), transform, crs)
         logger.info(f"  ✓ Hillshade terminé ({time.time()-t0:.1f}s){gpu_tag}")
         return output
@@ -240,8 +248,6 @@ def generate_lrm(dem_file, basename, vis_dir, resolution):
         _save_tif(output, lrm.astype(np.float32), transform, crs)
         logger.info(f"  ✓ LRM terminé ({time.time()-t0:.1f}s){gpu_tag}")
         return output
-        logger.info(f"  ✓ LRM terminé ({time.time()-t0:.1f}s){gpu_tag}")
-        return output
     except Exception as e:
         logger.error(f"  ✗ Erreur LRM: {e}", exc_info=True)
         return None
@@ -279,13 +285,18 @@ def generate_svf(dem_file, basename, vis_dir, resolution):
             ddx, ddy = dx[d_idx], dy[d_idx]
             horizon = xp.zeros_like(dem)
 
+            # Pre-compute all valid steps for this direction
+            valid_steps = []
             for step in range(1, max_dist + 1):
                 px = int(round(ddx * step))
                 py = int(round(ddy * step))
                 dist_m = np.sqrt((ddx * step * res) ** 2 + (ddy * step * res) ** 2)
                 if dist_m < res * 0.5:
                     continue
+                valid_steps.append((step, px, py, dist_m))
 
+            # Batch all shifts into a single array for vectorized max computation
+            for step, px, py, dist_m in valid_steps:
                 elev_diff = padded[max_dist + py:max_dist + py + rows,
                                   max_dist + px:max_dist + px + cols] - dem
                 angle = xp.arctan2(elev_diff, dist_m)
@@ -447,55 +458,6 @@ def generate_tpi(dem_file, basename, vis_dir, resolution):
         return None
 
 
-# ============================================================
-# Depression / hydrology
-# ============================================================
-
-def generate_depressions(dem_file, basename, vis_dir, resolution):
-    """Depression detection using hydrological sink filling — GPU if available."""
-    gpu_tag = " [GPU]" if HAS_GPU else ""
-    logger.info(f"  → Détection dépressions (hydrologique){gpu_tag}...")
-    t0 = time.time()
-    output = vis_dir / f"{basename}_depressions.tif"
-
-    try:
-        dem_np, transform, crs = _read_dem(dem_file)
-        dem = to_gpu(dem_np)
-
-        from scipy.ndimage import generate_binary_structure
-        struct = generate_binary_structure(2, 2)
-
-        dem_filled = xp.copy(dem)
-        nodata_mask = xp.isnan(dem_filled)
-        dem_filled[nodata_mask] = xp.nanmax(dem) + 1000
-
-        changed = True
-        iterations = 0
-        max_iter = 100
-
-        while changed and iterations < max_iter:
-            neighbor_min = xp_minimum_filter(dem_filled, footprint=struct)
-            sinks = (dem_filled < neighbor_min) & ~nodata_mask
-
-            if not xp.any(sinks):
-                break
-
-            new_dem = xp.maximum(dem_filled, neighbor_min)
-            new_dem[nodata_mask] = xp.nan
-            changed = bool(xp.any(new_dem != dem_filled))
-            dem_filled = new_dem
-            iterations += 1
-
-        depressions = to_cpu(dem_filled - dem)
-        depressions[to_cpu(nodata_mask)] = np.nan
-        depressions = np.where(depressions > 0.01, depressions, 0)
-
-        _save_tif(output, depressions, transform, crs)
-        logger.info(f"  ✓ Dépressions terminé ({time.time()-t0:.1f}s){gpu_tag}")
-        return output
-    except Exception as e:
-        logger.error(f"  ✗ Erreur dépressions: {e}", exc_info=True)
-        return None
 
 
 # ============================================================
@@ -680,8 +642,63 @@ def generate_wavelet(dem_file, basename, vis_dir, resolution):
 # Flow accumulation
 # ============================================================
 
+def _d8_accumulate_numba(flow_dir, nodata_mask, rows, cols):
+    """JIT-compiled D8 flow accumulation loop.
+
+    Uses numba for ~100x speedup over pure Python loop.
+    Falls back to pure Python if numba is unavailable.
+    """
+    try:
+        from numba import njit
+
+        @njit(cache=True)
+        def _accumulate(flow_dir, nodata_mask, rows, cols):
+            dx8 = np.array([1, 1, 0, -1, -1, -1, 0, 1], dtype=np.int8)
+            dy8 = np.array([0, 1, 1, 1, 0, -1, -1, -1], dtype=np.int8)
+
+            flow_acc = np.ones((rows, cols), dtype=np.float32)
+
+            # Sort cells by elevation (high to low) — walk downhill
+            # We use the fact that flow_dir already encodes steepest descent
+            # Process from highest to lowest elevation
+            for r in range(rows):
+                for c in range(cols):
+                    if nodata_mask[r, c]:
+                        flow_acc[r, c] = 0.0
+                        continue
+
+            # Iterative accumulation: process cells in top-down order
+            # Multiple passes until convergence
+            for _pass in range(10):
+                changed = 0
+                for r in range(rows):
+                    for c in range(cols):
+                        if nodata_mask[r, c]:
+                            continue
+                        d = flow_dir[r, c]
+                        if d < 0:
+                            continue
+                        nr = r + dy8[d]
+                        nc = c + dx8[d]
+                        if 0 <= nr < rows and 0 <= nc < cols and not nodata_mask[nr, nc]:
+                            old_acc = flow_acc[nr, nc]
+                            flow_acc[nr, nc] += flow_acc[r, c]
+                            if flow_acc[nr, nc] != old_acc:
+                                changed += 1
+                if changed == 0:
+                    break
+
+            return flow_acc
+
+        return _accumulate(flow_dir, nodata_mask, rows, cols)
+
+    except ImportError:
+        # Fallback: pure Python
+        return None
+
+
 def generate_flow(dem_file, basename, vis_dir, resolution):
-    """Flow accumulation using D8 algorithm — sink filling on GPU, accumulation on CPU."""
+    """Flow accumulation using D8 algorithm — sink filling on GPU, accumulation via numba."""
     gpu_tag = " [GPU]" if HAS_GPU else ""
     logger.info(f"  → Accumulation de flux D8{gpu_tag}...")
     t0 = time.time()
@@ -711,10 +728,10 @@ def generate_flow(dem_file, basename, vis_dir, resolution):
         dem_filled[nodata_mask_gpu] = xp.nan
         dem_filled_np = to_cpu(dem_filled)
 
-        # D8 slope + accumulation — CPU (sequential by nature)
-        dx8 = [1, 1, 0, -1, -1, -1, 0, 1]
-        dy8 = [0, 1, 1, 1, 0, -1, -1, -1]
-        dist8 = [1.0, np.sqrt(2), 1.0, np.sqrt(2), 1.0, np.sqrt(2), 1.0, np.sqrt(2)]
+        # D8 slope — vectorized
+        dx8 = np.array([1, 1, 0, -1, -1, -1, 0, 1], dtype=np.int32)
+        dy8 = np.array([0, 1, 1, 1, 0, -1, -1, -1], dtype=np.int32)
+        dist8 = np.array([1.0, np.sqrt(2), 1.0, np.sqrt(2), 1.0, np.sqrt(2), 1.0, np.sqrt(2)])
 
         flow_dir = np.full((rows, cols), -1, dtype=np.int8)
         max_slope = np.zeros((rows, cols), dtype=np.float64)
@@ -732,21 +749,30 @@ def generate_flow(dem_file, basename, vis_dir, resolution):
             flow_dir[better] = d
             max_slope[better] = slope[better]
 
-        flat_dem = dem_filled_np[~nodata_mask].flatten()
-        valid_indices = np.where(~nodata_mask.flatten())[0]
-        sort_order = valid_indices[np.argsort(-flat_dem)]
+        # D8 accumulation — try numba first, fallback to Python
+        result = _d8_accumulate_numba(flow_dir, nodata_mask.astype(np.bool_), rows, cols)
 
-        flow_acc = np.ones((rows, cols), dtype=np.float32)
-        flow_acc[nodata_mask] = 0
+        if result is not None:
+            flow_acc = result
+            logger.info(f"    Accumulation D8 via numba")
+        else:
+            # Pure Python fallback (slow for large DEMs)
+            logger.info(f"    Accumulation D8 via Python (installez numba pour accélérer)")
+            flat_dem = dem_filled_np[~nodata_mask].flatten()
+            valid_indices = np.where(~nodata_mask.flatten())[0]
+            sort_order = valid_indices[np.argsort(-flat_dem)]
 
-        for idx in sort_order:
-            r, c = divmod(idx, cols)
-            d = flow_dir[r, c]
-            if d < 0:
-                continue
-            nr, nc = r + dy8[d], c + dx8[d]
-            if 0 <= nr < rows and 0 <= nc < cols and not nodata_mask[nr, nc]:
-                flow_acc[nr, nc] += flow_acc[r, c]
+            flow_acc = np.ones((rows, cols), dtype=np.float32)
+            flow_acc[nodata_mask] = 0
+
+            for idx in sort_order:
+                r, c = divmod(idx, cols)
+                d = flow_dir[r, c]
+                if d < 0:
+                    continue
+                nr, nc = r + dy8[d], c + dx8[d]
+                if 0 <= nr < rows and 0 <= nc < cols and not nodata_mask[nr, nc]:
+                    flow_acc[nr, nc] += flow_acc[r, c]
 
         flow_log = np.log1p(flow_acc)
         _save_tif(output, flow_log, transform, crs)
diff --git a/run.sh b/run.sh
index cd0c13c..9540727 100755
--- a/run.sh
+++ b/run.sh
@@ -9,6 +9,7 @@
 #   -v              Mode verbeux (timestamps + niveaux)
 #   --debug         Mode debug (détails internes fichier:ligne)
 #   -f / --force    Régénérer tous les fichiers même si existants
+#   --keep-tif       Conserver les fichiers TIFF intermédiaires
 #   --force-classification
 #                   Reclassifier le sol même si le fichier .las existe déjà
 #   --ground-classification {auto,smrf,pmf,csf}
@@ -33,6 +34,7 @@ if [ $# -eq 0 ]; then
     echo "  -f / --force    Régénérer tous les fichiers même si les WebP existent"
     echo "  --force-classification"
     echo "                  Reclassifier le sol même si le fichier .las existe"
+    echo "  --keep-tif      Conserver les fichiers TIFF intermédiaires"
     echo "  --ground-classification {auto,smrf,pmf,csf}"
     echo "                  Méthode de classification du sol (défaut: auto)"
     echo "  --file NOM...   Traiter un ou plusieurs fichiers LAZ (nom complet sans .laz/.las)"
@@ -64,6 +66,7 @@ FORCE_FLAG=""
 FILE_ARGS=""
 GROUND_METHOD=""
 FORCE_CLASSIFY_FLAG=""
+KEEP_TIF_FLAG=""
 
 # Parse arguments manually (more robust than getopts for mixed short/long options)
 while [ $# -gt 0 ]; do
@@ -76,6 +79,7 @@ while [ $# -gt 0 ]; do
         --debug) VERBOSE_FLAG="--debug"; shift ;;
         --force) FORCE_FLAG="--force"; shift ;;
         --force-classification) FORCE_CLASSIFY_FLAG="--force-classification"; shift ;;
+        --keep-tif) KEEP_TIF_FLAG="--keep-tif"; shift ;;
         --ground-classification) GROUND_METHOD="$2"; shift 2 ;;
         --ground-classification=*) GROUND_METHOD="${1#--ground-classification=}"; shift ;;
         --file) shift; while [ $# -gt 0 ] && [[ ! "$1" =~ ^- ]]; do FILE_ARGS="$FILE_ARGS $1"; shift; done ;;
@@ -93,6 +97,7 @@ while [ $# -gt 0 ]; do
             echo "  -f / --force    Régénérer tous les fichiers même si les WebP existent"
             echo "  --force-classification"
             echo "                  Reclassifier le sol même si le fichier .las existe"
+            echo "  --keep-tif      Conserver les fichiers TIFF intermédiaires"
             echo "  --ground-classification {auto,smrf,pmf,csf}"
             echo "                  Méthode de classification du sol (défaut: auto)"
             echo "  --file NOM...   Traiter un ou plusieurs fichiers LAZ (nom complet sans .laz/.las)"
@@ -153,13 +158,14 @@ echo "  GPU        : $([ -n "$GPU_FLAG" ] && echo 'OUI' || echo 'non')"
 echo "  Verbeux    : $([ -n "$VERBOSE_FLAG" ] && echo 'OUI' || echo 'non')"
 echo "  Force      : $([ -n "$FORCE_FLAG" ] && echo 'OUI' || echo 'non')"
 echo "  Force classif.: $([ -n "$FORCE_CLASSIFY_FLAG" ] && echo 'OUI' || echo 'non')"
+echo "  Keep TIFF    : $([ -n "$KEEP_TIF_FLAG" ] && echo 'OUI' || echo 'non')"
 echo "  Classification sol : $([ -n "$GROUND_METHOD" ] && echo "$GROUND_METHOD" || echo 'auto')"
 if [ -n "$FILE_ARGS" ]; then
     echo "  Fichiers   :${FILE_ARGS}"
 fi
 echo "============================================"
 
-CMD_ARGS="-o /data/output -r $RESOLUTION -w $WORKERS $VERBOSE_FLAG $FORCE_FLAG $FORCE_CLASSIFY_FLAG"
+CMD_ARGS="-o /data/output -r $RESOLUTION -w $WORKERS $VERBOSE_FLAG $FORCE_FLAG $FORCE_CLASSIFY_FLAG $KEEP_TIF_FLAG"
 if [ -n "$GROUND_METHOD" ]; then
     CMD_ARGS="$CMD_ARGS --ground-classification $GROUND_METHOD"
 fi