Quantitative seismic interpretation of the Lower Cretaceous reservoirs in the Valdemar Field, Danish North Sea

A quantitative seismic interpretation study is presented for the Lower Cretaceous Tuxen reservoir in the Valdemar Field, which is associated with heterogeneous and complex geology. Our objective is to better outline the reservoir quality variations of the Tuxen reservoir across the Valdemar Field. Seismic pre-stack data and well logs from two appraisal wells form the basis of this study. The workflow used includes seismic and rock physics forward modelling, attribute analysis, a coloured inversion, and a Bayesian pre-stack inversion for litho-fluid classification. Based on log data, the rock physics properties of the Tuxen interval reveal that the seismic signal is more governed by porosity than water-saturation changes at near-offset (or small angle). The coloured and Bayesian inversion results were generally consistent with well-log observations at the reservoir level and conformed to interpreted horizons. Although the available data have some limitations and the geological setting is complex, the results implied more promising reservoir quality in some areas than others. Hence, the results may offer useful information for delineating the best reservoir zones for further field development and selecting appropriate production strategies.

The prospectivity of Lower Cretaceous reservoirs in the Central North Sea has received relatively little attention due to its complex and heterogeneous geology. The Lower Cretaceous consists of alternating chalk, marly chalk, marlstone and claystone layers, and is often associated with tight reservoirs. Nevertheless, in the Valdemar Field in the Danish North Sea, hydrocarbons have been produced from Lower Cretaceous reservoirs in recent decades. The recovery factor, however, has been relatively low and uneven across the field, and the distribution of reservoir quality is poorly outlined and understood (Jakobsen et al. 2004(Jakobsen et al. , 2005. Hence, field investigations using quantitative seismic interpretation tools can yield supplementary reservoir information to assess and reduce associated risks during field development and reservoir management (Avseth et al. 2005;Dvorkin et al. 2014).
A key challenge in further unlocking the potential of the Lower Cretaceous reservoir is the seismic image from the target interval characterized by low-resolution, interfering noise and variable reflectivity. This is partly due to interfering multiples and converted waves propagated from the overlying high-velocity chalk package (Vidalie et al. 2014;Montazeri et al. 2018), which is characterized by strongly contrasting layers and structural complexities such as faults, fractures and slumping. Lateral reservoir heterogeneity and overburden variabilities also make it more challenging to seismically interpret the reservoir targets.
From previous studies of the Valdemar Field, we find the Jakobsen et al. (2004) study where well logs were correlated with core samples to present a detailed zonation of the Lower Cretaceous reservoirs. Calvert et al. (2018) interpreted hydrocarbon production effects based on a 4D monitoring survey from Valdemar and the impact on the trajectory of a planned production well. Vidalie et al. (2014) presented a regional seismic stratigraphic analysis of the mudstones within the Valhall Formation in the Valdemar and Adda fields. Bredesen et al. (2020) presented a rock physics feasibility study to investigate the elastic properties of Lower Cretaceous reservoirs in the Valdemar Field based on well-log data.
In this paper, we present a quantitative seismic interpretation study to outline the reservoir quality variations of the Lower Cretaceous Tuxen reservoir across the Bo and Bo South areas in the Valdemar Field. Based on well-log data, the rock physics relationship between the elastic properties, such as the acoustic impedance (AI) and Poisson's ratio (or V P /V S ), and changes in porosity, water saturation and lithology-fluid classes (LFCs) are studied for the Tuxen reservoir. Our objective is to investigate whether these reservoir variations can be seismically distinguished based on seismic amplitude v. offset/angle (AVO/AVA) data. We approach this problem by deriving seismic AVA attributes and inversion results that have the potential to improve the reservoir resolution from standard post-stack seismic images. First, we review the geological setting and data coverage of the Valdemar Field. Then, we use well-log data to investigate the rock physics relationships between reservoir properties and seismic properties, perform a rock physics sensitivity study of relative changes in porosity and water saturation, and investigate the tuning effects from varying overburden thicknesses. Next, we investigate variations in the seismic AVA intercept and gradient attributes that may be associated with reservoir changes (Dvorkin et al. 2014;Simm et al. 2014). Similarly, a coloured inversion (Lancaster and Whitcombe 2000;Blouin and Gloaguen 2017) is evaluated before an absolute Bayesian pre-stack inversion is presented that shows how we can distinguish litho-fluid classes representing various Tuxen reservoir units associated with different reservoir quality (Ndingwan et al. 2018;Kolbjørnsen et al. 2020). Finally, we discuss some of the shortcomings and uncertainties of the main results.

Regional setting
The Valdemar Field is located in the central part of the Danish Central Graben, north of the Salt Dome Province (Fig. 1). The area is characterized by the north-south-elongated North Jens-Bo inversion anticline (Vejbaek and Andersen 2002;Japsen et al. 2003;Hansen et al. 2021) with two structural highs separated by a saddle point. The field can be subdivided into three main areas: (1) North Jens, (2) Bo and (3) Bo South. The reservoir section is divided into the Sola Formation and the Tuxen Formation within the Lower Cretaceous Cromer Knoll Group (Jakobsen et al. 2004) (Fig. 2).
The reservoirs were deposited from Early Cretaceous structural terraces, saddles or minor highs that provided good conditions for the accumulation of economically attractive thick reservoir  successions (Ineson 1993). The Sola and Tuxen formations form multilayered, heterogeneous reservoirs composed of chalk, marly chalk, organic-rich marlstone, and marlstone deposited primarily by suspension settling of pelagic carbonate and hemipelagic mud. Bedding is typically on a decimetre-metre scale and is cyclic at given levels. In the Tuxen Formation, which is the reservoir target in this study, the Munk Marl (MM) bed (Fig. 2) divides into upper and lower reservoir units, where the higher reservoir quality is typically observed within the upper part. More details on the Lower Cretaceous facies development and reservoir properties in the Danish area can be found in Ineson (1993), Copestake et al. (2003), Jakobsen et al. (2004Jakobsen et al. ( , 2005 and van Buchem et al. (2017).
The Valdemar Field was discovered in 1977, with production commencing in 1993 and 2007 from the North Jens and Bo areas, respectively (Calvert et al. 2018). The field today is the fourth largest oil-producing field managed by the Danish Underground Consortium (DUC), and is produced from the North Jens and Bo areas; Bo South is not developed. The current burial depths of the top reservoir are c. 2200 and 2500 m in the North Jens and Bo South areas, respectively. Consequently, the reservoir quality downgrades southwards. The North Jens area is prognosed to recover twice as many hydrocarbons as the Bo area. The Tuxen Formation is the most prominent reservoir section, with a thickness of around 50 m. The porosity-depth relationship for the Tuxen reservoir indicates that porosities less than 15% are to be expected at depths greater than 2500-3000 m (Jakobsen et al. 2005). In the Valdemar Field, the Lower Cretaceous reservoirs are typically characterized by low matrix permeability, even for high porosities, but the presence of natural fractures tend to influence the flow properties by enhancing the effective permeability (Glad et al. 2020).

Data coverage and conditioning
The seismic reflection data used in this study is from a time-lapse (4D) survey collected in 2016, Kirchhoff pre-stack depth-migrated (PSDM) data (Calvert et al. 2018). The data acquisition used a 50 m separation for the source shots and eight streamers separated 100 m apart, each with 240 groups at 12.5 m intervals, leading to a maximum 3 km offset. The geometry is set-up with a bin size of 12.5 m in the inline and cross-line directions, and a with a time sampling of 4 ms. The initial processing steps were: • 3 Hz low-cut filter; • attenuation of swell, impulse and linear noise; • designature to zero phase and de-bubble; • despike; • first radon multiple removal; • bulk amplitude time and phase match to baseline survey; • tidal statics; • inverse Q-compensation ( phase only); • TTI Kirchhoff pre-stack depth migration; • stacking velocity analysis; • second radon multiple removal; • normal moveout correction; and • wide-angle mute.
The data have reverse polarity where a negative sample corresponds to an increase in AI.
We focus on the Bo and Bo South areas where the reservoir quality distribution is not well known (Fig. 1). Figure 3 shows a seismic cross-section intersecting the BO-2X and BO-3X wells in the southern Valdemar Field, including interpretations of Lower Cretaceous units. The Lower Cretaceous interval is characterized by a low signal/noise ratio and poor continuity of reflection events. Hence, choosing a robust interpretation strategy is challenging. Key horizons for the Lower Cretaceous interval were obtained from a previous regional interpretation study (Smit et al. 2021). The Lower Cretaceous interval is bounded by the Base Chalk and Base Cretaceous Unconformity (BCU) events.
To obtain a reliable seismic pre-stack amplitude analysis or inversion, it is crucial to quality check the data and, subsequently, to reduce any errors encountered. There are several potential datarelated problems that can corrupt the amplitude analysis and inversions, such as non-flat reflection events, multiples mimicking primary reflection events ( particularly on near offsets), deviating frequencies at near and far offsets, random noise, amplitude variations with offset or angle not related to geology, and limited frequency bandwidth. For the Valdemar data, the following key post-migration data-processing steps were applied to further condition the data for amplitude v. offset (AVO) analysis and inversion: parabolic and linear radon; residual moveout correction; gather flattening; angle-dependent spectral balancing; and edgeand coherency-enhancing anisotropic diffusion filter. An example of quality checking the post-migration processing is shown in Figure 4. It shows a cross-correlation in a 100 ms window across the Top Tuxen horizon between near (5°-15°) and far (25°-35°) stacks for the same arbitrary line as in Figure 3, before and after conditioning. If the various offset stacks were identical (no moveout problem or noise interference), the cross-correlation would be a peak at time zero. Hence, the conditioning has certainly improved the near-far stack correlation at the reservoir level.

Seismic modelling, tuning and rock physics sensitivity analysis
A suite of well logs next to synthetic and real zero-phase pre-stack gathers along the BO-2X well is shown in Figure 5 covering the upper part of the Lower Cretaceous. An increase in AI is represented by negative amplitudes, shown in red colours. As a first approach, the synthetic gather was generated by convolving Zoeppritz reflection coefficients with the zero-phase statistical wavelet in Figure 5g extracted from the Lower Cretaceous interval with the frequency spectrum in Figure 5h. Some visual discrepancies in amplitude and travel time between the synthetics and real gathers are observed. This is partly because the statistical wavelet parameters, such as phase rotation and scaling, were not calibrated at this stage. A deterministic wavelet from a well-tie procedure for the inversion  is presented later in this paper. Nevertheless, the Top Tuxen event can be recognized by the blue peak that appears over a wider time interval in the real gather compared to the synthetic gather. Figure 6 shows the same plot for BO-3X, where the real gather appears noisier. The elastic contrasts and seismic response for the Tuxen reservoir are different in the two wells. In BO-2X, the Top Tuxen exhibits a distinct decrease in AI and V P /V S ; whereas in BO-3X, the Top Tuxen yields more modest elastic contrasts. The spatially varying elastic contrasts of the Lower Cretaceous yield differing seismic responses that make it challenging to resolve and map the Tuxen reservoir across the Valdemar Field. Hence, the regional interpretations of the Lower Cretaceous surfaces are mainly driven by manual interpretations guided by well ties (Fig. 3) rather than following a distinctive seismic signature (Smit et al. 2021). Figure 7a-d shows 1D blocky models for the AI and V P /V S logs based on the in situ reservoir conditions for the BO-2X and BO-3X  wells, and the corresponding synthetic gathers generated with the statistical wavelet (Fig. 7e). Top Tuxen is defined at two-way time (TWT) = 1 s. The synthetic responses derived directly from the well-log samples (Figs 5 and 6) and from the 1D blocky models are generally in a good agreement. To investigate how varying reservoir properties and pore-fluid saturations affect the seismic response of the Tuxen reservoir, a rock physics modelling sensitivity study was performed. The main reservoir variations within the Tuxen interval are in the porosity (f) and the water saturation (S w ), for example, above and below the Munk Marl (MM) bed (Figs 5 and 6), whereas the shale volume is relatively constant. We used the 1D blocky model from BO-2X in Figure 7 and implemented the iso-frame (IF) rock physics model (Fabricius et al. 2005(Fabricius et al. , 2007 as proposed by Bredesen et al. (2020) (see Appendix A for details) to calculate and replace the AI and V P /V S properties within the Tuxen interval as functions of varying porosity and water saturation. Figure 8a-e shows blocky logs and the corresponding synthetic gathers for three different Tuxen reservoir scenarios. From the initial scenario, assuming f ¼ 30% and S w ¼ 100%, the Df scenario represents a relative 50% porosity reduction to f ¼ 15% and keeps the remaining model parameters fixed. Likewise, DS w assumes a 50% water-saturation reduction to S w ¼ 50%. The shale volume was set to 10% in all scenarios. Figure 8g shows the corresponding amplitude v. angle (AVA) for Top Tuxen for the three scenarios, and Figure 8h shows the corresponding AVA intercept and gradient attributes from a two-term approximation (Shuey 1985) fitted to the 5°-30°range. Relative to the initial scenario, the Df scenario changes the AI contrast substantially at the Top Tuxen interface from a decrease to an increase, which produces an opposite polarity response at the Top Tuxen reflection event (Fig. 8d) and changes the AVA intercept sign (Fig. 8g, h). Correspondingly, the DS w scenario yields a decrease in both AI and V P /V S , resulting in an increasing AVA intercept and gradient. These results indicate that in a relative sense, porosity has a more significant impact on the seismic response than water-saturation changes.
The Tuxen Formation is c. 50 m thick in both wells, whereas the Sola Formation is 17 and 30 m thick in BO-2X and BO-3X, respectively. Hence, we investigated tuning effects from thickness variations of the Sola Formation. Figure 9a shows a wedge model for the Sola Formation where the thickness varies between 0 and 60 m along the x-axis based on the blocky model from BO-2X (Fig. 7). The amplitude values extracted along the Top and Base Tuxen events are shown in Figure 9b, where the base of the model represents the constant properties of the Valhall Formation. The wavelet side lobes for the Top Sola event start to interfere with the Top Tuxen event from around 30 m, and maximum tuning occurs at 21 m thickness (red dot), which is close to the Sola thickness observed in BO-2X. However, the strongest interference is only −16% from the amplitude at 60 m thickness, and the zero-phase statistical wavelet has symmetrical side lobes of c. 6 ms length with c. 1/5 amplitude relative to the centre trough. Comparing this result with the AVA characteristics in Figure 8, the DS w and Df scenarios yield, respectively, +84% and −268% changes in the zero-angle  amplitude, which equals the AVA intercept, relative to the initial scenario. Therefore, the modelled tuning effects have a relatively small impact on the seismic amplitudes. In addition, because the Sola Formation is thicker in BO-3X than in BO-2X, the tuning effects will probably be even less of an issue downflank towards Bo South.

Seismic AVA analysis
In this section, variations in seismic amplitude v. angle (AVA) attributes are investigated that could be linked to reservoir changes. In the blocky model in Figure 7, the Top Tuxen interface exhibits a decrease in both AI and V P /V S . The corresponding AVA behaviour for the two wells are plotted in Figure 10, where positive amplitudes are seen at the nears that increase with the angle, corresponding to positive intercept and gradient AVA attributes (Shuey 1985). To calculate the AVA intercept and gradient from the seismic pre-stack data, a two-term AVA approximation (Shuey 1985) was fitted to the data within a 5°-30°angle range to the pre-stack data in the angle domain. Figure 11a and b shows the AVA intercept and gradient extracted along the Top Tuxen horizon. Figure 12 shows the corresponding AVA intercept v. gradient cross-plot from the Top Tuxen horizon, where the black dashed line represents the background trend with a 16°rotation angle.
Ultimately, we want to investigate deviations in the data away from the background trend of Top Tuxen across the study area. This is because variations in reservoir properties and pore fluids may cause the AVA data to plot some distance away from the background trend as AVA variabilities. Figure 11c shows the perpendicular distance away from the background trend as a measure of data variability. Based on the AVA modelling in Figure 8h, the data plotting more towards the upper-right corner in Figure 12 (i.e. increasing positive values for the perpendicular distance) could be associated with higher hydrocarbon saturations. However, this interpretation relies on a fairly crude assumption that the reservoir properties, such as porosity, will be fairly constant across the Top Tuxen horizon, which is not the case over larger distances. It is still interesting to see how the increasing perpendicular distance conforms with the structural high in the Bo area, where we would expect higher hydrocarbon saturations to occur. Similar anomalies are not seen in the Bo South area, and it is generally difficult to make confident interpretations of the Top Tuxen reservoir. Some possible explanations are that the interpreted Top Tuxen horizon is not accurately aligned with the top reservoir interface or the noise levels vary spatially across the study area. Also, notice the higher resolution of the intercept compared to the gradient, and that the intercept has a more prominent footprint in Figure 11c.

Seismic inversion
In this section, we go from interface to layer properties by evaluating: (1) a relative coloured inversion (CI) of the data for relative acoustic impedance (RAI) (i.e. the change in AI going from one layer to another); and (2) a Bayesian (or probabilistic) AVA inversion for litho-fluid classification. Seismic inversions can help to better balance the frequencies in the seismic data to improve resolution and to compensate for geophysical artefacts such as seismic tuning, providing a more detailed image for reservoir characterization.

Coloured inversion
A CI provides a quick and easily accessible approach for relative impedance inversions (Lancaster and Whitcombe 2000; Blouin and  Gloaguen 2017). The CI attempts to make the seismic data look more like acoustic layer properties as measured by acoustic well-log data. It does so by applying a CI operator that shapes the seismic spectrum to approach the well-log spectrum and rotates the phase by 90°. Because the CI does not integrate any prior information or a low-frequency model from the wells, it is entirely driven by the seismic data and, therefore, band limited. Figure 13a shows the RAI derived from applying the CI operator in Figure 13b to the AVA intercept attribute (Fig. 11a) along the arbitrary line. As the RAI is strongly correlated to porosity in the Tuxen reservoir (Fig. 8), the CI results may offer a useful porosity indicator. Figure 13c shows the frequency spectrum of the seismic intercept attribute (pink), the AI log from BO-2X within the Lower Cretaceous interval (grey) and the CI operator (brown). The CI spectrum shows a relative boost of low and high frequencies as it approaches the AI-log spectrum. The selection of CI parameters, such as the operator length and frequency bandwidth cutoffs, is a compromise of the seismic signal/noise ratio, where we try to avoid boosting low-and high-frequency noise. Interestingly, the CI implies a decrease in RAI (blue colour) between the Top Tuxen and Top Valhall horizons at the BO-2X location. The low-impedance layer follows the horizons southward and fades out downflank of the Bo structure. Around BO-3X, a similar anomaly appears above the Top Tuxen horizon that agrees with the AI reduction observed in BO-3X in the Sola Formation (Fig. 6). Notice that the relative porosity change when going from the Sola Formation into the Tuxen Formation is much lower in BO-3X (Fig. 6) than in BO-2X (Fig. 5), which can explain why there is no distinctive AI contrast at Top Tuxen in BO-3X. Figure 13d shows a RAI map along the Top Tuxen horizon with a small time shift applied to capture the CI results within the Tuxen reservoir interval. Notice the improvement in resolution of the Tuxen reservoir when compared to the input AVA intercept in Figure 11a. The strongest negative RAI values are predicted towards the Bo structural high and north towards the saddle point, and apparently correlates to lower burial depths with more preserved porosities.
To approach a more quantitative correlation between the RAI and the porosity, we used the calibrated iso-frame (IF) rock physics model (Appendix A). A two-term approximation (Shuey 1985) was used to model the AVA intercept with fixed properties for the top layer defined by mean properties from the Sola interval in BO-2X (V P = 2844 m s −1 , V S = 1413 m s −1 and density = 2.375 g cm −3 ). For the bottom layer, we used the IF rock physics model with varying porosity and with fixed S w = 0.4 and a shale volume equal to 0.1. Figure 14a shows the modelled AVA intercept as a function  of porosity used to subsequently correlate it with the RAI data. Figure 14b shows the calculated mean RAI between the Top Tuxen and Top Valhall horizons. Consequently, the RAI maps in Figures  13d and 14b deviate slightly as the former represents the Top Tuxen surface values, whereas the latter represents the Tuxen interval. Figure 14c shows the corresponding porosity derived from the mean RAI in Figure 14b. At first glance, high-porosity distributions are observed at the Bo structure as it correlates with the mean decrease in RAI within the Tuxen interval. Outside the Bo structural high and southwards, the mean RAI yields an increase, which is also seen in the RAI section in Figure 13a. However, the quantitative accuracy of the derived porosities in Figure 14c depends on several assumptions. For example, using fixed layer properties and rock physics parameters will not be representable for all areas across the whole field. Also, the RAI needs to be properly scaled by calibrating a scaling factor such that it can be correlated with absolute rock physics properties (Jensen et al. 2016). Nevertheless, the maps can be used as a rough porosity screening indicator for the Tuxen interval.

Bayesian pre-stack inversion for litho-fluid classification
The CI results yield a relative AI inversion across layer boundaries and therefore they cannot be used to interpret the absolute reservoir properties for specific layers. Furthermore, as we are dealing with noisy and low-resolution seismic data from complex geology, it makes sense to incorporate some prior information to guide the inversion where the data contain little reservoir information. This calls for investigation of more rigorous pre-stack inversion methods. Therefore, we adopted the Bayesian AVA inversion technique as described by Buland et al. (2008) and Kolbjørnsen et al. (2016Kolbjørnsen et al. ( , 2020. Figure 15 shows a schematic illustration of the inversion where, in the lower left, Markov chains (Mukerji et al. 2001;de Figueiredo et al. 2019) generate a set of possible vertical configurations of litho-fluid classes (LFCs) as illustrated by different colours. Then, in the middle, the elastic properties of the LFCs are extracted from well-log data or calculated from a rock physics model (Avseth et al. 2005) and, as illustrated to the right, subsequently convolved with a wavelet to generate synthetic seismic. The inversion tests many different LFC configurations (e.g. 754 facies combinations were employed in this study) by comparing their synthetic seismic response with the seismic prestack data trace by trace . If a good match is obtained between the synthetic response and the seismic data, a high probability value is assigned to that particular LFC configuration.
The inversion can incorporate geological information into the prior model from nearby well data, or from regional or local knowledge. For example, we assigned ±15 and ±10 ms horizon uncertainties to the Top Tuxen and Top Valhall horizons, respectively, in accordance with the confidence level from the seismic interpretation. As opposed to pre-stack elastic inversions using a low-frequency background model, the underlying Markov chain formulation yields a prior model with a more blocky rather than a smooth character . Each LFC 'block' (Fig. 15) corresponds to the seismic sampling, which is 4 ms in this study. Hence, the inversion is generally suitable for thin layers in blocky geology (Kolbjørnsen et al. 2020), given sufficient elastic contrast between the layers. However, in geological settings where the velocities are more compaction-driven, rock physics depth trends should ideally be incorporated into the prior model (Rimstad and Omre 2010).
Angle stacks of 5°-15°(near), 15°-25°(mid) and 25°-35°(far) were used as input to the inversion. The inversion results were optimized at the well locations by visual inspection where signal/ noise (S/N) scalars of 3, 4 and 2 were calibrated for the near, mid and far stacks, respectively. The S/N works as data weights for the various angle stacks, where near-mid stacks are typically considered the most reliable because near angles often contain residual multiples with larger bandwidth than primaries, and far angles may be influenced by normal moveout (NMO) stretch and differential attenuation. Figure 16a-e shows a well tie for the BO-2X well, with blocked AI and V P /V S ratio logs, and interpolated near, mid and far angle stacks for synthetic, real and residuals. To obtain a robust inversion result, it is important to obtain a consistent wavelet for the Valdemar Field for the seismic modelling (Fig. 15). A better well tie was obtained using a deterministic wavelet (White 1980;Walden and White 1998), which gave correlation coefficients of 82 and 56 within the Tuxen interval in BO-2X and BO-3X, respectively. The variability in correlation coefficients implies that the inversion results may be more unstable and uncertain in the Bo South area than in the Bo area. By visual comparison of the synthetic and real gathers, we see that there are some discrepancies and misalignments as also reflected by the residual gather, which can be due, for example, to remnant noise or time-depth relationships or wavelet inaccuracies. Figure 16f and g shows both the deterministic and statistical wavelets and the frequency spectra. Compared to the statistical wavelet, the deterministic wavelet suggests a small phase rotation to improve the well tie and contains somewhat lower frequencies at the high end of the spectrum. Smaller errors in wavelet phase and bandwidth will generally not significantly influence the inversion output but larger errors can produce apparent time shifts and inaccurate inversion models.  However, the wavelet estimation was generally challenging due to the limited length of well-log data (300 ms in BO-2X) required for good frequency fidelity and interfering events from strong contrasts in overburden (e.g. Base Chalk).
For each LFC, we define statistical mean values, standard deviations and correlations between pairs of the P-velocity, Svelocity and density properties extracted from different depth zones in the well logs. Based on these statistical parameters, we estimated Gaussian distributions that represent the rock physics likelihood models for each LFC, as illustrated in the middle of Figure 15 with a rock physics cross-plot with different LFC clusters. In this study, we are particularly interested in whether we can distinguish the upper and lower reservoir units of the Tuxen separated by the Munk Marl (MM) bed (Figs 5 and 6), as these are associated with different reservoir quality across the Valdemar Field. Therefore, we defined the following three different Tuxen reservoir scenarios: • Tuxen A represents a high-quality reservoir scenario with high porosities (c. 25-35%) and oil saturations (c. 70%), as observed from the upper Tuxen in BO-2X (Fig. 5).
We chose to define the Tuxen C LFC over the whole Tuxen interval because of little elastic separation across the MM bed in BO-3X and to keep a lower number of LFCs as input into the inversion to minimize the non-uniqueness of the problem and to boost the calculation speed. Figure 17 shows a AI v. V P /V S cross-plot with (a) the Gaussian distributions and (b) the Gaussian probability density functions (PDFs) for all the LFCs used in this study. The PDFs in Figure 17 represents the rock physics likelihood models used to calculate the posterior probabilities for each LFC (as illustrated in Fig. 15). Note that Tuxen A represents a slight decrease in AI and V P /V S from Tuxen B, whereas Tuxen A and B are separable from Tuxen C mainly by a lower AI. Sufficient elastic separation between the different LFCs is fundamental for the inversion to be able to distinguish seismically between the different LFCs as it is linked to the seismic response. Therefore, the more overlap we observe between two different LFCs in Figure 17, the more difficult it is for the inversion to distinguish seismically between them. Besides the different LFCs for the Tuxen interval, we also defined LFCs for the Sola and Valhall formations, the Munk Marl bed (only 1-2 m thick; see the MM markers in Figs 5 and 6), and two different sections of the lower part of the overlying chalk package. Figures 18 and 19 shows inversion QC plots along BO-2X and BO-3X, respectively, presented by the LFC (a) prior and (b) posterior probabilities and (c) the depth intervals used to define the LFCs in Figure 17, together with inverted and well log (d) AI and (e) V P /V S ratio. The prior model shows smooth boundary transitions as a result of assigning horizon uncertainties. Furthermore, as part of the prior geological constraints, a conditional facies transition was set such that Tuxen B and C, above Tuxen A, are impermissible to constrain the range of solutions. The primary outputs of the inversion are the posterior probabilities of the various LFC configurations. The secondary outputs used for inversion QC are the inverted AI and V P /V S properties calculated from a weighted mean based on the LFC posterior probability distributions. The normalized root mean square errors (NRMSEs) are estimated as a measure of the fit between the inverted and log values for the AI and V P /V S ratio. Within the Tuxen interval, the NRMSE for AI are low in both wells, whereas the NRMSE for the V P /V S ratio is slightly lower in BO-2X than in BO-3X. In the Upper Tuxen unit, we see high posterior probabilities of Tuxen A, as we would expect, which then decrease down through the Lower Tuxen unit. However, we do not see a similar high posterior probability within the Lower Tuxen unit for Tuxen B. This is because the inversion struggles with distinguishing Tuxen A from Tuxen B due to low elastic separation between the two LFCs (Fig. 17). Moreover, higher NRMSEs of the AI and V P /V S ratio are generally seen above Top Tuxen. This misfit may partly be due to local strong noise values in the seismic data around the Base Chalk interface at the well locations. Figure 20a shows the corresponding posterior probability for the high-quality reservoir scenario, Tuxen A, along the arbitrary line. Figure 20b shows in map view the maximum probabilities selected within a small window around the Tuxen horizon. The inversion suggests some higher probabilities in the northern part from the saddle point and towards BO-2X, and at the top of the structural high in the Bo area. Finally, Figure 21 shows the corresponding LFCs with the highest probabilities along (a) the arbitrary section and (b) the Top Tuxen horizon in map view. Tuxen A dominates the Bo area and conforms well with the structural high and elevation contours, whereas downflank areas are governed by Tuxen C, which is geologically reasonable.

Discussions
The reservoir quality distribution of the Lower Cretaceous unit in the Valdemar Field is poorly delineated, making further field development challenging. The aim of this study was therefore to unravel the seismic-reservoir relationship using modern seismic attributes and inversion tools. In this section, we discuss some associated uncertainties, shortcomings and possible improvements to our quantitative seismic interpretation results.
In general, the seismic data quality covering the Lower Cretaceous in the Valdemar Field and nearby fields is not optimal. In particular, multiples, converted waves and migration artefacts from the structurally complex overlying chalk may contaminate primaries from the Lower Cretaceous (Montazeri et al. 2018). Despite all processing efforts to improve the data, the regional chalk represents a limitation for the seismic image quality of sub-chalk targets (Klinkby et al. 2002). Figure 22 shows the  40 Hz component from spectral decomposition (Partyka et al. 1999) extracted from the Base Chalk horizon (Fig. 3) with a −40 ms time shift into the Upper Cretaceous unit. Some distinctive NE-SWpropagating faults are revealed in the centre and southwards to BO-3X, whereas a more chaotic dense network of faults and possible fractured (or damaged) zones propagates from the north and bypasses the BO-2X location. The velocity contrasts across fault offsets can lead to focusing or defocusing of seismic wavefronts, and erroneous amplitudes may be introduced by an inaccurate velocity model. This may explain why we see stronger amplitudes for the Top Tuxen event at the Bo structure that fade out downflank (Fig. 3) where the overburden seems to be less affected by faults and fractures. The posterior probabilities for the Tuxen A LFC are strongly correlated to these amplitudes (Fig. 20). Hence, it may be that the seismic amplitudes and inversion results are to some degree influenced by these overburden artefacts.
In the Bayesian AVA inversion, a single deterministic wavelet was applied to each angle stack due to the angle-dependent spectral balancing applied to the pre-stack data (Perez and Marfurt 2007). In the well-tie procedure (Fig. 16), we noted that small changes in wavelet parameters, such as the selection of a time extraction window and the number of traces around the wells to be included, had a big impact on the estimated wavelet and well-tie correlation. In BO-2X, the elastic logs cover the Upper Cretaceous unit and the upper part of the Lower Cretaceous, which is characterized by internal high-energy amplitudes from structural complexities and strong contrasts. When estimating a wavelet from a zone around the Tuxen reservoir, we ideally want to avoid including such highenergy amplitudes as it can bias the effective wavelet from our target interval. However, because the Tuxen reservoir is located close to the Base Chalk and the time extraction window needs to be large enough to capture low-frequency signals (Naeini et al. 2017), it was not possible to avoid including the Upper Cretaceous in the extraction window. Despite having longer elastic logs available in BO-3X, it was even more difficult to obtain a good well tie at this position, which may be related to lateral variations in coherent noise. As the Bayesian inversion results are generally very sensitive to changes in the input wavelet, any uncertainties from the well-tie procedure will eventually also affect the uncertainties of the inversion. One approach for investigating and fine-tuning wavelet   Figure 18. parameters (such as wavelet scaling) is to examine whether the variations in inverted AI and V P /V S are too oscillatory or too static compared to the logs (Figs 18 and 19) (Zhang et al. 2013). The oscillating (or dynamic) range in elastic properties can also be affected in their calculation by the weighted averaging of the LFC posterior probability distributions. The inversion in this study focused on calibrating input parameters for a best possible result within the Tuxen reservoir, where we obtained the lowest NRMSE.
Variations in seismic amplitudes and inversion results between Bo and Bo South can also be explained by differences in elastic contrasts for the Tuxen reservoir between the two areas. For example, the AI in BO-3X (Fig. 6) shows a much smaller contrast for the Top Tuxen than in BO-2X (Fig. 5). Moreover, the elastic contrasts between the upper and lower Tuxen units in BO-3X are very small, which makes it difficult to separate them seismically based on the Bayesian inversion approach (Kolbjørnsen et al. 2020). The Tuxen reservoir is buried around 155 m deeper in BO-3X compared to BO-2X, which is also reflected by lower porosities that are negatively correlated with AI (Fig. 8). Also, the porosity contrast to the overlying Sola Formation is much smaller in BO-3X. Therefore, the potential for obtaining unambiguous and wellresolved seismic reservoir characterization of the Tuxen reservoir in the Bo South area will generally be lower compared with the Bo area.
As input to the Bayesian inversion, we used local Gaussian rock physics priors estimated from the BO-2X and BO-3X AI and V P /V S logs for different LFCs within specific depth intervals. However, the selection of LFCs may not be sufficient to represent all possible geological and fluid scenarios present in the Valdemar Field. An alternative could be to use the calibrated iso-frame model presented in Appendix A and perform stochastic simulations for different 'what-if' scenarios (Avseth et al. 2005). However, the main goal of this study was to see if we could segregate various units of the Tuxen reservoir with different reservoir-quality scenarios from well-log data and to keep the inversion simple. Furthermore, a rock physics puzzle that requires further research is to investigate how the impact of clay variations and fractures affect the seismic properties and possible anisotropy effects in the Lower Cretaceous.

Conclusions
A quantitative seismic interpretation study is presented for the Lower Cretaceous Tuxen reservoir in the Bo and Bo South areas located in the Valdemar Field, Danish North Sea. The workflow consists of seismic and rock physics forward modelling, seismic AVA attribute analysis, a relative AI inversion, and an absolute Bayesian AVA inversion for litho-fluid classification. The forward modelling showed that the seismic signal is more sensitive to porosity than water-saturation changes, causing AI contrasts mainly sensed on at the near AVA amplitudes (or the AVA intercept). The seismic AVA intercept revealed some variations for the Tuxen reservoir around the Bo structure associated with high porosities. Applying a CI operator to the AVA intercept apparently improved the resolution at the reservoir level and provided a possible porosity indication map for the Tuxen interval. Finally, the Bayesian AVA inversion predicted higher probabilities around the Bo area for the litho-fluid class representing the upper Tuxen unit associated with high porosity and low water saturation. The results generally agreed with well-log observations at the reservoir level and the anomalies conform to the Bo structure. In general, this type of Bayesian inversion can be useful for boosting reservoir resolution, incorporating geological information and classifying the seismic data into litho-fluid classes in a probabilistic manner. We can expect the inversion to provide good results if the different litho-fluid classes exhibit separable elastic properties, the seismic data have been properly conditioned to remove interfering noise, and representable wavelets and prior models for the target interval can be estimated. Despite not all these circumstances being satisfactory fulfilled in this study, the results can offer useful insights for further field development and reservoir management under the conditions given by the available dataset and geological complexities.
Elastically isotropic rocks were assumed in this study, such that seismic properties could be derived from the bulk modulus (K ), shear modulus (μ) and density (ρ). These were then used to calculate the acoustic impedance (AI), and the ratio between P-wave velocity (V P ) and S-wave velocity (V S ) (i.e. V P /V S ) via: For the modelling of the effective K and μ properties of the Tuxen reservoir, the iso-frame (IF) rock physics model was used (Fabricius et al. 2005(Fabricius et al. , 2007. Mineralogical and petrophysical analysis from core data from the Tuxen reservoir were used to assess the relative volume fractions of solid components of the rock (Table A1), specified by large calcite (c l : microfossils and carbonate cement) and large silicates (quartz and feldspar clasts (q l ), and clay clasts (k l )), as well as the amount of small calcite (c s ) and silicates (quartz (q s ) and clay (k s )) in suspension. The IF model assumes an effective frame (composed of the large calcite and silicate components) with spherical pores. The pores are filled with an effective pore fluid with small calcite and silicate particles in suspension, where their contribution to the effective stiffening of the frame is specified by the IF value (IF = 1 implies that all suspended materials have entered the effective frame). An effective pore-fluid mixture of water and light oil according to the model of Wood (1955) was assumed. The bulk modulus for the effective suspension for a given porosity (f) and water saturation (S w ) was calculated by (Fabricius et al. 2007):