diff --git a/CHANGELOG.md b/CHANGELOG.md
index 4436444..d347da2 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -7,6 +7,22 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ## [Unreleased]
 
+### Added
+
+- Optional terminal soft-clip artifact model for `simulate`. Real libraries
+  soft-clip a few bases off one or both read ends on a meaningful fraction of
+  reads (end-repair fill-in, damaged or non-templated ends); Twist EM-seq, for
+  example, shows ~8% of reads clipped, strongly 5'-biased and modal at 6–10 bp.
+  Four new flags control it: `--clip-5p-rate` and `--clip-3p-rate` (both
+  default `0.0`, i.e. disabled) give independent per-end probabilities, and
+  `--clip-length-mean` (default `8`) / `--clip-length-max` (default `20`) shape
+  a shared truncated-geometric clip-length distribution. Clipped bases are
+  replaced with random sequence and the soft-clip is recorded in the
+  golden-BAM CIGAR (advancing the alignment start when a forward read's 5' end
+  is clipped), so ground truth stays exact. The model is protocol-agnostic
+  (independent of `--methylation-mode`); a disabled or omitted model draws no
+  randomness, leaving existing seeded runs byte-identical.
+
 ## [0.3.0] - 2026-06-13
 
 ### Added
diff --git a/CLAUDE.md b/CLAUDE.md
index 186cfc5..2b897ce 100644
--- a/CLAUDE.md
+++ b/CLAUDE.md
@@ -40,6 +40,7 @@ Holodeck is a single-crate Rust project with a binary (`holodeck`) and library (
 | `vcf/methylation.rs` | MT/MB FORMAT field read/write; CpG classification across haplotypes |
 | `output/cpg_truth.rs` | Coverage-weighted per-CpG truth bedGraph (from simulated reads) |
 | `output/methylation_bedgraph.rs` | Closed-form population-fraction methylation bedGraph |
+| `clip.rs` | Terminal soft-clip artifact model (per-end 5'/3' clip sampling) |
 | `fragment.rs` | Fragment extraction, reverse complement, adapter padding |
 | `read.rs` | Read pair generation combining fragments + error model + naming |
 | `error_model/mod.rs` | ErrorModel trait + apply_errors free function |
diff --git a/README.md b/README.md
index 70276f9..611a806 100644
--- a/README.md
+++ b/README.md
@@ -125,6 +125,10 @@ holodeck simulate -r ref.fa -v methylated.vcf.gz -o output_taps \
 | `--min-error-rate` | 0.001 | Error rate at start of reads |
 | `--max-error-rate` | 0.01 | Error rate at end of reads |
 | `--max-n-frac` | 0.02 | Reject reads with >this fraction of bases from ambiguous reference positions (see [Ambiguous reference bases](#ambiguous-reference-bases)) |
+| `--clip-5p-rate` | 0.0 | Probability a read's 5' end carries a terminal soft-clip artifact (see [Terminal soft-clip artifacts](#terminal-soft-clip-artifacts)) |
+| `--clip-3p-rate` | 0.0 | Probability a read's 3' end carries a terminal soft-clip artifact, independent of 3' adapter read-through |
+| `--clip-length-mean` | 8 | Mean injected clip length (truncated geometric); only used when a clip rate is non-zero |
+| `--clip-length-max` | 20 | Maximum injected clip length; only used when a clip rate is non-zero |
 | `--methylation-mode` | none | Methylation chemistry: `em-seq` (or `bisulfite`) or `taps`. Presence enables methylation simulation |
 | `--methylation-conversion-rate` | 0.999 | Chemistry efficiency for normally-converting molecules (probability that the converting class of C converts to T) |
 | `--methylation-failure-rate` | 0.01 | Fraction of molecules that are whole-molecule conversion failures (convert at `1 − conversion-rate`). Requires `--methylation-mode` |
@@ -155,6 +159,14 @@ Real references contain a mix of `A`/`C`/`G`/`T`, large stretches of `N` (assemb
 
 **Known limitation:** requested `--coverage` is computed from raw contig/BED lengths, not from the non-ambiguous territory. For a reference like hs38DH (~5% N), rejection is noise and coverage lands where you'd expect. For simulations targeted at heavily-N contigs (or with `--max-n-frac 0.0` in N-dense regions), effective coverage will be slightly below the requested value; a warning is logged if the resampling budget is exhausted.
 
+### Terminal soft-clip artifacts
+
+By default holodeck's reads align essentially end-to-end: aside from VCF-derived indels and 3' adapter read-through, the simulated CIGAR is one long match block. Real libraries are messier — end-repair fill-in, damaged or non-templated read ends, and similar effects make aligners soft-clip a few bases off one or both ends of a meaningful fraction of reads. Twist EM-seq, for example, shows ~8% of reads carrying a soft clip, strongly 5'-biased and modal at 6–10 bp. The mechanism is protocol-agnostic, so the model is independent of `--methylation-mode` and useful for both methylation and non-methylation simulation.
+
+`--clip-5p-rate` and `--clip-3p-rate` turn this on (both default to `0.0`, i.e. disabled). For each read, the 5' and 3' ends independently receive a clip with their configured probabilities; a clipped end's length is drawn from a truncated geometric distribution with mean `--clip-length-mean`, capped at `--clip-length-max`. The 5' and 3' rates are separate because real libraries are markedly asymmetric (the 5' end clips several times more often), while the length distribution is shared across both ends.
+
+When an end is clipped, holodeck replaces those terminal bases with random sequence — so a downstream aligner re-derives the soft-clip rather than forcing a match — and records the soft-clip in the golden-BAM CIGAR (and advances the alignment start when a forward read's 5' end is clipped). Ground truth therefore stays exact: the aligned block still maps to the correct reference span. A disabled (or omitted) model draws no randomness, so existing seeded runs remain byte-identical.
+
 ### Golden BAM
 
 Pass `--golden-bam` to write a perfect-truth BAM alongside the FASTQ output. Every record carries the alignment the simulator generated it from — contig, position, strand, CIGAR — at MAPQ 60, in unsorted (generation) order. Sort downstream if a coordinate-sorted BAM is needed.
diff --git a/src/clip.rs b/src/clip.rs
new file mode 100644
index 0000000..5d54fd5
--- /dev/null
+++ b/src/clip.rs
@@ -0,0 +1,224 @@
+//! Terminal soft-clip artifact model.
+//!
+//! Real sequencing libraries frequently produce reads whose extreme 5' or 3'
+//! ends diverge from the reference enough that an aligner soft-clips them
+//! (end-repair fill-in, damaged or otherwise non-templated read ends). Twist
+//! EM-seq, for example, shows ~8% of reads carrying a soft clip, strongly
+//! 5'-biased and modal at 6-10 bp.
+//!
+//! Holodeck's simulated reads otherwise align essentially end-to-end, so this
+//! module injects that artifact: for a configurable fraction of reads it
+//! marks the first/last `K` bases of the *aligned* portion as soft-clipped.
+//! The caller corrupts those bases with random sequence (so a downstream
+//! aligner re-derives the clip) and records the soft-clip in the golden BAM
+//! CIGAR, keeping ground truth exact.
+//!
+//! The mechanism is protocol-agnostic — it is independent of methylation
+//! chemistry and useful for any library type.
+
+use rand::Rng;
+
+/// Configuration for the terminal soft-clip artifact model.
+///
+/// The 5' and 3' ends clip independently, each with its own probability,
+/// because real libraries are markedly asymmetric (the 5' end clips several
+/// times more often than the 3' end). The clip-length distribution is shared
+/// across both ends, as observed lengths do not differ meaningfully between
+/// them.
+#[derive(Debug, Clone, Copy)]
+pub struct TerminalClipConfig {
+    /// Probability in `[0.0, 1.0]` that a read's 5' end receives a clip.
+    pub rate_5p: f64,
+    /// Probability in `[0.0, 1.0]` that a read's 3' end receives a clip.
+    pub rate_3p: f64,
+    /// Mean clip length in bases. Lengths are drawn from a truncated
+    /// geometric distribution with this mean.
+    pub length_mean: usize,
+    /// Maximum clip length in bases; sampled lengths are clamped to this.
+    pub length_max: usize,
+}
+
+impl TerminalClipConfig {
+    /// Whether either end can produce a clip. When this is `false`,
+    /// [`Self::sample_clips`] consumes no randomness and returns `(0, 0)`, so
+    /// the simulator's output is byte-identical to having no clip model at
+    /// all.
+    #[must_use]
+    pub fn is_enabled(&self) -> bool {
+        self.rate_5p > 0.0 || self.rate_3p > 0.0
+    }
+
+    /// Sample `(clip_5p, clip_3p)` lengths for one read whose aligned
+    /// (genomic) portion is `aligned_len` bases long.
+    ///
+    /// Each end independently clips with its configured probability; a chosen
+    /// clip's length is drawn from a truncated geometric with mean
+    /// `length_mean`, clamped to `[1, length_max]`. The combined clip is then
+    /// capped so at least one aligned base survives (an alignment needs a
+    /// match block), favouring the 5' end when both are drawn and the budget
+    /// is tight.
+    ///
+    /// Returns `(0, 0)` without drawing any randomness when the model is
+    /// disabled or `aligned_len < 2`.
+    pub fn sample_clips(&self, aligned_len: usize, rng: &mut impl Rng) -> (usize, usize) {
+        // No clip is representable without leaving an aligned base behind, and
+        // a disabled model must not perturb the RNG stream.
+        if !self.is_enabled() || aligned_len < 2 {
+            return (0, 0);
+        }
+
+        let mut clip_5p = self.sample_end(self.rate_5p, rng);
+        let mut clip_3p = self.sample_end(self.rate_3p, rng);
+
+        // Cap so at least one aligned base remains. Trim the 3' end first
+        // (the rarer, less characteristic artifact), then the 5' end.
+        let budget = aligned_len - 1;
+        if clip_5p + clip_3p > budget {
+            clip_3p = clip_3p.min(budget.saturating_sub(clip_5p));
+            clip_5p = clip_5p.min(budget);
+        }
+
+        (clip_5p, clip_3p)
+    }
+
+    /// Draw a clip length for a single end: `0` if the end does not clip,
+    /// otherwise a truncated-geometric length in `[1, length_max]`. Consumes
+    /// no randomness when `rate <= 0.0`.
+    fn sample_end(&self, rate: f64, rng: &mut impl Rng) -> usize {
+        if rate <= 0.0 {
+            return 0;
+        }
+        if rng.random::<f64>() >= rate {
+            return 0;
+        }
+        self.sample_length(rng)
+    }
+
+    /// Sample a clip length from a truncated geometric distribution with mean
+    /// [`Self::length_mean`], clamped to `[1, length_max]`.
+    fn sample_length(&self, rng: &mut impl Rng) -> usize {
+        let max = self.length_max.max(1);
+        let mean = self.length_mean.max(1);
+        if mean <= 1 {
+            // Degenerate distribution: always the shortest representable clip.
+            return 1;
+        }
+        #[expect(clippy::cast_precision_loss, reason = "small clip means")]
+        let p = 1.0 / mean as f64;
+        let u: f64 = rng.random(); // [0, 1)
+        // Inverse-CDF of the geometric (trials >= 1): k = floor(ln(1-u)/ln(1-p)) + 1.
+        let k = ((1.0 - u).ln() / (1.0 - p).ln()).floor();
+        #[expect(
+            clippy::cast_possible_truncation,
+            clippy::cast_sign_loss,
+            reason = "k is finite, non-negative after floor"
+        )]
+        let k = k as usize + 1;
+        k.clamp(1, max)
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use rand::SeedableRng;
+    use rand::rngs::SmallRng;
+
+    fn cfg(rate_5p: f64, rate_3p: f64) -> TerminalClipConfig {
+        TerminalClipConfig { rate_5p, rate_3p, length_mean: 8, length_max: 20 }
+    }
+
+    #[test]
+    fn disabled_config_is_not_enabled() {
+        assert!(!cfg(0.0, 0.0).is_enabled());
+        assert!(cfg(0.1, 0.0).is_enabled());
+        assert!(cfg(0.0, 0.1).is_enabled());
+    }
+
+    #[test]
+    fn disabled_config_returns_zero_and_consumes_no_randomness() {
+        // A disabled model must leave the RNG stream untouched so existing
+        // simulations stay byte-identical.
+        let mut rng_a = SmallRng::seed_from_u64(7);
+        let clips = cfg(0.0, 0.0).sample_clips(100, &mut rng_a);
+        assert_eq!(clips, (0, 0));
+        let after: u64 = rng_a.random();
+
+        let mut rng_b = SmallRng::seed_from_u64(7);
+        let direct: u64 = rng_b.random();
+        assert_eq!(after, direct, "disabled clip model must not advance the RNG");
+    }
+
+    #[test]
+    fn short_alignment_cannot_clip() {
+        let mut rng = SmallRng::seed_from_u64(1);
+        assert_eq!(cfg(1.0, 1.0).sample_clips(0, &mut rng), (0, 0));
+        assert_eq!(cfg(1.0, 1.0).sample_clips(1, &mut rng), (0, 0));
+    }
+
+    #[test]
+    fn certain_5p_clip_always_clips_5p_only() {
+        let mut rng = SmallRng::seed_from_u64(42);
+        for _ in 0..200 {
+            let (c5, c3) = cfg(1.0, 0.0).sample_clips(150, &mut rng);
+            assert!(c5 >= 1, "rate_5p=1.0 must always clip the 5' end");
+            assert_eq!(c3, 0, "rate_3p=0.0 must never clip the 3' end");
+            assert!(c5 <= 20, "clip must respect length_max");
+        }
+    }
+
+    #[test]
+    fn certain_3p_clip_always_clips_3p_only() {
+        let mut rng = SmallRng::seed_from_u64(42);
+        for _ in 0..200 {
+            let (c5, c3) = cfg(0.0, 1.0).sample_clips(150, &mut rng);
+            assert_eq!(c5, 0);
+            assert!((1..=20).contains(&c3));
+        }
+    }
+
+    #[test]
+    fn clip_never_consumes_whole_alignment() {
+        let mut rng = SmallRng::seed_from_u64(99);
+        // length_max larger than the alignment forces the cap to engage.
+        let config =
+            TerminalClipConfig { rate_5p: 1.0, rate_3p: 1.0, length_mean: 8, length_max: 50 };
+        for aligned in 2..30 {
+            for _ in 0..50 {
+                let (c5, c3) = config.sample_clips(aligned, &mut rng);
+                assert!(c5 + c3 <= aligned - 1, "must leave >=1 aligned base (aligned={aligned})");
+            }
+        }
+    }
+
+    #[test]
+    fn observed_rate_tracks_configured_rate() {
+        let mut rng = SmallRng::seed_from_u64(2024);
+        let mut clipped = 0;
+        let n = 20_000;
+        for _ in 0..n {
+            let (c5, _) = cfg(0.3, 0.0).sample_clips(150, &mut rng);
+            if c5 > 0 {
+                clipped += 1;
+            }
+        }
+        let observed = f64::from(clipped) / f64::from(n);
+        assert!((observed - 0.3).abs() < 0.02, "expected ~0.30 clip rate, got {observed}");
+    }
+
+    #[test]
+    fn mean_length_is_in_the_right_ballpark() {
+        let mut rng = SmallRng::seed_from_u64(2025);
+        let config =
+            TerminalClipConfig { rate_5p: 1.0, rate_3p: 0.0, length_mean: 8, length_max: 1000 };
+        let mut sum = 0usize;
+        let n = 50_000;
+        for _ in 0..n {
+            let (c5, _) = config.sample_clips(2000, &mut rng);
+            sum += c5;
+        }
+        #[expect(clippy::cast_precision_loss, reason = "test arithmetic")]
+        let mean = sum as f64 / f64::from(n);
+        assert!((mean - 8.0).abs() < 1.0, "expected mean ~8, got {mean}");
+    }
+}
diff --git a/src/commands/simulate.rs b/src/commands/simulate.rs
index 3bcb102..a7f0a41 100644
--- a/src/commands/simulate.rs
+++ b/src/commands/simulate.rs
@@ -123,6 +123,34 @@ pub struct Simulate {
     #[arg(long, default_value_t = 0.02, value_name = "FLOAT")]
     pub max_n_frac: f64,
 
+    /// Probability in `[0.0, 1.0]` that a read's 5' end carries a terminal
+    /// soft-clip artifact (end-repair fill-in, damaged or non-templated read
+    /// ends). The clipped bases are replaced with random sequence so a
+    /// downstream aligner soft-clips them, and the golden BAM records the clip
+    /// in the CIGAR so ground truth stays exact. Real Twist EM-seq shows ~8%
+    /// of reads 5'-clipped. Defaults to `0.0` (disabled). Protocol-agnostic:
+    /// independent of `--methylation-mode`.
+    #[arg(long, default_value_t = 0.0, value_name = "FLOAT")]
+    pub clip_5p_rate: f64,
+
+    /// Probability in `[0.0, 1.0]` that a read's 3' end carries a terminal
+    /// soft-clip artifact, independent of (and in addition to) 3' adapter
+    /// read-through. Real Twist EM-seq shows ~2% of reads 3'-clipped. Defaults
+    /// to `0.0` (disabled).
+    #[arg(long, default_value_t = 0.0, value_name = "FLOAT")]
+    pub clip_3p_rate: f64,
+
+    /// Mean length in bases of an injected terminal soft-clip. Lengths are
+    /// drawn from a truncated geometric distribution centered here and capped
+    /// by `--clip-length-max`. Only takes effect when a clip rate is non-zero.
+    #[arg(long, default_value_t = 8, value_name = "INT")]
+    pub clip_length_mean: usize,
+
+    /// Maximum length in bases of an injected terminal soft-clip. Only takes
+    /// effect when a clip rate is non-zero.
+    #[arg(long, default_value_t = 20, value_name = "INT")]
+    pub clip_length_max: usize,
+
     /// Enable methylation chemistry simulation. `em-seq` (or `bisulfite`)
     /// converts unmethylated cytosines to thymine and preserves methylated
     /// ones (matches both bisulfite and em-seq protocols). `taps` is the
@@ -223,6 +251,22 @@ impl Simulate {
         if !self.max_n_frac.is_finite() || !(0.0..=1.0).contains(&self.max_n_frac) {
             bail!("--max-n-frac must be in [0.0, 1.0]");
         }
+        if !self.clip_5p_rate.is_finite() || !(0.0..=1.0).contains(&self.clip_5p_rate) {
+            bail!("--clip-5p-rate must be in [0.0, 1.0]");
+        }
+        if !self.clip_3p_rate.is_finite() || !(0.0..=1.0).contains(&self.clip_3p_rate) {
+            bail!("--clip-3p-rate must be in [0.0, 1.0]");
+        }
+        // Validate clip-length parameters only when the model is enabled, so a
+        // disabled model (both rates zero) never rejects on unused values.
+        if self.clip_5p_rate > 0.0 || self.clip_3p_rate > 0.0 {
+            if self.clip_length_max == 0 {
+                bail!("--clip-length-max must be > 0 when a terminal clip rate is set");
+            }
+            if self.clip_length_mean > self.clip_length_max {
+                bail!("--clip-length-mean must be <= --clip-length-max");
+            }
+        }
         if let Some(rate) = self.methylation_conversion_rate
             && (!rate.is_finite() || !(0.0..=1.0).contains(&rate))
         {
@@ -349,6 +393,18 @@ impl Simulate {
 
         let error_model =
             IlluminaErrorModel::new(self.read_length, self.min_error_rate, self.max_error_rate);
+
+        // Terminal soft-clip artifact model. Passed as `Some` only when at
+        // least one end can clip; a disabled model is equivalent to `None`
+        // (it draws no randomness) but keeping it `None` makes intent clear.
+        let clip_config = (self.clip_5p_rate > 0.0 || self.clip_3p_rate > 0.0).then_some(
+            crate::clip::TerminalClipConfig {
+                rate_5p: self.clip_5p_rate,
+                rate_3p: self.clip_3p_rate,
+                length_mean: self.clip_length_mean,
+                length_max: self.clip_length_max,
+            },
+        );
         let frag_dist = Normal::new(self.fragment_mean as f64, self.fragment_stddev as f64)
             .map_err(|e| anyhow::anyhow!("Invalid fragment distribution parameters: {e}"))?;
 
@@ -457,6 +513,7 @@ impl Simulate {
                 targets.as_ref(),
                 total_reads,
                 &error_model,
+                clip_config.as_ref(),
                 &frag_dist,
                 adapter_r1.as_bytes(),
                 adapter_r2.as_bytes(),
@@ -599,6 +656,7 @@ impl Simulate {
         targets: Option<&TargetRegions>,
         total_reads: u64,
         error_model: &IlluminaErrorModel,
+        clip_config: Option<&crate::clip::TerminalClipConfig>,
         frag_dist: &Normal<f64>,
         adapter_r1: &[u8],
         adapter_r2: &[u8],
@@ -803,6 +861,7 @@ impl Simulate {
                 // YS:Z tag — skip the per-mate clone when no golden BAM is
                 // requested.
                 self.golden_bam,
+                clip_config,
                 rng,
             ) else {
                 // Too many ambiguity-resolved bases in this read pair; resample.
@@ -926,6 +985,10 @@ mod tests {
             min_error_rate: 0.001,
             max_error_rate: 0.01,
             max_n_frac: 0.02,
+            clip_5p_rate: 0.0,
+            clip_3p_rate: 0.0,
+            clip_length_mean: 8,
+            clip_length_max: 20,
             methylation_mode: None,
             methylation_conversion_rate: None,
             methylation_failure_rate: None,
diff --git a/src/error_model/mod.rs b/src/error_model/mod.rs
index 7033dc1..a1cb01a 100644
--- a/src/error_model/mod.rs
+++ b/src/error_model/mod.rs
@@ -90,7 +90,7 @@ pub fn apply_errors(
 ///
 /// Given a DNA base (A, C, G, T), returns one of the other three bases
 /// chosen uniformly at random using a single RNG draw and a lookup table.
-fn random_different_base(base: u8, rng: &mut impl Rng) -> u8 {
+pub(crate) fn random_different_base(base: u8, rng: &mut impl Rng) -> u8 {
     // For each input base, the three possible substitute bases.
     const ALTS: [u8; 12] = [
         b'C', b'G', b'T', // A -> C, G, T
diff --git a/src/lib.rs b/src/lib.rs
index 7494fc7..23eb87d 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -5,6 +5,7 @@
 compile_error!("holodeck requires a 64-bit or wider platform");
 
 pub mod bed;
+pub mod clip;
 pub mod commands;
 pub mod error_model;
 pub mod fasta;
diff --git a/src/read.rs b/src/read.rs
index 79a4298..861046f 100644
--- a/src/read.rs
+++ b/src/read.rs
@@ -9,6 +9,7 @@ use noodles::sam::alignment::record::cigar::op::{Kind, Op};
 use noodles::sam::alignment::record_buf::Cigar;
 use rand::Rng;
 
+use crate::clip::TerminalClipConfig;
 use crate::error_model::{self, ErrorModel, ReadEnd};
 use crate::fragment::{
     Fragment, extract_read_bases, lowercase_fraction, reverse_complement, uppercase_in_place,
@@ -40,8 +41,8 @@ pub struct ReadPair {
     pub r1_truth: TruthAlignment,
     /// Truth alignment for R2 (present only for paired-end mode).
     pub r2_truth: Option<TruthAlignment>,
-    /// Truth CIGAR for R1, reflecting haplotype variants and adapter
-    /// soft-clipping.
+    /// Truth CIGAR for R1, reflecting haplotype variants plus adapter
+    /// read-through and terminal-artifact soft-clips.
     pub r1_cigar: Cigar,
     /// Truth CIGAR for R2. `None` for single-end reads.
     pub r2_cigar: Option<Cigar>,
@@ -61,7 +62,7 @@ struct MateOutput {
     /// Pre-conversion bases (5'→3' read orientation). `None` when methylation
     /// chemistry was not applied or capture was not requested.
     pre_conversion: Option<Vec<u8>>,
-    /// CIGAR including any adapter soft-clip.
+    /// CIGAR including any adapter and terminal-artifact soft-clips.
     cigar: Cigar,
     /// Truth alignment for this mate.
     truth: TruthAlignment,
@@ -116,6 +117,34 @@ fn apply_fragment_chemistry(
     (bases, conversion_failed)
 }
 
+/// Overwrite the soft-clipped terminal bases of a read with divergent
+/// sequence.
+///
+/// `bases` is in read 5'→3' orientation, so the 5' clip covers the first
+/// `clip_5p` bases and the 3' clip covers the `clip_3p` genomic bases ending
+/// at `genomic` (the adapter pad, if any, sits beyond `genomic` and is left
+/// untouched). Each clipped base is replaced with a base *different* from the
+/// one already there (reusing the error model's
+/// [`random_different_base`](error_model::random_different_base)), so the ends
+/// are guaranteed to diverge from the reference and a downstream aligner
+/// re-derives the soft-clip rather than forcing matches. Drawing uniformly
+/// from all four bases would, for a short clip, frequently regenerate the
+/// original base and leave the artifact invisible.
+fn corrupt_clipped_bases(
+    bases: &mut [u8],
+    genomic: usize,
+    clip_5p: usize,
+    clip_3p: usize,
+    rng: &mut impl Rng,
+) {
+    for b in &mut bases[..clip_5p] {
+        *b = error_model::random_different_base(*b, rng);
+    }
+    for b in &mut bases[genomic - clip_3p..genomic] {
+        *b = error_model::random_different_base(*b, rng);
+    }
+}
+
 /// Build one mate (R1 or R2) from already-extracted bases.
 ///
 /// `bases` is post-chemistry, post-uppercase, in read 5'→3' orientation
@@ -124,9 +153,11 @@ fn apply_fragment_chemistry(
 /// this function: errors mutate `bases` and advance `rng`, so rejection
 /// must happen first to avoid desynchronising the RNG stream.
 ///
-/// Pipeline (per-mate): apply_errors → CIGAR + truth.
+/// Pipeline (per-mate): apply_errors → terminal clip → CIGAR + truth.
 /// Chemistry runs once at fragment scale upstream — see
-/// [`apply_fragment_chemistry`].
+/// [`apply_fragment_chemistry`]. `clip_config`, when present and enabled,
+/// injects terminal soft-clip artifacts; a disabled or absent config draws no
+/// randomness and leaves the read end-to-end aligned.
 #[allow(clippy::too_many_arguments)]
 fn build_mate(
     fragment: &Fragment,
@@ -137,6 +168,7 @@ fn build_mate(
     pre_conversion: Option<Vec<u8>>,
     adapter_bases: usize,
     model: &impl ErrorModel,
+    clip_config: Option<&TerminalClipConfig>,
     rng: &mut impl Rng,
 ) -> MateOutput {
     let frag_len = fragment.bases.len();
@@ -153,15 +185,34 @@ fn build_mate(
 
     let (n_errors, qualities) = error_model::apply_errors(model, &mut bases, end, rng);
 
-    let cigar = cigar_from_ref_positions(positions, adapter_bases, is_negative_strand);
+    // Terminal soft-clip artifacts: sample per-end clip lengths over the
+    // aligned portion, corrupt those bases, and fold the clips into the CIGAR
+    // and alignment start so ground truth stays exact. A disabled/absent
+    // config yields (0, 0) without touching the RNG.
+    let (clip_5p, clip_3p) = clip_config.map_or((0, 0), |c| c.sample_clips(genomic, rng));
+    if clip_5p > 0 || clip_3p > 0 {
+        corrupt_clipped_bases(&mut bases, genomic, clip_5p, clip_3p, rng);
+    }
 
-    let ref_pos = if fragment.ref_positions.is_empty() {
-        0
-    } else if is_negative_strand {
-        fragment.ref_positions[right_start] + 1
+    // Read 5'→3' maps to ascending reference positions for a forward read but
+    // to descending for a reverse read (stored reverse-complemented in BAM),
+    // so a read-5' clip trims the front of `positions` on the forward strand
+    // and the back on the reverse strand.
+    let (front_trim, back_trim) =
+        if is_negative_strand { (clip_3p, clip_5p) } else { (clip_5p, clip_3p) };
+    let kept = &positions[front_trim..positions.len() - back_trim];
+
+    // Adapter read-through sits at the read's 3' end — the left of the record
+    // for a reverse read, the right for a forward read. Terminal clips add to
+    // the same record ends as their position trims.
+    let (lead_clip, trail_clip) = if is_negative_strand {
+        (front_trim + adapter_bases, back_trim)
     } else {
-        fragment.ref_positions[0] + 1
+        (front_trim, back_trim + adapter_bases)
     };
+    let cigar = cigar_from_ref_positions(kept, lead_clip, trail_clip);
+
+    let ref_pos = if kept.is_empty() { 0 } else { kept[0] + 1 };
 
     #[expect(clippy::cast_possible_truncation, reason = "fragment length fits in u32")]
     let fragment_length = frag_len as u32;
@@ -214,6 +265,9 @@ fn build_mate(
 /// * `capture_pre_conversion` — When `true` AND `methylation` is `Some`,
 ///   capture the pre-conversion bases of each mate for the `YS:Z` golden-BAM
 ///   tag. Has no effect when `methylation` is `None`.
+/// * `clip_config` — Optional terminal soft-clip artifact model. Each mate
+///   independently samples 5'/3' clips; a disabled or absent config leaves
+///   reads end-to-end aligned and draws no randomness.
 /// * `rng` — Random number generator.
 #[allow(clippy::too_many_arguments)] // Orchestrator for the read-pair pipeline
 pub fn generate_read_pair(
@@ -229,6 +283,7 @@ pub fn generate_read_pair(
     simple_names: bool,
     methylation: Option<&MethylationConfig>,
     capture_pre_conversion: bool,
+    clip_config: Option<&TerminalClipConfig>,
     rng: &mut impl Rng,
 ) -> Option<ReadPair> {
     let frag_len = fragment.bases.len();
@@ -297,6 +352,7 @@ pub fn generate_read_pair(
         r1_pre_conversion,
         adapter_bases,
         model,
+        clip_config,
         rng,
     );
 
@@ -336,6 +392,7 @@ pub fn generate_read_pair(
         r2_pre_conversion,
         adapter_bases,
         model,
+        clip_config,
         rng,
     );
 
@@ -365,8 +422,8 @@ pub fn generate_read_pair(
     })
 }
 
-/// Compute a CIGAR from a slice of ascending reference positions, plus
-/// optional adapter soft-clipping.
+/// Compute a CIGAR from a slice of ascending reference positions, with
+/// explicit leading and trailing soft-clip lengths.
 ///
 /// Positions must be in ascending order (forward strand). Consecutive
 /// positions incrementing by 1 produce M ops, same position produces I ops,
@@ -376,26 +433,32 @@ pub fn generate_read_pair(
 /// reference order from the leftmost aligned position, so even negative-strand
 /// reads use ascending positions.
 ///
-/// Adapter bases are appended as a soft-clip (S) operation. For
-/// negative-strand reads (`negative_strand = true`) the adapter is sequenced
-/// at the 3' end of the read but sits at the left (5') end of the stored BAM
-/// record after reverse-complementing, so the S op is placed at the *start*
-/// of the CIGAR. For forward-strand reads the S op is placed at the *end*.
+/// `lead_clip` and `trail_clip` are soft-clip lengths in BAM-record
+/// orientation (left and right of the alignment, respectively). The caller is
+/// responsible for folding both adapter read-through and terminal-artifact
+/// clips into these counts according to strand — for a forward-strand read the
+/// 3' adapter sits at the right (trailing), while for a reverse-strand read it
+/// sits at the left (leading) after reverse-complementing.
+///
+/// When `positions` is empty the read has no aligned bases, so the two clip
+/// counts collapse into a single soft-clip op (a CIGAR cannot carry two `S`
+/// ops with nothing between them).
 #[must_use]
-pub fn cigar_from_ref_positions(
-    positions: &[u32],
-    adapter_bases: usize,
-    negative_strand: bool,
-) -> Cigar {
+pub fn cigar_from_ref_positions(positions: &[u32], lead_clip: usize, trail_clip: usize) -> Cigar {
     let mut ops: Vec<Op> = Vec::new();
 
     if positions.is_empty() {
-        if adapter_bases > 0 {
-            ops.push(Op::new(Kind::SoftClip, adapter_bases));
+        let total = lead_clip + trail_clip;
+        if total > 0 {
+            ops.push(Op::new(Kind::SoftClip, total));
         }
         return Cigar::from(ops);
     }
 
+    if lead_clip > 0 {
+        ops.push(Op::new(Kind::SoftClip, lead_clip));
+    }
+
     let mut match_run: usize = 1; // First base is always a match.
     let mut ins_run: usize = 0;
 
@@ -442,15 +505,8 @@ pub fn cigar_from_ref_positions(
         ops.push(Op::new(Kind::Match, match_run));
     }
 
-    // Adapter bases as soft-clip. Negative-strand reads are stored
-    // reverse-complemented in BAM, so the adapter (at the 3' end in read
-    // order) moves to the left (5') end of the record — a leading S op.
-    if adapter_bases > 0 {
-        if negative_strand {
-            ops.insert(0, Op::new(Kind::SoftClip, adapter_bases));
-        } else {
-            ops.push(Op::new(Kind::SoftClip, adapter_bases));
-        }
+    if trail_clip > 0 {
+        ops.push(Op::new(Kind::SoftClip, trail_clip));
     }
 
     Cigar::from(ops)
@@ -503,69 +559,77 @@ mod tests {
 
     #[test]
     fn test_cigar_all_match() {
-        let cigar = cigar_from_ref_positions(&[0, 1, 2, 3, 4], 0, false);
+        let cigar = cigar_from_ref_positions(&[0, 1, 2, 3, 4], 0, 0);
         assert_eq!(cigar_to_string(&cigar), "5M");
     }
 
     #[test]
     fn test_cigar_with_insertion() {
         // Positions 0,1,2,2,2,3,4: two inserted bases at ref pos 2.
-        let cigar = cigar_from_ref_positions(&[0, 1, 2, 2, 2, 3, 4], 0, false);
+        let cigar = cigar_from_ref_positions(&[0, 1, 2, 2, 2, 3, 4], 0, 0);
         assert_eq!(cigar_to_string(&cigar), "3M2I2M");
     }
 
     #[test]
     fn test_cigar_with_deletion() {
         // Gap from 2 to 5: 2 deleted ref bases.
-        let cigar = cigar_from_ref_positions(&[0, 1, 2, 5, 6], 0, false);
+        let cigar = cigar_from_ref_positions(&[0, 1, 2, 5, 6], 0, 0);
         assert_eq!(cigar_to_string(&cigar), "3M2D2M");
     }
 
     #[test]
-    fn test_cigar_with_adapter_softclip_forward() {
-        // Forward-strand: adapter is a trailing soft-clip.
-        let cigar = cigar_from_ref_positions(&[0, 1, 2], 2, false);
+    fn test_cigar_with_trailing_softclip() {
+        // Forward-strand adapter (or 3' clip): trailing soft-clip.
+        let cigar = cigar_from_ref_positions(&[0, 1, 2], 0, 2);
         assert_eq!(cigar_to_string(&cigar), "3M2S");
     }
 
     #[test]
-    fn test_cigar_with_adapter_softclip_negative_strand() {
-        // Negative-strand: adapter moves to a leading soft-clip after RC.
-        let cigar = cigar_from_ref_positions(&[0, 1, 2], 2, true);
+    fn test_cigar_with_leading_softclip() {
+        // Negative-strand adapter (or read-3' clip after RC): leading soft-clip.
+        let cigar = cigar_from_ref_positions(&[0, 1, 2], 2, 0);
         assert_eq!(cigar_to_string(&cigar), "2S3M");
     }
 
     #[test]
-    fn test_cigar_all_adapter() {
-        // All-adapter read: placement doesn't matter, but it still round-trips.
-        let cigar = cigar_from_ref_positions(&[], 5, false);
+    fn test_cigar_with_both_leading_and_trailing_softclip() {
+        // A terminal clip on one end plus adapter on the other yields a clip
+        // at both ends of the record.
+        let cigar = cigar_from_ref_positions(&[5, 6, 7], 2, 3);
+        assert_eq!(cigar_to_string(&cigar), "2S3M3S");
+    }
+
+    #[test]
+    fn test_cigar_all_softclip() {
+        // No aligned bases: leading and trailing clips collapse to one S op.
+        let cigar = cigar_from_ref_positions(&[], 2, 3);
         assert_eq!(cigar_to_string(&cigar), "5S");
     }
 
     #[test]
     fn test_cigar_with_insertion_and_deletion() {
         // Insertion at pos 2 (two extra bases), then deletion of 2 ref bases.
-        let cigar = cigar_from_ref_positions(&[0, 1, 2, 2, 5, 6], 0, false);
+        let cigar = cigar_from_ref_positions(&[0, 1, 2, 2, 5, 6], 0, 0);
         assert_eq!(cigar_to_string(&cigar), "3M1I2D2M");
     }
 
     #[test]
-    fn test_cigar_with_adapter_and_deletion() {
-        let cigar = cigar_from_ref_positions(&[0, 1, 4, 5], 3, false);
+    fn test_cigar_with_trailing_clip_and_deletion() {
+        let cigar = cigar_from_ref_positions(&[0, 1, 4, 5], 0, 3);
         assert_eq!(cigar_to_string(&cigar), "2M2D2M3S");
     }
 
     #[test]
     fn test_cigar_single_base() {
-        let cigar = cigar_from_ref_positions(&[42], 0, false);
+        let cigar = cigar_from_ref_positions(&[42], 0, 0);
         assert_eq!(cigar_to_string(&cigar), "1M");
     }
 
     #[test]
     fn test_cigar_high_positions() {
         // Negative-strand R2: positions start from a high offset (ascending),
-        // no adapter — CIGAR is identical to forward strand.
-        let cigar = cigar_from_ref_positions(&[100, 101, 102, 103, 104], 0, true);
+        // no clips — CIGAR is identical to forward strand.
+        let cigar = cigar_from_ref_positions(&[100, 101, 102, 103, 104], 0, 0);
         assert_eq!(cigar_to_string(&cigar), "5M");
     }
 
@@ -579,7 +643,7 @@ mod tests {
 
         let pair = generate_read_pair(
             &fragment, "chr1", 1, 10, true, b"ADAPTER", b"ADAPTER", 1.0, &model, false, None,
-            false, &mut rng,
+            false, None, &mut rng,
         )
         .expect("no ambiguous bases — should not reject");
 
@@ -599,7 +663,7 @@ mod tests {
 
         let pair = generate_read_pair(
             &fragment, "chr1", 5, 10, false, b"ADAPTER", b"ADAPTER", 1.0, &model, false, None,
-            false, &mut rng,
+            false, None, &mut rng,
         )
         .unwrap();
 
@@ -618,7 +682,7 @@ mod tests {
 
         let pair = generate_read_pair(
             &fragment, "chr1", 1, 5, true, b"TTTTT", b"GGGGG", 1.0, &model, false, None, false,
-            &mut rng,
+            None, &mut rng,
         )
         .unwrap();
 
@@ -626,6 +690,177 @@ mod tests {
         assert_eq!(cigar_to_string(pair.r2_cigar.as_ref().unwrap()), "3S2M");
     }
 
+    #[test]
+    fn test_terminal_clip_5p_forward_shifts_pos_and_leading_clips() {
+        use crate::clip::TerminalClipConfig;
+        // length_mean = 1 makes the clip length deterministic (always 1 bp),
+        // rate 1.0 makes the clip certain — so the test asserts exact CIGAR/POS.
+        let clip = TerminalClipConfig { rate_5p: 1.0, rate_3p: 0.0, length_mean: 1, length_max: 5 };
+        let fragment = test_fragment(b"ACGTACGTACGTACGTACGT", 100); // forward, 20 bp.
+        let model = IlluminaErrorModel::new(20, 0.0, 0.0);
+        let mut rng = SmallRng::seed_from_u64(7);
+
+        let pair = generate_read_pair(
+            &fragment,
+            "chr1",
+            1,
+            20,
+            true,
+            b"ADAPTER",
+            b"ADAPTER",
+            1.0,
+            &model,
+            false,
+            None,
+            false,
+            Some(&clip),
+            &mut rng,
+        )
+        .unwrap();
+
+        // R1 is forward strand: read 5' is the left of the record, so the clip
+        // is leading and the alignment start advances past the clipped base.
+        assert_eq!(cigar_to_string(&pair.r1_cigar), "1S19M");
+        assert_eq!(pair.r1_truth.position, 102, "5' clip on a forward read must advance POS by 1");
+        // R2 is reverse strand: read 5' is the right of the record, so its own
+        // 5' clip is trailing and POS is unchanged.
+        assert_eq!(cigar_to_string(pair.r2_cigar.as_ref().unwrap()), "19M1S");
+        assert_eq!(pair.r2_truth.as_ref().unwrap().position, 101);
+    }
+
+    #[test]
+    fn test_terminal_clip_bases_diverge_from_reference() {
+        use crate::clip::TerminalClipConfig;
+        // With the error model disabled, the only mutation to R1's genomic
+        // bases is the clip. Every clipped 5' base must differ from the
+        // corresponding original fragment base, so the clip is guaranteed
+        // visible to a downstream aligner (not a coincidental match).
+        let clip = TerminalClipConfig { rate_5p: 1.0, rate_3p: 0.0, length_mean: 6, length_max: 8 };
+        let original = b"ACGTACGTACGTACGTACGT";
+        let fragment = test_fragment(original, 100);
+        let model = IlluminaErrorModel::new(20, 0.0, 0.0);
+        let mut rng = SmallRng::seed_from_u64(3);
+
+        let pair = generate_read_pair(
+            &fragment,
+            "chr1",
+            1,
+            20,
+            false,
+            b"ADAPTER",
+            b"ADAPTER",
+            1.0,
+            &model,
+            false,
+            None,
+            false,
+            Some(&clip),
+            &mut rng,
+        )
+        .unwrap();
+
+        // R1 is forward and only the 5' end clips, so the CIGAR is `<k>S<rest>M`.
+        // Recover the clip length from the leading soft-clip op rather than
+        // assuming a fixed sample.
+        let ops = pair.r1_cigar.as_ref();
+        assert_eq!(ops[0].kind(), Kind::SoftClip, "expected a leading 5' soft-clip");
+        let clip_len = ops[0].len();
+        assert!(clip_len >= 1);
+        for (i, (&got, &orig)) in
+            pair.read1.bases[..clip_len].iter().zip(&original[..clip_len]).enumerate()
+        {
+            assert_ne!(got, orig, "clipped 5' base {i} must differ from the reference base");
+        }
+        // The aligned interior is untouched.
+        assert_eq!(&pair.read1.bases[clip_len..], &original[clip_len..]);
+    }
+
+    #[test]
+    fn test_terminal_clip_3p_forward_is_trailing() {
+        use crate::clip::TerminalClipConfig;
+        let clip = TerminalClipConfig { rate_5p: 0.0, rate_3p: 1.0, length_mean: 1, length_max: 5 };
+        let fragment = test_fragment(b"ACGTACGTACGTACGTACGT", 100);
+        let model = IlluminaErrorModel::new(20, 0.0, 0.0);
+        let mut rng = SmallRng::seed_from_u64(7);
+
+        let pair = generate_read_pair(
+            &fragment,
+            "chr1",
+            1,
+            20,
+            true,
+            b"ADAPTER",
+            b"ADAPTER",
+            1.0,
+            &model,
+            false,
+            None,
+            false,
+            Some(&clip),
+            &mut rng,
+        )
+        .unwrap();
+
+        // Forward read, 3' clip → trailing, POS unchanged.
+        assert_eq!(cigar_to_string(&pair.r1_cigar), "19M1S");
+        assert_eq!(pair.r1_truth.position, 101);
+    }
+
+    #[test]
+    fn test_terminal_clip_disabled_is_identical_to_none() {
+        use crate::clip::TerminalClipConfig;
+        // A disabled config (both rates 0) must produce byte-identical output
+        // to passing no config at all — same bases, qualities, CIGAR, POS.
+        let disabled =
+            TerminalClipConfig { rate_5p: 0.0, rate_3p: 0.0, length_mean: 8, length_max: 20 };
+        let fragment = test_fragment(b"ACGTACGTACGTACGTACGT", 100);
+        let model = IlluminaErrorModel::new(20, 0.01, 0.05);
+
+        let mut rng_none = SmallRng::seed_from_u64(123);
+        let none = generate_read_pair(
+            &fragment,
+            "chr1",
+            1,
+            20,
+            true,
+            b"ADAPTER",
+            b"ADAPTER",
+            1.0,
+            &model,
+            false,
+            None,
+            false,
+            None,
+            &mut rng_none,
+        )
+        .unwrap();
+
+        let mut rng_dis = SmallRng::seed_from_u64(123);
+        let dis = generate_read_pair(
+            &fragment,
+            "chr1",
+            1,
+            20,
+            true,
+            b"ADAPTER",
+            b"ADAPTER",
+            1.0,
+            &model,
+            false,
+            None,
+            false,
+            Some(&disabled),
+            &mut rng_dis,
+        )
+        .unwrap();
+
+        assert_eq!(none.read1.bases, dis.read1.bases);
+        assert_eq!(none.read1.qualities, dis.read1.qualities);
+        assert_eq!(cigar_to_string(&none.r1_cigar), cigar_to_string(&dis.r1_cigar));
+        assert_eq!(none.r1_truth.position, dis.r1_truth.position);
+        assert_eq!(none.read2.as_ref().unwrap().bases, dis.read2.as_ref().unwrap().bases);
+    }
+
     #[test]
     fn test_simple_name_mode() {
         let fragment = test_fragment(b"ACGT", 0);
@@ -633,7 +868,8 @@ mod tests {
         let mut rng = SmallRng::seed_from_u64(42);
 
         let pair = generate_read_pair(
-            &fragment, "chr1", 42, 4, true, b"A", b"A", 1.0, &model, true, None, false, &mut rng,
+            &fragment, "chr1", 42, 4, true, b"A", b"A", 1.0, &model, true, None, false, None,
+            &mut rng,
         )
         .unwrap();
 
@@ -647,7 +883,8 @@ mod tests {
         let mut rng = SmallRng::seed_from_u64(42);
 
         let pair = generate_read_pair(
-            &fragment, "chr1", 1, 10, true, b"A", b"A", 1.0, &model, false, None, false, &mut rng,
+            &fragment, "chr1", 1, 10, true, b"A", b"A", 1.0, &model, false, None, false, None,
+            &mut rng,
         )
         .unwrap();
 
@@ -667,7 +904,7 @@ mod tests {
 
         let pair = generate_read_pair(
             &fragment, "chr1", 1, 10, true, b"ADAPTER", b"ADAPTER", 0.5, &model, false, None,
-            false, &mut rng,
+            false, None, &mut rng,
         );
 
         assert!(pair.is_none(), "all-lowercase fragment should be rejected at threshold 0.5");
@@ -683,7 +920,7 @@ mod tests {
 
         let pair = generate_read_pair(
             &fragment, "chr1", 1, 10, true, b"ADAPTER", b"ADAPTER", 0.5, &model, false, None,
-            false, &mut rng,
+            false, None, &mut rng,
         )
         .expect("0.3 < 0.5 — should accept");
 
@@ -730,6 +967,7 @@ mod tests {
             false,
             Some(&mc),
             true,
+            None,
             &mut rng,
         )
         .unwrap();
@@ -799,6 +1037,7 @@ mod tests {
             false,
             Some(&mc),
             true,
+            None,
             &mut rng,
         )
         .unwrap();
@@ -871,6 +1110,7 @@ mod tests {
             false,
             Some(&mc),
             false,
+            None,
             &mut rng,
         )
         .unwrap();
@@ -925,6 +1165,7 @@ mod tests {
             false,
             Some(&mc),
             false,
+            None,
             &mut rng,
         )
         .unwrap();
@@ -969,6 +1210,7 @@ mod tests {
             false,
             None,
             false,
+            None,
             &mut rng_a,
         );
         assert!(rejected.is_none(), "expected rejection when R2 is all-lowercase");
@@ -985,6 +1227,7 @@ mod tests {
             false,
             None,
             false,
+            None,
             &mut rng_a,
         )
         .unwrap();
@@ -1004,6 +1247,7 @@ mod tests {
             false,
             None,
             false,
+            None,
             &mut rng_b,
         )
         .unwrap();
diff --git a/tests/test_simulate.rs b/tests/test_simulate.rs
index 27e1ff9..4ccac5f 100644
--- a/tests/test_simulate.rs
+++ b/tests/test_simulate.rs
@@ -12,7 +12,7 @@ use std::path::PathBuf;
 
 use helpers::{
     PileupColumn, TestEnv, VcfVariant, count_fastq_records, non_repetitive_seq, pileup_bases,
-    read_gzipped, run_simulate, simple_env,
+    read_bam_records, read_gzipped, run_simulate, simple_env,
 };
 use std::collections::HashMap;
 
@@ -3506,3 +3506,167 @@ fn test_simulate_reproducible_ambiguity_resolution() {
     let r1_b = read_gzipped(&PathBuf::from(format!("{}.r1.fastq.gz", out_b.display())));
     assert_eq!(r1_a, r1_b, "seeded runs should be byte-identical");
 }
+
+/// Sum of soft-clip lengths across a record's CIGAR.
+fn total_softclip_bp(record: &RecordBuf) -> usize {
+    record
+        .cigar()
+        .as_ref()
+        .iter()
+        .filter(|op| op.kind() == CigarKind::SoftClip)
+        .map(|op| op.len())
+        .sum()
+}
+
+#[test]
+fn test_terminal_clip_appears_in_golden_bam() {
+    // With a high 5' clip rate and no adapter read-through (fragments well
+    // longer than the read), the golden BAM should carry soft-clips that did
+    // not exist before the model. The default run (no clip flags) must have
+    // essentially none, so the difference is attributable to the model.
+    let seq = non_repetitive_seq(5000);
+    let env = TestEnv::new(&[("chr1", &seq)]);
+
+    let run = |out: &std::path::Path, clip_args: &[&str]| {
+        let mut args = vec![
+            "simulate",
+            "-r",
+            env.fasta_path.to_str().unwrap(),
+            "-o",
+            out.to_str().unwrap(),
+            "--coverage",
+            "30",
+            "--read-length",
+            "100",
+            "--fragment-mean",
+            "400",
+            "--fragment-stddev",
+            "20",
+            "--golden-bam",
+            "--seed",
+            "7",
+        ];
+        args.extend_from_slice(clip_args);
+        let (ok, _, stderr) = run_simulate(&args);
+        assert!(ok, "simulate failed: {stderr}");
+        read_bam_records(&PathBuf::from(format!("{}.golden.bam", out.display())))
+    };
+
+    let baseline = run(&env.dir.path().join("baseline"), &[]);
+    let clipped =
+        run(&env.dir.path().join("clipped"), &["--clip-5p-rate", "0.5", "--clip-3p-rate", "0.0"]);
+
+    let baseline_clipped_reads = baseline.iter().filter(|r| total_softclip_bp(r) > 0).count();
+    let clipped_clipped_reads = clipped.iter().filter(|r| total_softclip_bp(r) > 0).count();
+
+    // Fragments (mean 400) are far longer than the 100bp reads, so adapter
+    // read-through is essentially impossible: the baseline should be clip-free.
+    assert_eq!(
+        baseline_clipped_reads, 0,
+        "baseline (no clip model) unexpectedly has soft-clipped records"
+    );
+
+    // ~50% 5' clip rate over both mates → a large fraction of records clipped.
+    let frac = clipped_clipped_reads as f64 / clipped.len() as f64;
+    assert!(
+        (0.40..=0.60).contains(&frac),
+        "expected ~50% of records soft-clipped, got {frac} ({clipped_clipped_reads}/{})",
+        clipped.len()
+    );
+
+    // Every soft-clip must respect the configured maximum length per end.
+    for r in &clipped {
+        assert!(
+            total_softclip_bp(r) <= 20,
+            "soft-clip exceeds --clip-length-max: {}",
+            helpers_cigar(r)
+        );
+    }
+}
+
+#[test]
+fn test_terminal_clip_preserves_truth_alignment_start() {
+    // A clip injected on the genomic read must keep the recorded alignment
+    // exact: the un-clipped reference span still aligns end-to-end (only M/D
+    // ops between the soft-clips), so no mismatch is introduced into the
+    // aligned block. We verify by checking that every clipped record's aligned
+    // portion is a clean M run against the reference at its POS.
+    let seq = non_repetitive_seq(5000);
+    let env = TestEnv::new(&[("chr1", &seq)]);
+    let out = env.dir.path().join("clip_truth");
+
+    let (ok, _, stderr) = run_simulate(&[
+        "simulate",
+        "-r",
+        env.fasta_path.to_str().unwrap(),
+        "-o",
+        out.to_str().unwrap(),
+        "--coverage",
+        "20",
+        "--read-length",
+        "100",
+        "--fragment-mean",
+        "400",
+        "--fragment-stddev",
+        "20",
+        "--golden-bam",
+        "--clip-5p-rate",
+        "0.5",
+        "--clip-3p-rate",
+        "0.3",
+        "--seed",
+        "11",
+    ]);
+    assert!(ok, "simulate failed: {stderr}");
+
+    let records = read_bam_records(&PathBuf::from(format!("{}.golden.bam", out.display())));
+    let mut saw_lead = false;
+    let mut saw_trail = false;
+    for r in &records {
+        let ops: Vec<_> = r.cigar().as_ref().iter().collect();
+        // The aligned interior (between any leading/trailing soft-clips) must
+        // consist only of M ops here (no variants in this reference), so the
+        // soft-clip is the only non-M op and ground truth is internally
+        // consistent.
+        let n = ops.len();
+        for (i, op) in ops.iter().enumerate() {
+            if op.kind() == CigarKind::SoftClip {
+                assert!(i == 0 || i == n - 1, "soft-clip must be terminal: {}", helpers_cigar(r));
+                if i == 0 {
+                    saw_lead = true;
+                } else {
+                    saw_trail = true;
+                }
+            } else {
+                assert_eq!(
+                    op.kind(),
+                    CigarKind::Match,
+                    "aligned interior should be all-M in a variant-free reference: {}",
+                    helpers_cigar(r)
+                );
+            }
+        }
+    }
+    assert!(saw_lead, "expected at least one leading soft-clip");
+    assert!(saw_trail, "expected at least one trailing soft-clip");
+}
+
+/// Format a record's CIGAR for assertion messages.
+fn helpers_cigar(record: &RecordBuf) -> String {
+    use std::fmt::Write;
+    record.cigar().as_ref().iter().fold(String::new(), |mut s, op| {
+        let c = match op.kind() {
+            CigarKind::Match => 'M',
+            CigarKind::Insertion => 'I',
+            CigarKind::Deletion => 'D',
+            CigarKind::SoftClip => 'S',
+            CigarKind::HardClip => 'H',
+            CigarKind::Skip => 'N',
+            CigarKind::Pad => 'P',
+            CigarKind::SequenceMatch => '=',
+            CigarKind::SequenceMismatch => 'X',
+        };
+        let _ = write!(s, "{}{c}", op.len());
+        s
+    })
+}