diff --git a/USING.md b/USING.md
index 0dae801426..6e2dbdf677 100644
--- a/USING.md
+++ b/USING.md
@@ -123,6 +123,7 @@ module org.foo.baz {
 ## JTS System Properties
 
 * `-Djts.overlay=ng` enables the use of OverlayNG in `Geometry` overlay methods. (*Note: in a future release this will become the default behaviour*) 
+* `-Djts.relate=ng` enables the use of RelateNG in `Geometry` topological predicate methods. (*Note: in a future release this will become the default behaviour*) 
 
 ## JTS Tools
 
diff --git a/modules/core/src/main/java/org/locationtech/jts/geom/Geometry.java b/modules/core/src/main/java/org/locationtech/jts/geom/Geometry.java
index ba0d353fb2..104a885c8e 100644
--- a/modules/core/src/main/java/org/locationtech/jts/geom/Geometry.java
+++ b/modules/core/src/main/java/org/locationtech/jts/geom/Geometry.java
@@ -713,10 +713,7 @@ public boolean disjoint(Geometry g) {
    *      Returns <code>false</code> if both <code>Geometry</code>s are points
    */
   public boolean touches(Geometry g) {
-    // short-circuit test
-    if (! getEnvelopeInternal().intersects(g.getEnvelopeInternal()))
-      return false;
-    return relate(g).isTouches(getDimension(), g.getDimension());
+    return GeometryRelate.touches(this, g);
   }
 
   /**
@@ -771,18 +768,8 @@ public boolean intersects(Geometry g) {
     if (g.isRectangle()) {
       return RectangleIntersects.intersects((Polygon) g, this);
     }
-    if (isGeometryCollection() || g.isGeometryCollection()) {
-      for (int i = 0 ; i < getNumGeometries() ; i++) {
-        for (int j = 0 ; j < g.getNumGeometries() ; j++) {
-          if (getGeometryN(i).intersects(g.getGeometryN(j))) {
-            return true;
-          }
-        }
-      }
-      return false;
-    }
-    // general case
-    return relate(g).isIntersects();
+
+    return GeometryRelate.intersects(this, g);
   }
 
   /**
@@ -845,7 +832,7 @@ public boolean crosses(Geometry g) {
    * @see Geometry#coveredBy
    */
   public boolean within(Geometry g) {
-    return g.contains(this);
+    return GeometryRelate.within(this, g);
   }
 
   /**
@@ -876,25 +863,13 @@ public boolean within(Geometry g) {
    * @see Geometry#covers
    */
   public boolean contains(Geometry g) {
-    // optimization - lower dimension cannot contain areas
-    if (g.getDimension() == 2 && getDimension() < 2) {
-      return false;
-    }
-    // optimization - P cannot contain a non-zero-length L
-    // Note that a point can contain a zero-length lineal geometry,
-    // since the line has no boundary due to Mod-2 Boundary Rule
-    if (g.getDimension() == 1 && getDimension() < 1 && g.getLength() > 0.0) {
-      return false;
-    }
-    // optimization - envelope test
-    if (! getEnvelopeInternal().contains(g.getEnvelopeInternal()))
-      return false;
+
     // optimization for rectangle arguments
     if (isRectangle()) {
       return RectangleContains.contains((Polygon) this, g);
     }
     // general case
-    return relate(g).isContains();
+    return GeometryRelate.contains(this, g);
   }
 
   /**
@@ -919,10 +894,7 @@ public boolean contains(Geometry g) {
    *@return        <code>true</code> if the two <code>Geometry</code>s overlap.
    */
   public boolean overlaps(Geometry g) {
-    // short-circuit test
-    if (! getEnvelopeInternal().intersects(g.getEnvelopeInternal()))
-      return false;
-    return relate(g).isOverlaps(getDimension(), g.getDimension());
+    return GeometryRelate.overlaps(this, g);
   }
 
   /**
@@ -960,24 +932,7 @@ public boolean overlaps(Geometry g) {
    * @see Geometry#coveredBy
    */
   public boolean covers(Geometry g) {
-    // optimization - lower dimension cannot cover areas
-    if (g.getDimension() == 2 && getDimension() < 2) {
-      return false;
-    }
-    // optimization - P cannot cover a non-zero-length L
-    // Note that a point can cover a zero-length lineal geometry
-    if (g.getDimension() == 1 && getDimension() < 1 && g.getLength() > 0.0) {
-      return false;
-    }
-    // optimization - envelope test
-    if (! getEnvelopeInternal().covers(g.getEnvelopeInternal()))
-      return false;
-    // optimization for rectangle arguments
-    if (isRectangle()) {
-    	// since we have already tested that the test envelope is covered
-      return true;
-    }
-    return relate(g).isCovers();
+    return GeometryRelate.covers(this, g);
   }
 
   /**
@@ -1010,7 +965,7 @@ public boolean covers(Geometry g) {
    * @see Geometry#covers
    */
   public boolean coveredBy(Geometry g) {
-    return g.covers(this);
+    return GeometryRelate.coveredBy(this, g);
   }
 
   /**
@@ -1037,7 +992,7 @@ public boolean coveredBy(Geometry g) {
    * @see IntersectionMatrix
    */
   public boolean relate(Geometry g, String intersectionPattern) {
-    return relate(g).matches(intersectionPattern);
+    return GeometryRelate.relate(this, g, intersectionPattern);
   }
 
   /**
@@ -1048,9 +1003,7 @@ public boolean relate(Geometry g, String intersectionPattern) {
    *      boundaries and exteriors of the two <code>Geometry</code>s
    */
   public IntersectionMatrix relate(Geometry g) {
-    checkNotGeometryCollection(this);
-    checkNotGeometryCollection(g);
-    return RelateOp.relate(this, g);
+    return GeometryRelate.relate(this, g);
   }
 
   /**
@@ -1101,10 +1054,7 @@ public boolean equals(Geometry g) {
    */
   public boolean equalsTopo(Geometry g)
   {
-    // short-circuit test
-    if (! getEnvelopeInternal().equals(g.getEnvelopeInternal()))
-      return false;
-    return relate(g).isEquals(getDimension(), g.getDimension());
+    return GeometryRelate.equalsTopo(this, g);
   }
 
   /**
diff --git a/modules/core/src/main/java/org/locationtech/jts/geom/GeometryRelate.java b/modules/core/src/main/java/org/locationtech/jts/geom/GeometryRelate.java
new file mode 100644
index 0000000000..615e7d1ab8
--- /dev/null
+++ b/modules/core/src/main/java/org/locationtech/jts/geom/GeometryRelate.java
@@ -0,0 +1,217 @@
+/*
+ * Copyright (c) 2020 Martin Davis.
+ *
+ * All rights reserved. This program and the accompanying materials
+ * are made available under the terms of the Eclipse Public License 2.0
+ * and Eclipse Distribution License v. 1.0 which accompanies this distribution.
+ * The Eclipse Public License is available at http://www.eclipse.org/legal/epl-v20.html
+ * and the Eclipse Distribution License is available at
+ *
+ * http://www.eclipse.org/org/documents/edl-v10.php.
+ */
+package org.locationtech.jts.geom;
+
+import org.locationtech.jts.operation.relate.RelateOp;
+import org.locationtech.jts.operation.relateng.RelateNG;
+import org.locationtech.jts.operation.relateng.RelatePredicate;
+
+/**
+ * Internal class which encapsulates the runtime switch to use RelateNG.
+ * <p>
+ * This class allows the {@link Geometry} predicate methods to be 
+ * switched between the original {@link RelateOp} algorithm 
+ * and the modern {@link RelateNG} codebase
+ * via a system property <code>jts.relate</code>.
+ * <ul>
+ * <li><code>jts.relate=old</code> - (default) use original RelateOp algorithm
+ * <li><code>jts.relate=ng</code> - use RelateNG
+ * </ul>
+ * 
+ * @author mdavis
+ *
+ */
+class GeometryRelate 
+{
+  public static String RELATE_PROPERTY_NAME = "jts.relate";
+  
+  public static String RELATE_PROPERTY_VALUE_NG = "ng";
+  public static String RELATE_PROPERTY_VALUE_OLD = "old";
+  
+  /**
+   * Currently the old relate implementation is the default
+   */
+  public static boolean RELATE_NG_DEFAULT = false;
+
+  private static boolean isRelateNG = RELATE_NG_DEFAULT;
+
+  static {
+    setRelateImpl(System.getProperty(RELATE_PROPERTY_NAME));
+  }
+  
+  /**
+   * This function is provided primarily for unit testing.
+   * It is not recommended to use it dynamically, since 
+   * that may result in inconsistent overlay behaviour.
+   * 
+   * @param relateImplCode the code for the overlay method (may be null)
+   */
+  static void setRelateImpl(String relateImplCode) {
+    if (relateImplCode == null) 
+      return;
+    // set flag explicitly since current value may not be default
+    isRelateNG = RELATE_NG_DEFAULT;
+    
+    if (RELATE_PROPERTY_VALUE_NG.equalsIgnoreCase(relateImplCode) )
+      isRelateNG = true;
+  }
+  
+  static boolean intersects(Geometry a, Geometry b)
+  {
+    if (isRelateNG) {
+      return RelateNG.relate(a, b, RelatePredicate.intersects());
+    }
+    if (a.isGeometryCollection() || b.isGeometryCollection()) {
+      for (int i = 0 ; i < a.getNumGeometries() ; i++) {
+        for (int j = 0 ; j < b.getNumGeometries() ; j++) {
+          if (a.getGeometryN(i).intersects(b.getGeometryN(j))) {
+            return true;
+          }
+        }
+      }
+      return false;
+    }
+    return RelateOp.relate(a, b).isIntersects();
+  }
+
+  static boolean contains(Geometry a, Geometry b)
+  {
+    if (isRelateNG) {
+      return RelateNG.relate(a, b, RelatePredicate.contains());
+    }
+    // optimization - lower dimension cannot contain areas
+    if (b.getDimension() == 2 && a.getDimension() < 2) {
+      return false;
+    }
+    // optimization - P cannot contain a non-zero-length L
+    // Note that a point can contain a zero-length lineal geometry,
+    // since the line has no boundary due to Mod-2 Boundary Rule
+    if (b.getDimension() == 1 && a.getDimension() < 1 && b.getLength() > 0.0) {
+      return false;
+    }
+    // optimization - envelope test
+    if (! a.getEnvelopeInternal().contains(b.getEnvelopeInternal()))
+      return false;
+    return RelateOp.relate(a, b).isContains();
+  }
+
+  static boolean covers(Geometry a, Geometry b)
+  {
+    if (isRelateNG) {
+      return RelateNG.relate(a, b, RelatePredicate.covers());
+    }
+    // optimization - lower dimension cannot cover areas
+    if (b.getDimension() == 2 && a.getDimension() < 2) {
+      return false;
+    }
+    // optimization - P cannot cover a non-zero-length L
+    // Note that a point can cover a zero-length lineal geometry
+    if (b.getDimension() == 1 && a.getDimension() < 1 && b.getLength() > 0.0) {
+      return false;
+    }
+    // optimization - envelope test
+    if (! a.getEnvelopeInternal().covers(b.getEnvelopeInternal()))
+      return false;
+    // optimization for rectangle arguments
+    if (a.isRectangle()) {
+      // since we have already tested that the test envelope is covered
+      return true;
+    }
+    return RelateOp.relate(a, b).isCovers();  
+  }
+  
+  static boolean coveredBy(Geometry a, Geometry b)
+  {
+    if (isRelateNG) {
+      return RelateNG.relate(a, b, RelatePredicate.coveredBy());
+    }
+    return covers(b, a);  
+  }
+  
+  static boolean crosses(Geometry a, Geometry b)
+  {
+    if (isRelateNG) {
+      return RelateNG.relate(a, b, RelatePredicate.crosses());
+    }
+    // short-circuit test
+    if (! a.getEnvelopeInternal().intersects(b.getEnvelopeInternal()))
+      return false;
+    return RelateOp.relate(a, b).isCrosses(a.getDimension(), b.getDimension()); 
+  }
+  
+  static boolean disjoint(Geometry a, Geometry b)
+  {
+    if (isRelateNG) {
+      return RelateNG.relate(a, b, RelatePredicate.disjoint());
+    }
+    return ! intersects(a, b);
+  }
+  
+  static boolean equalsTopo(Geometry a, Geometry b)
+  {
+    if (isRelateNG) {
+      return RelateNG.relate(a, b, RelatePredicate.equalsTopo());
+    }
+    if (! a.getEnvelopeInternal().equals(b.getEnvelopeInternal()))
+      return false;
+    return RelateOp.relate(a, b).isEquals(a.getDimension(), b.getDimension());  
+  }
+  
+  static boolean overlaps(Geometry a, Geometry b)
+  {
+    if (isRelateNG) {
+      return RelateNG.relate(a, b, RelatePredicate.overlaps());
+    }
+    if (! a.getEnvelopeInternal().intersects(b.getEnvelopeInternal()))
+      return false;
+    return RelateOp.relate(a, b).isOverlaps(a.getDimension(), b.getDimension());
+  }
+  
+  static boolean touches(Geometry a, Geometry b)
+  {
+    if (isRelateNG) {
+      return RelateNG.relate(a, b, RelatePredicate.touches());
+    }
+    if (! a.getEnvelopeInternal().intersects(b.getEnvelopeInternal()))
+      return false;
+    return RelateOp.relate(a, b).isTouches(a.getDimension(), b.getDimension());
+  }
+
+  static boolean within(Geometry a, Geometry b)
+  {
+    if (isRelateNG) {
+      return RelateNG.relate(a, b, RelatePredicate.within());
+    }
+    return contains(b, a);
+  }
+
+  static IntersectionMatrix relate(Geometry a, Geometry b)
+  {
+    if (isRelateNG) {
+      return RelateNG.relate(a, b);
+    }
+    Geometry.checkNotGeometryCollection(a);
+    Geometry.checkNotGeometryCollection(b);
+    return RelateOp.relate(a, b);
+  }
+  
+  static boolean relate(Geometry a, Geometry b, String intersectionPattern)
+  {
+    if (isRelateNG) {
+      return RelateNG.relate(a, b, intersectionPattern);
+    }
+    Geometry.checkNotGeometryCollection(a);
+    Geometry.checkNotGeometryCollection(b);
+    return RelateOp.relate(a, b).matches(intersectionPattern);
+  }
+  
+}