diff --git a/R/interface.R b/R/interface.R
index 41e3f4235..595800aed 100644
--- a/R/interface.R
+++ b/R/interface.R
@@ -79,6 +79,15 @@ population <- function(name, time, N, parent = NULL, map = FALSE,
       check_split_time(time, parent)
   }
 
+  parent_remove <- attr(parent, "remove")
+  if (!is.null(parent) && parent_remove != -1) {
+    direction <- if (time > parent$time) "forward" else "backward"
+    if (direction == "forward" && parent_remove <= time)
+      stop("Parent population will be removed before the split event", call. = FALSE)
+    else if (direction == "backward" && parent_remove >= time)
+      stop("Parent population will be removed before the split event", call. = FALSE)
+  }
+
   if (!is.logical(map) && !inherits(map, "slendr_map"))
     stop("A simulation landscape must be an object of the class slendr_map", call. = FALSE)
 
diff --git a/tests/testthat/test-geneflow.R b/tests/testthat/test-geneflow.R
index c2fe9bf98..040f98f3e 100644
--- a/tests/testthat/test-geneflow.R
+++ b/tests/testthat/test-geneflow.R
@@ -81,7 +81,7 @@ test_that("populations must be already created for a gene flow to happen (forwar
 # that the issue has been solved
 # https://github.com/bodkan/slendr/issues/132
 test_that("gene_flow() behaves as expected in some concrete situations", {
-  anc <- population("ancestor", time = 9e6, remove=8e6, N=100)
+  anc <- population("ancestor", time = 9e6, remove=6e6, N=100)
   popA <- population("A", time = 8e6, parent=anc, N=100)
   popB <- population("B", time = 7e6, parent=anc, N=100)
 
diff --git a/tests/testthat/test-population.R b/tests/testthat/test-population.R
index baf445273..4bdbaa1e0 100644
--- a/tests/testthat/test-population.R
+++ b/tests/testthat/test-population.R
@@ -23,6 +23,43 @@ test_that("non-integer population split time is rounded", {
   expect_true(attr(pop, "history")[[1]]$N == 101)
 })
 
+test_that("parent cannot be scheduled for removal before a daughter splits (forward)", {
+  error_msg <- "Parent population will be removed"
+
+  parent <- population("parent", time = 1, N = 1, remove = 50)
+  expect_error(population("daughter", time = 100, N = 1, parent = parent), error_msg)
+  expect_error(population("daughter", time = 50, N = 1, parent = parent), error_msg)
+  expect_s3_class(population("daughter", time = 30, N = 1, parent = parent), "slendr_pop")
+
+  parent <- population("parent", time = 100, N = 1, remove = 150)
+  expect_error(population("daughter", time = 200, N = 1, parent = parent), error_msg)
+  expect_error(population("daughter", time = 150, N = 1, parent = parent), error_msg)
+  expect_s3_class(daughter <- population("daughter", time = 120, N = 1, parent = parent), "slendr_pop")
+
+  model <- compile_model(list(parent, daughter), simulation_length = 300, generation_time = 10)
+
+  # successful model definition in slendr is one thing, but let's make sure the
+  # simulation themselves really run
+  expect_s3_class(slim(model, sequence_length = 1, recombination_rate = 0), "slendr_ts")
+  expect_s3_class(msprime(model, sequence_length = 1, recombination_rate = 0), "slendr_ts")
+})
+
+test_that("parent cannot be scheduled for removal before a daughter splits (backward)", {
+  error_msg <- "Parent population will be removed"
+
+  parent <- population("parent", time = 1000, N = 1, remove = 500)
+  expect_error(population("daughter", time = 100, N = 1, parent = parent), error_msg)
+  expect_error(population("daughter", time = 500, N = 1, parent = parent), error_msg)
+  expect_s3_class(daughter <- population("daughter", time = 800, N = 1, parent = parent), "slendr_pop")
+
+  model <- compile_model(list(parent, daughter), generation_time = 10)
+
+  # successful model definition in slendr is one thing, but let's make sure the
+  # simulation themselves really run
+  expect_s3_class(slim(model, sequence_length = 1, recombination_rate = 0), "slendr_ts")
+  expect_s3_class(msprime(model, sequence_length = 1, recombination_rate = 0), "slendr_ts")
+})
+
 #
 # population resizes (step)
 #