#10827 Fault reactivation: extract temperatures using well log.

ApplicationLibCode/FileInterface/RifFaultReactivationModelExporter.cpp

@@ -850,7 +850,8 @@ std::shared_ptr<RimFaultReactivationDataAccess>
     if ( eCase == nullptr ) return nullptr;
 
     // extract data for each timestep
-    auto dataAccess = std::make_shared<RimFaultReactivationDataAccess>( eCase, rimModel.geoMechCase(), rimModel.selectedTimeStepIndexes() );
+    auto dataAccess =
+        std::make_shared<RimFaultReactivationDataAccess>( rimModel, eCase, rimModel.geoMechCase(), rimModel.selectedTimeStepIndexes() );
     dataAccess->extractModelData( *rimModel.model() );
     return dataAccess;
 }

ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccess.cpp

@@ -39,19 +39,23 @@
 #include "RimFaultReactivationDataAccessorTemperature.h"
 #include "RimFaultReactivationDataAccessorVoidRatio.h"
 #include "RimFaultReactivationEnums.h"
+#include "RimFaultReactivationModel.h"
 
 //--------------------------------------------------------------------------------------------------
 ///
 //--------------------------------------------------------------------------------------------------
-RimFaultReactivationDataAccess::RimFaultReactivationDataAccess( RimEclipseCase*            thecase,
+RimFaultReactivationDataAccess::RimFaultReactivationDataAccess( const RimFaultReactivationModel& model,
-                                                                RimGeoMechCase*            geoMechCase,
+                                                                RimEclipseCase*                  thecase,
-                                                                const std::vector<size_t>& timeSteps )
+                                                                RimGeoMechCase*                  geoMechCase,
+                                                                const std::vector<size_t>&       timeSteps )
     : m_timeSteps( timeSteps )
 {
     double porePressureGradient = 1.0;
-    m_accessors.push_back( std::make_shared<RimFaultReactivationDataAccessorPorePressure>( thecase, porePressureGradient ) );
+    double topTemperature       = model.seabedTemperature();
+    double seabedDepth          = -model.seaBedDepth();
+    m_accessors.push_back( std::make_shared<RimFaultReactivationDataAccessorPorePressure>( thecase, porePressureGradient, seabedDepth ) );
     m_accessors.push_back( std::make_shared<RimFaultReactivationDataAccessorVoidRatio>( thecase, 0.0001 ) );
-    m_accessors.push_back( std::make_shared<RimFaultReactivationDataAccessorTemperature>( thecase ) );
+    m_accessors.push_back( std::make_shared<RimFaultReactivationDataAccessorTemperature>( thecase, topTemperature, seabedDepth ) );
     if ( geoMechCase )
     {
         std::vector<RimFaultReactivation::Property> properties = { RimFaultReactivation::Property::YoungsModulus,
@@ -71,7 +75,8 @@ RimFaultReactivationDataAccess::RimFaultReactivationDataAccess( RimEclipseCase*
                                                                          RimFaultReactivation::Property::LateralStressComponentY };
         for ( auto property : stressProperties )
         {
-            m_accessors.push_back( std::make_shared<RimFaultReactivationDataAccessorStress>( geoMechCase, property, porePressureGradient ) );
+            m_accessors.push_back(
+                std::make_shared<RimFaultReactivationDataAccessorStress>( geoMechCase, property, porePressureGradient, seabedDepth ) );
         }
     }
 }

ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccess.h

@@ -32,6 +32,7 @@ class RimEclipseCase;
 class RimGeoMechCase;
 class RimFaultReactivationDataAccessor;
 class RigFaultReactivationModel;
+class RimFaultReactivationModel;
 
 //==================================================================================================
 ///
@@ -40,7 +41,10 @@ class RigFaultReactivationModel;
 class RimFaultReactivationDataAccess
 {
 public:
-    RimFaultReactivationDataAccess( RimEclipseCase* eclipseCase, RimGeoMechCase* geoMechCase, const std::vector<size_t>& timeSteps );
+    RimFaultReactivationDataAccess( const RimFaultReactivationModel& model,
+                                    RimEclipseCase*                  eclipseCase,
+                                    RimGeoMechCase*                  geoMechCase,
+                                    const std::vector<size_t>&       timeSteps );
     ~RimFaultReactivationDataAccess();
 
     void extractModelData( const RigFaultReactivationModel& model );

ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorPorePressure.cpp

@@ -42,9 +42,11 @@
 ///
 //--------------------------------------------------------------------------------------------------
 RimFaultReactivationDataAccessorPorePressure::RimFaultReactivationDataAccessorPorePressure( RimEclipseCase* eclipseCase,
-                                                                                            double          porePressureGradient )
+                                                                                            double          porePressureGradient,
+                                                                                            double          seabedDepth )
     : m_eclipseCase( eclipseCase )
     , m_defaultPorePressureGradient( porePressureGradient )
+    , m_seabedDepth( seabedDepth )
     , m_caseData( nullptr )
     , m_mainGrid( nullptr )
 {
@@ -75,9 +77,10 @@ void RimFaultReactivationDataAccessorPorePressure::updateResultAccessor()
     m_resultAccessor =
         RigResultAccessorFactory::createFromResultAddress( m_caseData, 0, RiaDefines::PorosityModelType::MATRIX_MODEL, m_timeStep, resVarAddress );
 
-    auto [wellPaths, extractors] = RimFaultReactivationDataAccessorWellLogExtraction::createEclipseWellPathExtractors( *m_model, *m_caseData );
+    auto [wellPaths, extractors] =
-    m_wellPaths                  = wellPaths;
+        RimFaultReactivationDataAccessorWellLogExtraction::createEclipseWellPathExtractors( *m_model, *m_caseData, m_seabedDepth );
-    m_extractors                 = extractors;
+    m_wellPaths  = wellPaths;
+    m_extractors = extractors;
 }
 
 //--------------------------------------------------------------------------------------------------

ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorPorePressure.h

@@ -37,7 +37,7 @@ class RigWellPath;
 class RimFaultReactivationDataAccessorPorePressure : public RimFaultReactivationDataAccessor
 {
 public:
-    RimFaultReactivationDataAccessorPorePressure( RimEclipseCase* eclipseCase, double porePressureGradient );
+    RimFaultReactivationDataAccessorPorePressure( RimEclipseCase* eclipseCase, double porePressureGradient, double seabedDepth );
     ~RimFaultReactivationDataAccessorPorePressure();
 
     bool isMatching( RimFaultReactivation::Property property ) const override;
@@ -58,6 +58,7 @@ private:
     RigEclipseCaseData*         m_caseData;
     const RigMainGrid*          m_mainGrid;
     double                      m_defaultPorePressureGradient;
+    double                      m_seabedDepth;
     cvf::ref<RigResultAccessor> m_resultAccessor;
 
     std::map<RimFaultReactivation::GridPart, cvf::ref<RigWellPath>>                m_wellPaths;

ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStress.cpp

@@ -48,10 +48,12 @@
 //--------------------------------------------------------------------------------------------------
 RimFaultReactivationDataAccessorStress::RimFaultReactivationDataAccessorStress( RimGeoMechCase*                geoMechCase,
                                                                                 RimFaultReactivation::Property property,
-                                                                                double                         gradient )
+                                                                                double                         gradient,
+                                                                                double                         seabedDepth )
     : m_geoMechCase( geoMechCase )
     , m_property( property )
     , m_gradient( gradient )
+    , m_seabedDepth( seabedDepth )
 {
     m_geoMechCaseData = geoMechCase->geoMechData();
 }
@@ -99,6 +101,7 @@ void RimFaultReactivationDataAccessorStress::updateResultAccessor()
         std::vector<cvf::Vec3d> wellPoints =
             RimFaultReactivationDataAccessorWellLogExtraction::generateWellPoints( faultTopPosition,
                                                                                    faultBottomPosition,
+                                                                                   m_seabedDepth,
                                                                                    faultNormal * distanceFromFault );
         m_faceAWellPath = new RigWellPath( wellPoints, RimFaultReactivationDataAccessorWellLogExtraction::generateMds( wellPoints ) );
         m_partIndexA    = geoMechPartCollection->getPartIndexFromPoint( wellPoints[1] );
@@ -109,6 +112,7 @@ void RimFaultReactivationDataAccessorStress::updateResultAccessor()
         std::vector<cvf::Vec3d> wellPoints =
             RimFaultReactivationDataAccessorWellLogExtraction::generateWellPoints( faultTopPosition,
                                                                                    faultBottomPosition,
+                                                                                   m_seabedDepth,
                                                                                    -faultNormal * distanceFromFault );
         m_faceBWellPath = new RigWellPath( wellPoints, RimFaultReactivationDataAccessorWellLogExtraction::generateMds( wellPoints ) );
         m_partIndexB    = geoMechPartCollection->getPartIndexFromPoint( wellPoints[1] );

ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStress.h

@@ -49,7 +49,10 @@ class RigGeoMechWellLogExtractor;
 class RimFaultReactivationDataAccessorStress : public RimFaultReactivationDataAccessor
 {
 public:
-    RimFaultReactivationDataAccessorStress( RimGeoMechCase* geoMechCase, RimFaultReactivation::Property property, double gradient );
+    RimFaultReactivationDataAccessorStress( RimGeoMechCase*                geoMechCase,
+                                            RimFaultReactivation::Property property,
+                                            double                         gradient,
+                                            double                         seabedDepth );
     ~RimFaultReactivationDataAccessorStress();
 
     bool isMatching( RimFaultReactivation::Property property ) const override;
@@ -76,6 +79,7 @@ private:
     RimGeoMechCase*                m_geoMechCase;
     RimFaultReactivation::Property m_property;
     double                         m_gradient;
+    double                         m_seabedDepth;
     RigGeoMechCaseData*            m_geoMechCaseData;
     RigFemScalarResultFrames*      m_s11Frames;
     RigFemScalarResultFrames*      m_s22Frames;

ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorTemperature.cpp

@@ -25,10 +25,13 @@
 #include "RigCaseCellResultsData.h"
 #include "RigEclipseCaseData.h"
 #include "RigEclipseResultAddress.h"
+#include "RigEclipseWellLogExtractor.h"
 #include "RigMainGrid.h"
 #include "RigResultAccessorFactory.h"
+#include "RigWellPath.h"
 
 #include "RimEclipseCase.h"
+#include "RimFaultReactivationDataAccessorWellLogExtraction.h"
 
 #include <cmath>
 #include <limits>
@@ -36,10 +39,14 @@
 //--------------------------------------------------------------------------------------------------
 ///
 //--------------------------------------------------------------------------------------------------
-RimFaultReactivationDataAccessorTemperature::RimFaultReactivationDataAccessorTemperature( RimEclipseCase* eclipseCase )
+RimFaultReactivationDataAccessorTemperature::RimFaultReactivationDataAccessorTemperature( RimEclipseCase* eclipseCase,
+                                                                                          double          seabedTemperature,
+                                                                                          double          seabedDepth )
     : m_eclipseCase( eclipseCase )
     , m_caseData( nullptr )
     , m_mainGrid( nullptr )
+    , m_seabedTemperature( seabedTemperature )
+    , m_seabedDepth( seabedDepth )
 {
     if ( m_eclipseCase )
     {
@@ -60,16 +67,17 @@ RimFaultReactivationDataAccessorTemperature::~RimFaultReactivationDataAccessorTe
 //--------------------------------------------------------------------------------------------------
 void RimFaultReactivationDataAccessorTemperature::updateResultAccessor()
 {
-    if ( m_caseData )
+    if ( !m_caseData ) return;
-    {
+
-        RigEclipseResultAddress resVarAddress( RiaDefines::ResultCatType::DYNAMIC_NATIVE, "TEMP" );
+    RigEclipseResultAddress resVarAddress( RiaDefines::ResultCatType::DYNAMIC_NATIVE, "SOIL" );
-        m_eclipseCase->results( RiaDefines::PorosityModelType::MATRIX_MODEL )->ensureKnownResultLoaded( resVarAddress );
+    m_eclipseCase->results( RiaDefines::PorosityModelType::MATRIX_MODEL )->ensureKnownResultLoaded( resVarAddress );
-        m_resultAccessor = RigResultAccessorFactory::createFromResultAddress( m_caseData,
+    m_resultAccessor =
-                                                                              0,
+        RigResultAccessorFactory::createFromResultAddress( m_caseData, 0, RiaDefines::PorosityModelType::MATRIX_MODEL, m_timeStep, resVarAddress );
-                                                                              RiaDefines::PorosityModelType::MATRIX_MODEL,
+
-                                                                              m_timeStep,
+    auto [wellPaths, extractors] =
-                                                                              resVarAddress );
+        RimFaultReactivationDataAccessorWellLogExtraction::createEclipseWellPathExtractors( *m_model, *m_caseData, m_seabedDepth );
-    }
+    m_wellPaths  = wellPaths;
+    m_extractors = extractors;
 }
 
 //--------------------------------------------------------------------------------------------------
@@ -92,11 +100,19 @@ double RimFaultReactivationDataAccessorTemperature::valueAtPosition( const cvf::
 {
     if ( ( m_mainGrid != nullptr ) && m_resultAccessor.notNull() )
     {
-        auto cellIdx = m_mainGrid->findReservoirCellIndexFromPoint( position );
+        CAF_ASSERT( m_extractors.find( gridPart ) != m_extractors.end() );
-        if ( cellIdx != cvf::UNDEFINED_SIZE_T )
+        auto extractor = m_extractors.find( gridPart )->second;
-        {
+
-            return m_resultAccessor->cellScalar( cellIdx );
+        CAF_ASSERT( m_wellPaths.find( gridPart ) != m_wellPaths.end() );
-        }
+        auto wellPath = m_wellPaths.find( gridPart )->second;
+
+        auto [values, intersections] =
+            RimFaultReactivationDataAccessorWellLogExtraction::extractValuesAndIntersections( *m_resultAccessor.p(), *extractor.p(), *wellPath );
+
+        auto [value, pos] =
+            RimFaultReactivationDataAccessorWellLogExtraction::calculateTemperature( intersections, values, position, m_seabedTemperature );
+
+        return value;
     }
 
     return std::numeric_limits<double>::infinity();

ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorTemperature.h

@@ -28,6 +28,8 @@ class RimEclipseCase;
 class RigEclipseCaseData;
 class RigMainGrid;
 class RigResultAccessor;
+class RigWellPath;
+class RigEclipseWellLogExtractor;
 
 //==================================================================================================
 ///
@@ -36,7 +38,7 @@ class RigResultAccessor;
 class RimFaultReactivationDataAccessorTemperature : public RimFaultReactivationDataAccessor
 {
 public:
-    RimFaultReactivationDataAccessorTemperature( RimEclipseCase* eclipseCase );
+    RimFaultReactivationDataAccessorTemperature( RimEclipseCase* eclipseCase, double seabedTemperature, double seabedDepth );
     ~RimFaultReactivationDataAccessorTemperature();
 
     bool isMatching( RimFaultReactivation::Property property ) const override;
@@ -51,8 +53,14 @@ public:
 private:
     void updateResultAccessor() override;
 
-    RimEclipseCase*             m_eclipseCase;
+    RimEclipseCase*     m_eclipseCase;
-    RigEclipseCaseData*         m_caseData;
+    RigEclipseCaseData* m_caseData;
-    const RigMainGrid*          m_mainGrid;
+    const RigMainGrid*  m_mainGrid;
+    double              m_seabedTemperature;
+    double              m_seabedDepth;
+
     cvf::ref<RigResultAccessor> m_resultAccessor;
+
+    std::map<RimFaultReactivation::GridPart, cvf::ref<RigWellPath>>                m_wellPaths;
+    std::map<RimFaultReactivation::GridPart, cvf::ref<RigEclipseWellLogExtractor>> m_extractors;
 };

ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorWellLogExtraction.cpp

@@ -63,13 +63,37 @@ RimFaultReactivationDataAccessorWellLogExtraction::~RimFaultReactivationDataAcce
 ///
 //--------------------------------------------------------------------------------------------------
 std::pair<double, cvf::Vec3d> RimFaultReactivationDataAccessorWellLogExtraction::calculatePorBar( const std::vector<cvf::Vec3d>& intersections,
-                                                                                                  std::vector<double>&           values,
+                                                                                                  std::vector<double>& values,
-                                                                                                  const cvf::Vec3d&              position,
+                                                                                                  const cvf::Vec3d&    position,
-                                                                                                  double                         gradient )
+                                                                                                  double               gradient )
 {
     // Fill in missing values
-    fillInMissingValues( intersections, values, gradient );
+    fillInMissingValuesWithGradient( intersections, values, gradient );
+    return findValueAndPosition( intersections, values, position );
+}
 
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+std::pair<double, cvf::Vec3d>
+    RimFaultReactivationDataAccessorWellLogExtraction::calculateTemperature( const std::vector<cvf::Vec3d>& intersections,
+                                                                             std::vector<double>&           values,
+                                                                             const cvf::Vec3d&              position,
+                                                                             double                         seabedTemperature )
+{
+    // Fill in missing values
+    fillInMissingValuesWithTopValue( intersections, values, seabedTemperature );
+    return findValueAndPosition( intersections, values, position );
+}
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+std::pair<double, cvf::Vec3d>
+    RimFaultReactivationDataAccessorWellLogExtraction::findValueAndPosition( const std::vector<cvf::Vec3d>& intersections,
+                                                                             const std::vector<double>&     values,
+                                                                             const cvf::Vec3d&              position )
+{
     // Linear interpolation between two points
     auto lerp = []( const cvf::Vec3d& start, const cvf::Vec3d& end, double t ) { return start + t * ( end - start ); };
 
@@ -99,7 +123,7 @@ std::pair<double, cvf::Vec3d> RimFaultReactivationDataAccessorWellLogExtraction:
 ///
 //--------------------------------------------------------------------------------------------------
 std::pair<int, int> RimFaultReactivationDataAccessorWellLogExtraction::findIntersectionsForTvd( const std::vector<cvf::Vec3d>& intersections,
-                                                                                                double                         tvd )
+                                                                                                double tvd )
 {
     int topIdx    = -1;
     int bottomIdx = -1;
@@ -109,7 +133,7 @@ std::pair<int, int> RimFaultReactivationDataAccessorWellLogExtraction::findInter
         {
             auto top    = intersections[i - 1];
             auto bottom = intersections[i];
-            if ( top.z() > tvd && bottom.z() < tvd )
+            if ( top.z() >= tvd && bottom.z() < tvd )
             {
                 topIdx    = static_cast<int>( i ) - 1;
                 bottomIdx = static_cast<int>( i );
@@ -153,11 +177,11 @@ std::pair<int, int> RimFaultReactivationDataAccessorWellLogExtraction::findOverb
 //--------------------------------------------------------------------------------------------------
 ///
 //--------------------------------------------------------------------------------------------------
-double RimFaultReactivationDataAccessorWellLogExtraction::computePorBarWithGradient( const std::vector<cvf::Vec3d>& intersections,
+double RimFaultReactivationDataAccessorWellLogExtraction::computeValueWithGradient( const std::vector<cvf::Vec3d>& intersections,
-                                                                                     const std::vector<double>&     values,
+                                                                                    const std::vector<double>&     values,
-                                                                                     int                            i1,
+                                                                                    int                            i1,
-                                                                                     int                            i2,
+                                                                                    int                            i2,
-                                                                                     double                         gradient )
+                                                                                    double                         gradient )
 {
     double tvdDiff = intersections[i2].z() - intersections[i1].z();
     return tvdDiff * gradient + values[i2];
@@ -166,50 +190,51 @@ double RimFaultReactivationDataAccessorWellLogExtraction::computePorBarWithGradi
 //--------------------------------------------------------------------------------------------------
 ///
 //--------------------------------------------------------------------------------------------------
-void RimFaultReactivationDataAccessorWellLogExtraction::fillInMissingValues( const std::vector<cvf::Vec3d>& intersections,
+void RimFaultReactivationDataAccessorWellLogExtraction::fillInMissingValuesWithGradient( const std::vector<cvf::Vec3d>& intersections,
-                                                                             std::vector<double>&           values,
+                                                                                         std::vector<double>&           values,
-                                                                             double                         gradient )
+                                                                                         double                         gradient )
 {
     CAF_ASSERT( intersections.size() == values.size() );
 
     auto calculatePorePressure = []( double depth, double gradient )
     { return RiaEclipseUnitTools::pascalToBar( gradient * 9.81 * depth * 1000.0 ); };
 
-    auto computeGradient = []( double depth1, double value1, double depth2, double value2 )
-    { return ( value2 - value1 ) / ( depth2 - depth1 ); };
-
     auto [lastOverburdenIndex, firstUnderburdenIndex] = findOverburdenAndUnderburdenIndex( values );
 
     // Fill in overburden values using gradient
     double topPorePressure = calculatePorePressure( std::abs( intersections[0].z() ), gradient );
-    double overburdenGradient =
+    insertOverburdenValues( intersections, values, lastOverburdenIndex, topPorePressure );
-        computeGradient( intersections[0].z(), topPorePressure, intersections[lastOverburdenIndex].z(), values[lastOverburdenIndex] );
-
-    for ( int i = 0; i < lastOverburdenIndex; i++ )
-    {
-        values[i] = computePorBarWithGradient( intersections, values, i, lastOverburdenIndex, -overburdenGradient );
-    }
 
     // Fill in underburden values using gradient
-    int    lastElementIndex    = static_cast<int>( values.size() ) - 1;
+    int    lastElementIndex   = static_cast<int>( values.size() ) - 1;
-    double bottomPorePressure  = calculatePorePressure( std::abs( intersections[lastElementIndex].z() ), gradient );
+    double bottomPorePressure = calculatePorePressure( std::abs( intersections[lastElementIndex].z() ), gradient );
-    double underburdenGradient = computeGradient( intersections[firstUnderburdenIndex].z(),
+    insertUnderburdenValues( intersections, values, firstUnderburdenIndex, bottomPorePressure );
-                                                  values[firstUnderburdenIndex],
-                                                  intersections[lastElementIndex].z(),
-                                                  bottomPorePressure );
-
-    for ( int i = lastElementIndex; i >= firstUnderburdenIndex; i-- )
-    {
-        values[i] = computePorBarWithGradient( intersections, values, i, firstUnderburdenIndex, -underburdenGradient );
-    }
 
     // Interpolate the missing values (should only be intra-reservoir by now)
-    std::vector<double> intersectionsZ;
+    std::vector<double> intersectionsZ = extractDepthValues( intersections );
-    for ( auto i : intersections )
+    RiaInterpolationTools::interpolateMissingValues( intersectionsZ, values );
-    {
+}
-        intersectionsZ.push_back( i.z() );
-    }
 
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+void RimFaultReactivationDataAccessorWellLogExtraction::fillInMissingValuesWithTopValue( const std::vector<cvf::Vec3d>& intersections,
+                                                                                         std::vector<double>&           values,
+                                                                                         double                         topValue )
+{
+    CAF_ASSERT( intersections.size() == values.size() );
+
+    auto [lastOverburdenIndex, firstUnderburdenIndex] = findOverburdenAndUnderburdenIndex( values );
+
+    // Fill in overburden values using gradient
+    insertOverburdenValues( intersections, values, lastOverburdenIndex, topValue );
+
+    // Fill in underburden values
+    double bottomValue = values[firstUnderburdenIndex - 1];
+    insertUnderburdenValues( intersections, values, firstUnderburdenIndex, bottomValue );
+
+    // Interpolate the missing values (should only be intra-reservoir by now)
+    std::vector<double> intersectionsZ = extractDepthValues( intersections );
     RiaInterpolationTools::interpolateMissingValues( intersectionsZ, values );
 }
 
@@ -218,10 +243,11 @@ void RimFaultReactivationDataAccessorWellLogExtraction::fillInMissingValues( con
 //--------------------------------------------------------------------------------------------------
 std::vector<cvf::Vec3d> RimFaultReactivationDataAccessorWellLogExtraction::generateWellPoints( const cvf::Vec3d& faultTopPosition,
                                                                                                const cvf::Vec3d& faultBottomPosition,
+                                                                                               double            seabedDepth,
                                                                                                const cvf::Vec3d& offset )
 {
     cvf::Vec3d faultTop = faultTopPosition + offset;
-    cvf::Vec3d seabed( faultTop.x(), faultTop.y(), 0.0 );
+    cvf::Vec3d seabed( faultTop.x(), faultTop.y(), seabedDepth );
     cvf::Vec3d faultBottom = faultBottomPosition + offset;
     cvf::Vec3d underburdenBottom( faultBottom.x(), faultBottom.y(), -10000.0 );
     return { seabed, faultTop, faultBottom, underburdenBottom };
@@ -250,7 +276,8 @@ std::vector<double> RimFaultReactivationDataAccessorWellLogExtraction::generateM
 //--------------------------------------------------------------------------------------------------
 std::pair<std::map<RimFaultReactivation::GridPart, cvf::ref<RigWellPath>>, std::map<RimFaultReactivation::GridPart, cvf::ref<RigEclipseWellLogExtractor>>>
     RimFaultReactivationDataAccessorWellLogExtraction::createEclipseWellPathExtractors( const RigFaultReactivationModel& model,
-                                                                                        RigEclipseCaseData&              eclipseCaseData )
+                                                                                        RigEclipseCaseData&              eclipseCaseData,
+                                                                                        double                           seabedDepth )
 {
     auto [faultTopPosition, faultBottomPosition] = model.faultTopBottom();
     auto   faultNormal                           = model.faultNormal();
@@ -265,6 +292,7 @@ std::pair<std::map<RimFaultReactivation::GridPart, cvf::ref<RigWellPath>>, std::
         std::vector<cvf::Vec3d> wellPoints =
             RimFaultReactivationDataAccessorWellLogExtraction::generateWellPoints( faultTopPosition,
                                                                                    faultBottomPosition,
+                                                                                   seabedDepth,
                                                                                    sign * faultNormal * distanceFromFault );
         cvf::ref<RigWellPath> wellPath =
             new RigWellPath( wellPoints, RimFaultReactivationDataAccessorWellLogExtraction::generateMds( wellPoints ) );
@@ -302,3 +330,61 @@ std::pair<std::vector<double>, std::vector<cvf::Vec3d>>
 
     return { values, intersections };
 }
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+double RimFaultReactivationDataAccessorWellLogExtraction::computeGradient( double depth1, double value1, double depth2, double value2 )
+{
+    return ( value2 - value1 ) / ( depth2 - depth1 );
+}
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+std::vector<double> RimFaultReactivationDataAccessorWellLogExtraction::extractDepthValues( const std::vector<cvf::Vec3d>& intersections )
+{
+    std::vector<double> intersectionsZ;
+    for ( auto i : intersections )
+    {
+        intersectionsZ.push_back( i.z() );
+    }
+    return intersectionsZ;
+}
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+void RimFaultReactivationDataAccessorWellLogExtraction::insertUnderburdenValues( const std::vector<cvf::Vec3d>& intersections,
+                                                                                 std::vector<double>&           values,
+                                                                                 int                            firstUnderburdenIndex,
+                                                                                 double                         bottomValue )
+{
+    int    lastElementIndex    = static_cast<int>( values.size() ) - 1;
+    double underburdenGradient = computeGradient( intersections[firstUnderburdenIndex].z(),
+                                                  values[firstUnderburdenIndex],
+                                                  intersections[lastElementIndex].z(),
+                                                  bottomValue );
+
+    for ( int i = lastElementIndex; i >= firstUnderburdenIndex; i-- )
+    {
+        values[i] = computeValueWithGradient( intersections, values, i, firstUnderburdenIndex, -underburdenGradient );
+    }
+}
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+void RimFaultReactivationDataAccessorWellLogExtraction::insertOverburdenValues( const std::vector<cvf::Vec3d>& intersections,
+                                                                                std::vector<double>&           values,
+                                                                                int                            lastOverburdenIndex,
+                                                                                double                         topValue )
+{
+    double overburdenGradient =
+        computeGradient( intersections[0].z(), topValue, intersections[lastOverburdenIndex].z(), values[lastOverburdenIndex] );
+
+    for ( int i = 0; i < lastOverburdenIndex; i++ )
+    {
+        values[i] = computeValueWithGradient( intersections, values, i, lastOverburdenIndex, -overburdenGradient );
+    }
+}

ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorWellLogExtraction.h

@@ -43,12 +43,18 @@ public:
     static std::pair<double, cvf::Vec3d>
         calculatePorBar( const std::vector<cvf::Vec3d>& intersections, std::vector<double>& values, const cvf::Vec3d& position, double gradient );
 
+    static std::pair<double, cvf::Vec3d> calculateTemperature( const std::vector<cvf::Vec3d>& intersections,
+                                                               std::vector<double>&           values,
+                                                               const cvf::Vec3d&              position,
+                                                               double                         seabedTemperature );
     static std::pair<std::map<RimFaultReactivation::GridPart, cvf::ref<RigWellPath>>,
                      std::map<RimFaultReactivation::GridPart, cvf::ref<RigEclipseWellLogExtractor>>>
-        createEclipseWellPathExtractors( const RigFaultReactivationModel& model, RigEclipseCaseData& eclipseCaseData );
+        createEclipseWellPathExtractors( const RigFaultReactivationModel& model, RigEclipseCaseData& eclipseCaseData, double seabedDepth );
 
-    static std::vector<cvf::Vec3d>
+    static std::vector<cvf::Vec3d> generateWellPoints( const cvf::Vec3d& faultTopPosition,
-        generateWellPoints( const cvf::Vec3d& faultTopPosition, const cvf::Vec3d& faultBottomPosition, const cvf::Vec3d& offset );
+                                                       const cvf::Vec3d& faultBottomPosition,
+                                                       double            seabedDepth,
+                                                       const cvf::Vec3d& offset );
 
     static std::vector<double> generateMds( const std::vector<cvf::Vec3d>& points );
 
@@ -59,10 +65,26 @@ public:
 protected:
     static std::pair<int, int> findIntersectionsForTvd( const std::vector<cvf::Vec3d>& intersections, double tvd );
     static std::pair<int, int> findOverburdenAndUnderburdenIndex( const std::vector<double>& values );
-    static double              computePorBarWithGradient( const std::vector<cvf::Vec3d>& intersections,
+    static double              computeValueWithGradient( const std::vector<cvf::Vec3d>& intersections,
-                                                          const std::vector<double>&     values,
+                                                         const std::vector<double>&     values,
-                                                          int                            i1,
+                                                         int                            i1,
-                                                          int                            i2,
+                                                         int                            i2,
-                                                          double                         gradient );
+                                                         double                         gradient );
-    static void fillInMissingValues( const std::vector<cvf::Vec3d>& intersections, std::vector<double>& values, double gradient );
+    static void fillInMissingValuesWithGradient( const std::vector<cvf::Vec3d>& intersections, std::vector<double>& values, double gradient );
+    static void fillInMissingValuesWithTopValue( const std::vector<cvf::Vec3d>& intersections, std::vector<double>& values, double topValue );
+
+    static std::pair<double, cvf::Vec3d>
+        findValueAndPosition( const std::vector<cvf::Vec3d>& intersections, const std::vector<double>& values, const cvf::Vec3d& position );
+
+    static double              computeGradient( double depth1, double value1, double depth2, double value2 );
+    static std::vector<double> extractDepthValues( const std::vector<cvf::Vec3d>& intersections );
+
+    static void insertUnderburdenValues( const std::vector<cvf::Vec3d>& intersections,
+                                         std::vector<double>&           values,
+                                         int                            firstUnderburdenIndex,
+                                         double                         bottomValue );
+    static void insertOverburdenValues( const std::vector<cvf::Vec3d>& intersections,
+                                        std::vector<double>&           values,
+                                        int                            lastOverburdenIndex,
+                                        double                         topValue );
 };

ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationModel.cpp

@@ -127,6 +127,7 @@ RimFaultReactivationModel::RimFaultReactivationModel()
 
     CAF_PDM_InitField( &m_waterDensity, "WaterDensity", 1030.0, "Water Density [kg/m3]" );
     CAF_PDM_InitField( &m_frictionAngleDeg, "FrictionAngle", 20.0, "Friction Angle [degree]" );
+    CAF_PDM_InitField( &m_seabedTemperature, "SeabedTemperature", 5.0, "Seabed Temperature [C]" );
 
     CAF_PDM_InitFieldNoDefault( &m_targets, "Targets", "Targets" );
     m_targets.uiCapability()->setUiEditorTypeName( caf::PdmUiTableViewEditor::uiEditorTypeName() );
@@ -480,6 +481,11 @@ void RimFaultReactivationModel::defineUiOrdering( QString uiConfigName, caf::Pdm
         propertiesGrp->add( &m_useGridElasticProperties );
         propertiesGrp->add( &m_useGridStress );
         propertiesGrp->add( &m_waterDensity );
+
+        propertiesGrp->add( &m_seabedTemperature );
+
+        bool useTemperatureFromGrid = m_useGridTemperature();
+        m_seabedTemperature.uiCapability()->setUiReadOnly( !useTemperatureFromGrid );
     }
 
     propertiesGrp->add( &m_frictionAngleDeg );
@@ -790,3 +796,11 @@ double RimFaultReactivationModel::waterDensity() const
 {
     return m_waterDensity;
 }
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+double RimFaultReactivationModel::seabedTemperature() const
+{
+    return m_seabedTemperature;
+}

ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationModel.h

@@ -128,6 +128,7 @@ public:
     double seaBedDepth() const;
     double waterDensity() const;
     double frictionAngleDeg() const;
+    double seabedTemperature() const;
 
     RimEclipseCase* eclipseCase() const;
     RimGeoMechCase* geoMechCase() const;
@@ -182,6 +183,7 @@ private:
 
     caf::PdmField<double> m_waterDensity;
     caf::PdmField<double> m_frictionAngleDeg;
+    caf::PdmField<double> m_seabedTemperature;
 
     caf::PdmField<size_t> m_startCellIndex;
     caf::PdmField<int>    m_startCellFace;